Chapter 5
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.109375beta.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.3330165 and the mean number of iteration is 2.10559hist(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)
Example 6:
Generate \(X\) from \(N(0,1)\) by using the AR method with \(g(x) = e^{-x}, \ x>0\). Note that the support of the normal distribution is not the same as the exponential. Therefore, we let \(Y = \lvert X \rvert\).
\[ f_Y(y) = \frac{2}{\sqrt{2\pi}}e^{-x^2/2}, \quad y>0. \]
Then
\[ h(x) = \sqrt{\frac{2}{\pi}}e^{x-x^2/2} \]
\[ h^{\prime}(x) = \sqrt{\frac{2}{\pi}}e^{x-x^2/2}(1-x) \overset{set}{=} 0 \]
When \(x=1\), \(h(x)\) has a maximum \(c = \sqrt{2/\pi}e^{1-1/2} = \sqrt{2e/\pi}\).
std.norm = function(n){
X = rep(NA,n)
iter = rep(NA,n)
f = function(x){sqrt(2/pi) * exp(-x^2/2)}
g = function(x){exp(-x)}
h = function(x){ f(x)/g(x) }
max.h = optim(1, h, method = "SANN", control = list(fnscale = -1))
c = max.h$value
for(j in 1:n){
Y = -log(runif(1)) ### simulate exp(1)
cri = h(Y)/c
i = 1
while( runif(1) > cri ){
Y = -log(runif(1)) ### simulate exp(1)
cri = h(Y)/c
i = i+1
}
X[j] = ifelse(runif(1) < 0.5, Y, -Y)
iter[j] = i
}
return(list("X" = X, "iter" = iter, "c" = c))
}
res = std.norm(100000)
hist(res$X, breaks = 50, probability = T)
lines(seq(-2.5,2.5,0.01), dnorm(seq(-2.5,2.5,0.01)), col = "blue", lwd = 2)
cat("The maximum of h is", res$c, "and the average iteration time is", mean(res$iter))## The maximum of h is 1.315489 and the average iteration time is 1.31047Next, consider \(N(\mu, \sigma^2)\) distribution.
norm.fun = function(n, mu, var){
Z = std.norm(n)$X
return(mu + Z * var)
}
res = norm.fun(100000, 5, 3)
hist(res, breaks = 50, probability = T)
lines(seq(-5,15,0.01), dnorm(seq(-5,15,0.01), 5, 3), col = "blue", lwd = 2)
Example 7:
Let \(X\) from the distribution of \(Gamma(2,1)\) conditional on \(X \ge 5\).
\[ f(x) = P(X = x \mid x \ge 5) = \frac{P(X=x, X \ge 5)}{P(X \ge 5)} = \frac{xe^{-x}}{\int_5^\infty xe^{-x} dx} = \frac{xe^{-x}}{ -xe^{-x}-e^{-x} \mid_5^\infty } = \frac{xe^{-x}}{6e^{-5}}, \ x \ge 5 \]
Let \(Y_1 \sim exp(2) \mid Y_1 \ge 5\), \(Y_2 \sim exp(1) \mid Y_2 \ge 5\), \(Y_3 \sim exp(3) \mid Y_3 \ge 5\). That is,
\[ g_1(x) = \frac{\frac{1}{2}e^{-x/2}}{e^{-5/2}}, \quad x \ge 5 \]
\[ g_2(x) = \frac{e^{-x}}{e^{-5}}, \quad x \ge 5 \]
\[ g_3(x) = \frac{\frac{1}{3}e^{-x/3}}{e^{-5/3}}, \quad x \ge 5 \]
Then,
\[ c_1 = \max_{x\ge5}\frac{f(x)}{g_1(x)} = \max_{x\ge5}\frac{xe^{-x/2}}{3e^{-5/2}} = 5/3 \]
\[ c_2 = \max_{x\ge5}\frac{f(x)}{g_2(x)} = \max_{x\ge5}\frac{x}{6} = \infty \]
\[ c_3 = \max_{x\ge5}\frac{f(x)}{g_3(x)} = \max_{x\ge5}\frac{xe^{-2x/3}}{2e^{-10/3}} = 5/2 \]
ex7.fun = function(n){
X = rep(NA,n)
iter = rep(NA,n)
f = function(x){return( x*exp(-x)/(6*exp(-5)) )}
g = function(x){return( dexp(x,1/2)/exp(-5/2) )}
h = function(x){return( f(x)/g(x) )}
max.h = optim(5, h, lower = 5, upper = 10, method = "Brent",
control = list(fnscale = -1))
c = max.h$value
for(j in 1:n){
Y = 5-2*log(runif(1))
cri = h(Y)/c
i = 1
U = runif(1)
while(U > cri){
Y = 5-2*log(runif(1))
cri = h(Y)/c
i = i+1
U = runif(1)
}
X[j] = Y
iter[j] = i
}
return(list("X"=X, "iter"=iter, "c"=c))
}
res = ex7.fun(100000)
hist(res$X, breaks = 50, probability = T, xlim = c(5,14))
x = seq(5,14,0.01)
lines(x, x*exp(-x)/(6*exp(-5)), col = "blue", lwd = 2)