Results 1 to 2 of 2

Math Help - rejection sampler

  1. #1
    Junior Member
    Joined
    Feb 2009
    Posts
    41

    Question 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
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Moo
    Moo is offline
    A Cute Angle Moo's Avatar
    Joined
    Mar 2008
    From
    P(I'm here)=1/3, P(I'm there)=t+1/3
    Posts
    5,618
    Thanks
    6
    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.
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Rejection Sampling
    Posted in the Advanced Statistics Forum
    Replies: 1
    Last Post: October 20th 2011, 06:35 AM
  2. Fast convergence for the random site Gibbs sampler
    Posted in the Advanced Statistics Forum
    Replies: 0
    Last Post: June 1st 2010, 02:04 PM
  3. Evaluating my sampler
    Posted in the Advanced Statistics Forum
    Replies: 0
    Last Post: May 4th 2010, 06:52 AM
  4. Rejection Region
    Posted in the Advanced Statistics Forum
    Replies: 0
    Last Post: November 30th 2009, 12:32 PM
  5. Showing the rejection region
    Posted in the Advanced Statistics Forum
    Replies: 0
    Last Post: November 16th 2009, 06:16 AM

Search Tags


/mathhelpforum @mathhelpforum