# Math Help - rejection sampler

1. ## rejection sampler

Q) Explain the principle behind the rejection sampler, and specifically outline how you would use this approach to generate Gamma variables with mean 10,000 and variance 1,000,000. Implement an algorithm to do this and provide appropriate graphical and statistical results to demonstrate that it is working correctly. Briefly comment on the efficiency of this method when
compared with the previous method.Provide working code for this rejection sampler in codes. R

i have the following code but every single value is getting rejected. Does anybody know why and help me fix it?
Thanks

maxd <- dgamma(99/0.01, 100,0.01) # the maximum density
x <- rexp(1000,0.01) # as exponention has the same support randomly select from it
u <- runif(1000) # as all samplers have the uniform distribution
accept <- ifelse(u < dgamma(x, 100,0.01)/maxd * dunif(u), 1, 0) # the ratio
sample <- x[accept == 1] # adding up all the accepted values but there are none
ks.test(sample, "pgamma", 100, 0.01) # Kolmogorov-Smirnov test which doesnt work as no values to do test on

2. Hello,

I think you should read something about the rejection method... You need to find another pdf, that you can simulate easily...

Here is what I did for simulating the Gamma distribution (a,1), with a given g : $g(x,a)=\frac{ae}{a+e}(x^{a-1}\bold{1}_{x\in (0,1)}+e^{-x}\bold{1}_{x\geq 1})$

and with $c=\frac{a+e}{ae\Gamma(a)}$

we simulate rv's with density g by the inversion method. (the generalized inverse function of the cdf of g is defined as F here)

Code:
f <- function(x,a){
y=1/gamma(a)*exp(-x)*x^(a-1)*(x>0)
y
}
e=exp(1)

g <- function(x,a){
y=(a*e)/(a+e)*(x^(a-1)*((x>0)&(x<1))+exp(-x)*(x>1))
y
}
g(-1,6)

F <- function(z,a){
y=((a+e)/e*z)^(1/a)*(z<e/(a+e))-log((a+e)/(a*e)*(1-z))*(z>=e/(a+e))
y
}

gamaux<-function(n,a){
x <- rep(0,n)
c=(a+e)/(a*e*gamma(a))
count<- 1
while (count< n+1){
v <- runif(1)
u=runif(1)
y=F(u,a)
if(v<=f(y,a)/(c*g(y,a))){
x[count] <- y
count <- count+ 1
}
}
return(x)
}
which will give n rv's following a gamma(a,1)

I hope this helps.