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 : $\displaystyle 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 $\displaystyle 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.