5
$\begingroup$

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

1 Answer 1

5
$\begingroup$

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

EDIT: Illustration

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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