# How to use rejection sampling to generate draws from Unit Exponential

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
}

• 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
Commented May 4, 2011 at 21:46
• Thanks for the link, but I am supposed to use the $\Gamma(2,1)$ as the envelope function.
– Will
Commented May 4, 2011 at 22:55
• One small thing, rather than constantly resizing your vector, start off with: x = numeric(100000) This will be much faster. Commented May 5, 2011 at 21:11

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

EDIT: Illustration

• +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
Commented May 4, 2011 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
Commented May 4, 2011 at 23:36
• @Will Plot $x \to pdf(\Gamma)(x+1)$ to see what's going on.
– whuber
Commented May 5, 2011 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
Commented May 5, 2011 at 16:10