# Note: this is only a toy example demonstrating how to program rejection sampling
# this is not a good sampling scheme for gamma distribution.
# Do not use it for serious applications that demand high efficiency.
#log of a function which is always above the Gamma density function
# alpha must be > 2
log_g_gamma <- function(x, alpha)
{
(alpha-1) * (log(alpha-1) - 1) - log( 1 + (x-(alpha-1))^2 / (2*alpha-1) )
}
## look at the approximating functionn
xvec <- seq (0, 3, by = 0.0001)
alpha <- 1.5
log_g_gamma_val <- log_g_gamma(xvec, alpha = alpha)
log_gamma_val <- dgamma (xvec, shape = alpha, log =T)
ylim <- range (log_g_gamma_val, log_gamma_val, finite = T)
plot (xvec, log_g_gamma_val, col = "black", type = "l", ylim = ylim )
points (xvec,log_gamma_val, col = "red",type = "l")
## look at the approximating functionn
xvec <- seq (0, 3, by = 0.0001)
alpha <- 2.1
log_g_gamma_val <- log_g_gamma(xvec, alpha = alpha)
log_gamma_val <- dgamma (xvec, shape = alpha, log =T)
ylim <- range (log_g_gamma_val, log_gamma_val, finite = T)
plot (xvec, log_g_gamma_val, col = "black", type = "l", ylim = ylim )
points (xvec,log_gamma_val, col = "red",type = "l")
## look at the approximating functionn
xvec <- seq (0, 3, by = 0.0001)
alpha <- 4.5
log_g_gamma_val <- log_g_gamma(xvec, alpha = alpha)
log_gamma_val <- dgamma (xvec, shape = alpha, log =T)
ylim <- range (log_g_gamma_val, log_gamma_val, finite = T)
plot (xvec, log_g_gamma_val, col = "black", type = "l", ylim = ylim )
points (xvec,log_gamma_val, col = "red",type = "l")
#sampling from Gamma distribution with rejection sampling
sample_gamma_rej <- function(n,alpha)
{ sample_gamma <- rep(0,n)
no.draw <- 0
for(i in 1:n)
{ rejected <- TRUE
while(rejected)
{ sample_gamma[i] <- rcauchy(1) * sqrt(2*alpha-1) + (alpha -1)
no.draw <- no.draw + 1
U <- runif(1)
rejected <- (log(U) > dgamma(sample_gamma[i],shape=alpha,log=TRUE) -
log_g_gamma(sample_gamma[i],alpha) )
}
}
attr(sample_gamma, "accept.rate") <- n/no.draw
sample_gamma
}
alpha <- 2.1
gammarn <- sample_gamma_rej (2000,alpha); attr (gammarn, "accept.rate")
## [1] 0.4878049
hist (gammarn)
qqplot(gammarn, rgamma (2000, alpha))
abline (a = 0, b=1)
alpha <- 2.5
gammarn <- sample_gamma_rej (2000,alpha); attr (gammarn, "accept.rate")
## [1] 0.3800114
hist (gammarn)
qqplot(gammarn, rgamma (2000, alpha))
abline (a = 0, b=1)
alpha <- 4.5
gammarn <- sample_gamma_rej (2000,alpha); attr (gammarn, "accept.rate")
## [1] 0.04686804
hist (gammarn)
qqplot(gammarn, rgamma (2000, alpha))
abline (a = 0, b=1)
We see that when alpha is larger than 2, the overall acceptance rate is very low
If you requires highly efficient gamma random numbers generators, read a paper in Computational Statistics and Data Analysis (2007) (http://www.sciencedirect.com/science/article/pii/S0167947306003616)
library (ars)
## a direct rejection sampling for truncated normal
sample_tnorm_drs <- function (n, lb = -Inf, ub = Inf)
{
x <- rep (0, n)
for (i in 1:n)
{
rej <- TRUE
while (rej)
{
x[i] <- rnorm (1)
if (x[i] >= lb & x[i] <= ub) rej <- FALSE
}
}
x
}
## sample from truncated normal using ars package
sample_tnorm_ars <- function (n, lb, ub)
{
logf <- function (x) dnorm (x, log = TRUE) ## define log density
fprima <- function (x) -x ## define derivative of log density
ars (n, f = logf, fprima = fprima,
x = c(lb, (lb + ub )/2, ub), # starting points
lb = TRUE, ub = TRUE, xlb = lb, xub = ub) # boundary of log density
}
n <- 1000
system.time (
rn_tnorm_ars <- sample_tnorm_ars (n, -5, -4)
)
## user system elapsed
## 0.012 0.001 0.013
system.time(
rn_tnorm_rej <- sample_tnorm_drs (n, -5, -4)
)
## user system elapsed
## 46.525 8.763 56.840
par (mfrow = c(1,3))
hist (rn_tnorm_ars, main = "Adaptive Rejection Sampling")
hist (rn_tnorm_rej, main = "Naive Rejection")
qqplot (rn_tnorm_ars, rn_tnorm_rej)
abline (a = 0, b = 1)
# Draw truncated normal sample on the far tail
n <- 1000
system.time (
rn_tnorm_ars2 <- sample_tnorm_ars (n, -50, -40)
)
## user system elapsed
## 0.023 0.001 0.023
system.time (
rn_tnorm_ars3 <- sample_tnorm_ars (n, 100, 110)
)
## user system elapsed
## 0.007 0.000 0.007
par (mfrow = c(1,2))
hist (rn_tnorm_ars2, main = "Adaptive Rejection Sampling")
hist (rn_tnorm_ars3, main = "Adaptive Rejection Sampling")