Random Number Generation Using Inverse Transformations

Generating random numbers that follow a specific probability distribution is a fundamental problem in statistical computing and simulation. While generating uniform random numbers is straightforward, obtaining random numbers from non-uniform distributions requires specialized techniques. Among the various methods available, the inversion method is one of the simplest and most widely used approaches.

The inversion method is particularly useful when the cumulative distribution function (CDF) of the desired random variable has a closed-form expression that can be easily inverted. The key idea behind this method is that if \(X\) is a continuous random variable with CDF \(F(\cdot)\), then the transformation \[X=F^{−1}(U)\] where \(U\sim U(0,1)\) (a uniform random variable between 0 and 1), generates random values \(X\) that follow the desired distribution. This technique is based on the fundamental principle that the CDF of any continuous random variable transforms a uniform random variable into a variable with the desired distribution.

Theorem

For any continuous random variable \(X\), its cumulative distribution function (CDF), follows a uniform distribution between 0 and 1.

Proof

Given that \(X\) is a continuous random variable with distribution function \(F(\cdot)\). Let \(Y=F(x), X\in \mathbb{R}\) We can see that the transformation \(Y= F(x), X\in \mathbb{R}\) is a one-to-one transformation. Hence we can use the Jacobian Tranformation to find the density function of the random variable \(Y\).

\[\begin{eqnarray} \frac{dY}{dX}&=& \frac{dF(x)}{dx}=f_x(x)\\f_Y(y)&=&f_X(x \text{ in terms of }y) \left|\frac{dx}{dy}\right|=f_x(x)\times \frac{1}{f_X(x)}=1,; 0\leq x\leq 1 \end{eqnarray}\]

Random numbers for the distribution having a compact form for the distribution function can be generated using the following steps

  1. Generate the uniform random variables

  2. Let the random variable \(X\) have a continuous and increasing distribution function \(F(.)\). By the above result, distribution function \(F(.) \sim U(0,1)\). Let \(F^{-1}\) denote the inverse of \(F\). Then find \(F^{-1}\) for every random number generated in step 1.

  3. We can use the same method in the case of discrete random variables, but the inverse being taken using the formula \[\begin{equation} F^{-1}(u)=\min \{x:F(x)>u\} \end{equation}\]

Generate random numbers from rectangular distribution range [10,20]

Let \(X\) be a random variable following rectangular distribution with parameters \(a\) and \(b\), ie., \(X \sim U(a,b)\). Then the distribution function is given by \[ F(x)=\begin{cases} 0 & x \leq a\\\frac{x-a}{b-a}, &a \leq x \leq b\\ 0 & x\geq b \end{cases}\]

Then the inverse between the range \(a\leq x\leq b\) is given by \[\begin{eqnarray*} u=F(x) &\implies& u= \frac{x-a}{b-a} \implies x-a= u(b-a) \implies x=a+u(b-a) \end{eqnarray*}\]

Now we shall consider the R function to generate abitrary number \(n\) of random numbers from rectangular distribution with parameters \(a\) and \(b\).

RectaRN=function(a,b,n){
  if (a >=b){
    print("Invalid Arguments, a should be less than b")
    }else if (n%%1!=0){
    print("Give the correct number of random numbers to be generated") 
    }else {
    u=runif(n)
    rand=a+(b-a)*u
    return(rand)
  }
}

Now we shall use this function

RectaRN(5,2,8)
## [1] "Invalid Arguments, a should be less than b"

An error occurs because the value of \(a\) is greater than \(b\). Now, let’s consider the case where the specified number of random numbers to be generated is a non-integer.

RectaRN(2,8,5.4)
## [1] "Give the correct number of random numbers to be generated"

Now we will generate 10 rectangular random numbers between 5 and 8.

RectaRN(5,8,10)
##  [1] 7.094923 7.524041 6.946635 7.156392 6.406077 6.793078 5.441658 5.462514
##  [9] 7.721091 5.060565

The same result shall be obatined by calling the function by arguments, specifying the arguments by name, it does not be in same order

RectaRN(b=8,a=5,10)
##  [1] 5.146661 5.567034 6.643321 6.024222 6.745984 7.156222 6.433045 5.920168
##  [9] 5.205208 7.736361

Next, we will plot a histogram of 1,000 uniformly distributed random numbers between 5 and 15, using 50 bins. Additionally, we will overlay the theoretical density on the plot for comparison.

a=5
b=15
set.seed(6666776)
rns=RectaRN(a,b,1000)
hist(rns,breaks = 50,probability=TRUE,
     col='lightblue',border='red',main=paste("Histogram of Reactangular distribution over(",a,",",b,")"),xlab="values",ylab="Density")
curve((x >= a & x <= b) * (1 / (b - a)),from = a,to=b,add = TRUE,col='green',lwd=2)
# Overlay the theoretical density
#curve(dunif(x, min = 5, max = 15), from = 5, to = 15, 
  #   col = "red", lwd = 2, add = TRUE)

# Add a legend
legend("topright", legend = c("Theoretical Density"), 
       col = "green", lwd = 2)

Now we will test whether the samples generated follows rectangular distribution using Kolmogrov-Smrinov test.

c=ks.test(rns,'punif',min=a,max=b)
c
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  rns
## D = 0.02111, p-value = 0.7644
## alternative hypothesis: two-sided

Since p-value, 0.7643598 is not small, we cannot reject the nul hypothesis that the data follows rectangular distribution over the range 5 and 15. Hence the data may be following rectangular distribution over the range 5 and 15.

Generating expenential random numbers from standard uniform random numbers.

The density function and the distribution function of the exponential random variable with rate \(\lambda\) or mean \(\dfrac{1}{\lambda}\) is given by \[f(x)=\lambda e^{-\lambda x}, x\geq 0 \text{ and } F(x)=1-e^{-\lambda x} \].

To generate an exponential random number using the inverse transformation method, follow these steps:

  1. Generate a uniform random variable \(U \sim U(0,1)\).

  2. Set \(X\) as the inverse of the CDF:

    \[ X = F^{-1}(U) \]

    Since \(F(x) = 1 - e^{-\lambda x}\), solving for \(x\):

    \[ U = 1 - e^{-\lambda x} \]

    Rearranging,

    \[ e^{-\lambda x} = 1 - U \]

    Taking the natural logarithm on both sides,

    \[ -\lambda x = \ln(1 - U) \]

    Finally,

    \[ X = -\frac{\ln(1 - U)}{\lambda} \]

    Since \(U \sim U(0,1)\), and \(1 - U\) also follows the same distribution, we often simplify this to:

    \[ X = -\frac{\ln U}{\lambda} \]

Thus, to generate an exponential random variable with rate \(\lambda\), use:

\[ X = -\frac{\ln U}{\lambda}, \quad U \sim U(0,1) \]

This method ensures that the generated random numbers follow the desired exponential distribution.

As an illustration, we will generate expoenential random numbers with rate \[\lambda=2\] or mean \[\frac{1}{2}\] in R.

#rate <- as.numeric(readline(prompt="Enter a rate of exponential random distribution: "))
#As we want the output we will assume value for rate.
#Let n be the number of random numbers to be generated

exprn=function(rate,n){
  if (rate<=0){
    print("Error!!!!! For exponential distribution, rate should be positive")
  } else {
    u=runif(n)
    exprand=(-1)*log(u)/rate
    return(exprand)
  }
}

Now we shall generate the exponential random numbers for specific values of rate.

exprn(-5,56)
## [1] "Error!!!!! For exponential distribution, rate should be positive"
rate=0.5
exp1=exprn(rate,100)
expdens=function(rate,x){
  return(rate*exp((-1)*rate*x))
}
hist(exp1,breaks = 50, col = 'lightblue',border = 'black',probability = TRUE,main = paste("Histogram of Exponential Random numbers with Rate",rate),xlab = "Values",ylab = "Density")

curve(expdens(rate,x),from=0, to=max(exp1),add=TRUE, col='green',lwd=2)

Now we shall verify the goodness of fit using Kolmogrov- Smronov Test.

# The ML estimator for the rate is used to estimate the parameter 
rateHat=1/mean(exp1)
rateHat
## [1] 0.4177688
ks.test(exp1,'pexp',rate=rateHat)
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  exp1
## D = 0.099111, p-value = 0.2797
## alternative hypothesis: two-sided

Since the p-value is large, we can not reject the null hypothesis that the data follows normal distribution.

Generating Cauchy random numbers from standard uniform random numbers.

The density function of Cauchy distribution with parameter \(\sigma\) is given by

\[f(x)=\dfrac{1}{\pi} \dfrac{\sigma}{x^2+\sigma^2},-\infty<x<\infty \]

The cumulative distribution function (CDF) is:

\[ F(x) = \frac{1}{\pi} \tan^{-1} \left( \frac{x }{\sigma} \right) + \frac{1}{2} \]

Now we shall consider the inverse transformation. \[\begin{aligned} U = F(X) &\implies U= \frac{1}{\pi} \tan^{-1} \left( \frac{X }{\sigma} \right) + \frac{1}{2} \\ & \implies U-\frac{1}{2}= \frac{1}{\pi} \tan^{-1} \left( \frac{X }{\sigma} \right) \\ &\implies \pi \left(U-\frac{1}{2}\right) = \tan^{-1} \left( \frac{X }{\sigma} \right) \\ &\implies \tan \left[ \pi \left(U-\frac{1}{2}\right)\right] = \frac{X }{\sigma} \\ &\implies X = \sigma \tan \left[ \pi \left(U-\frac{1}{2}\right)\right] \end{aligned}\]

Reference

  1. Ross, S. M. (2022). Simulation (6th ed.). Academic Press. https://doi.org/10.1016/C2019-0-02165-0

  2. Devroye, L. (1986). Non-uniform random variate generation. Springer.

  3. Gentle, J. E. (2003). Random number generation and Monte Carlo methods (2nd ed.). Springer