Tell me more ×
Cross Validated is a question and answer site for statisticians, data analysts, data miners and data visualization experts. It's 100% free, no registration required.

I'm working on some practice test problems, and one of them says to design a rejection sampling algorithm to produce draws from a unit exponential using draws from a Gamma(2,1).

I don't understand how this is possible, because I am under the impression that the "envelope function" g(x) needs to be scalable in such a manner that for some constant $M$, $Mg(x)\geq f(x)\; \forall x$.

I can't see any way to do this, as the Gamma(2,1) is going to have little mass around 0, while the exponential function has most of its mass around 0. What kind of transformation do I need to do to the Gamma function to allow it to function as an envelope?

Using R, I tried flipping it to make it an inverse gamma, but that won't adequately capture the probability mass close to 0 and a K-S test confirmed that the points I generated did not arise from a unit exponential.

edit: I will include my code in which I tried to use an inverse-gamma(2,1) as an envelope:

x <- c()
for(i in 1:100000)
{
  g <- runif(1, 0, 1) 
  h <- rigamma(1, 2, 1)
  M <- densigamma(h, 2, 1)
  crit <- dexp(h, 1)/M
  if(g < crit)
    x[i] = h
}
share|improve this question
Wikipedia provides (without explanation) an acceptance-rejection sampling algorithm for $\Gamma(\delta), 0 \lt \delta \lt 1$. Using it you could generate (say) two $\Gamma(1/2)$ values and add them. – whuber May 4 '11 at 21:46
Thanks for the link, but I am supposed to use the $\Gamma(2,1)$ as the envelope function. – Will May 4 '11 at 22:55
One small thing, rather than constantly resizing your vector, start off with: x = numeric(100000) This will be much faster. – csgillespie May 5 '11 at 21:11

1 Answer

up vote 4 down vote accepted

Try a location shift on the Gamma(2,1)

EDIT: Illustration

share|improve this answer
1  
+1 You got it in just as I was writing the same thing... :-). The expected number of $\Gamma(2)$ values needed per Exponential draw is about 3.407 with the optimal shift (of $\sqrt{1/2}$). – whuber May 4 '11 at 23:18
Thanks, but I am not sure if I understand how to implement a location shift. If I understand correctly, I just adjust the Gamma function by adding a constant, but I am not exactly sure where to put it in my code. – Will May 4 '11 at 23:36
@Will Plot $x \to pdf(\Gamma)(x+1)$ to see what's going on. – whuber May 5 '11 at 2:16
Thanks, I plotted the PDF, and it looks like the inverse of the unit exponential distribution. So if I take the negative of $\Gamma(2,1)(x+1)$ and add a constant (like 3), it should serve as an envelope for the unit exponential? – Will May 5 '11 at 16:10
Thanks for the plot! I understand it now. – Will May 5 '11 at 17:38

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.