1 The Inverse Transform Algorithm

Inverse transformation method:

Let a continuous random variable have distribution function $F$. A general method is based on the following proposition.

Proposition

Let $U$ be a uniform $(0,1)$ random variable. For any continuous distribution function $F$$, the random variable $X$ which is defined by

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

has distribution $F$.

proof:

\[ F_X(x) = P(X \le x) = P(F^{-1}(U) \le x) = P(U \le F(x)) = F(x) \]

Example 1:

$F(x)=x^n, \ 0<x<1$. Then $x = F^{-1}(u) = u^{1/n}, \ 0<u<1$.

Example 2:

If $X$ is a exponential distribution with rate 1, then the distribution function is

\[ F(x) = 1-e^{-x}, \ x>0. \]

Then,

\[ x = F^{-1}(u) = -\log(1-u), \ 0<u<1. \]

$\because U \sim unif(0,1) \ \text{and} \ 1-U \sim unif(0,1)$

$\therefore X=-\log(U)$

Example 3:

If $X$ is a exponential distribution with rate $\lambda$, then the distribution function is

\[ F(x) = 1-e^{-\lambda x}, \ x>0 \]

Then

\[ x = F^{-1}(u) = \frac{-1}{\lambda}\log(u) \]

and

\[ X = \frac{-1}{\lambda} \log(1-U) \ \text{or} \ X = \frac{-1}{\lambda} \log(U), \]

since both $U$ and $1-U$ are uniform $(0,1)$.

Therefore,

\[ Y=\frac{1}{\lambda} X, \ \text{where} \ Y \sim Exp(\lambda) \ \text{and} \ X \sim Exp(1) \]

Example 4:

If $X$ is a gamma distribution $(n,\lambda)$. then the distribution function is

\[ F(x) = \int_0^x \frac{ \lambda e^{-\lambda y} (\lambda x)^{n-1} }{(n-1)!} dy, \ x>0 \]

Let $X_1, \dots, X_n \overset{i.i.d.}{\sim} Exp(\lambda)$. Then $Z = \sum_{i=1}^n X_i \sim \Gamma(n, \lambda)$. Therefore,

\[ X_i = \frac{-1}{\lambda}\log(1-U_i), \ i=1,\dots,n \quad \Rightarrow \quad Z = \sum_{i=1}^n \frac{-1}{\lambda}\log(1-U_i) \]

2 The Rejection Method

Suppose that we want to generate an $X$$ random variable from density function $f(x)$. Suppose that we can easily generate a random variable, $Y$ , with density function $g(y)$. Then, we accept the generated value with a probability proportion to $f(y)/g(y)$,

\[ \frac{f(y)}{g(y)} \le c, \ \forall y \]

Example 5:

Use $\boldsymbol{U(0,1)}$ to draw $\boldsymbol{Beta(2,4)}$.

$X$ is a random variable from $Beta(2,4)$, and the pdf is

\[ f(x) = 20x(1-x)^3, \ 0<x<1 \]

Let us consider the rejection method with

\[ g(y) = 1, \ 0<y<1 \]

Let

\[ h(x) = \frac{f(x)}{g(x)} = 20x(1-x)^3 \le c \]

Then

\[ \frac{dh}{dx} = 20[ (1-x)^3 - 3x(1-x)^2 ] \overset{\text{set}}{=} 0 \]

$\Rightarrow \ x = 1/4 \mbox{ or } 1 \ \Rightarrow \ h(1/4)=135/64 \mbox{ and } h(1)=0$. Therefore, $c = 135/64$.

Alternatively, we can use optim() function to find maximum.

f.g = function(x){
  dbeta(x,2,4) / 1
}
x = seq(0, 1, 0.01)
plot(x, f.g(x), type = "l")

Next, we choose 0.2 as the initial value for finding an optimum.

det.c = optim(0.2, f.g, lower = 0, upper = 1, method = "Brent",
              control = list(fnscale = -1)) ### maximize
cat("c is", det.c$value) ### the location of the optimum

## c is 2.109375

beta.unif = function(n, a, b){
  X = rep(NA, n)
  iter = rep(NA, n)
  
  h = function(x){dbeta(x,a,b)}
  det.c = optim(0.2, h, lower = 0, upper = 1, method = "Brent",
              control = list(fnscale = -1))
  par = det.c$par
  c = det.c$value
  
  for(j in 1:n){
    Y = runif(1,0,1)
    cri = h(Y)/c
    i = 1
    U = runif(1,0,1)
    while(U > cri){
      Y = runif(1,0,1)
      cri = h(Y)/c
      i = i+1
      U = runif(1,0,1)
    }
    X[j] = Y
    iter[j] = i
  }
  return(list("X" = X,
              "iter" = iter,
              "par" = par,
              "c" = c))
}

res = beta.unif(100000,2,4)
cat("The sample mean of X is", mean(res$X), 
    "and the mean number of iteration is", mean(res$iter))

## The sample mean of X is 0.3330219 and the mean number of iteration is 2.10084

hist(res$X, breaks = 50, probability = T)
lines(seq(0,1,0.01), dbeta(seq(0,1,0.01),2,4), col = 'blue', lwd = 1.2)

Use another distribution to draw $\boldsymbol{Beta(6,4)}$

The target pdf is

\[ f(x) = \frac{\Gamma(10)}{\Gamma(6)\Gamma(4)}x^5(1-x)^3, \quad 0<x<1. \]

From example 1, we can easily draw samples from a distribution with $g(x)=5x^4, \ 0<x<1$. By the rejection method, let

\[ h(x) = f(x)/g(x) = \frac{9!}{5!3! \times 5} x(1-x)^3 \le c ,\quad 0<x<1 \]

\[ h^{\prime}(x) = \frac{7\times8\times9}{5}[(1-x)^3-3x(1-x)^2] \overset{set}{=}0 \]

$\Rightarrow \ x = 1/4 \mbox{ and } c = 1701/160$

f = function(x){
  dbeta(x,6,4)
}
g = function(x){
  5 * x^4
}
h = function(x){
  f(x)/g(x)
}

plot(seq(0,1,0.01), h(seq(0,1,0.01)), type = 'l',
     xlab = "X", ylab = "h(x)")

beta = function(n, a, b){
  X = rep(NA,n)
  iter = rep(NA,n)
  
  f = function(x){dbeta(x,a,b)}
  g = function(x){(a-1)*x^{a-2}}
  h = function(x){f(x)/g(x)}
  
  det.c = optim(0.2, h, lower = 0, upper = 1, method = "Brent",
                control = list(fnscale = -1))
  maximizer = det.c$par
  c = det.c$value
  
  for(j in 1:n){
    Y = runif(1)^(1/(a-1))
    cri = h(Y)/c
    U = runif(1)
    i = 1
    while (U > cri) {
      Y = runif(1)^(1/(a-1))
      cri = h(Y)/c
      U = runif(1)
      i = i+1
    }
    X[j] = Y
    iter[j] = i
  }
  return(list("X" = X,
              "iter" = iter,
              "maximizer" = maximizer,
              "c" = c))
}

res = beta(100000, 6, 4)
hist(res$X, breaks = 50, probability = T)
lines(seq(0,1,0.01), dbeta(seq(0,1,0.01), 6, 4), col = "blue", lwd = 2)