[R] generating a gamma random variable

Faheem Mitha faheem at email.unc.edu
Sun Oct 21 07:50:18 CEST 2001

Dear R People,

This question has nothing to do with R directly, but it is a simulation
question. I need to generate a random variable distributed as
gamma(\alpha,\beta), with the additional proviso that it must be a
function of random variable(s) which do not depend on \alpha, \beta. In
this case, I have \alpha = (T-1)/2, where T is a positive integer.

So, it seems reasonable to first simulate from a chi-squared distribution
with T-1 degrees of freedom. Then multiply by the appropriate scale
factor.

There seem to be at least two different ways to simulate from a
chi-squared distribution with n degrees of freedom.

1) Add together the sum of squares of n standard normal random variates.

2) Adding together n/2 independent exponential (mean 2) variates or
(n-1)/2 independent exponential variates plus a standard normal, depending
whether n is even or odd respectively.

Method 2) involves simulating roughly half as many random variates, but I
am wondering if it is really more efficient. How does simulating from an
exponential distribution compare in "expense" to simulating from a
standard normal?

I'm wondering whether this would be the best way. Does anyone else have a
better idea? One minor downside of the above method is that it requires
using a varying number (depending on T) of random variables, to obtain one
random variate. It would be preferable if this number was fixed (ideally
1).

One more question: how does one debug such an algorithm once coded up? Ie.
what is an efficient way of testing whether some random variates are
indeed distributed as Gamma(\alpha,\beta) for some \alpha, \beta?

Please cc me if you reply, I'm not subscribed to the list. Thanks in