[R] generating a gamma random variable
tlumley at u.washington.edu
Sun Oct 21 22:07:26 CEST 2001
On Sun, 21 Oct 2001, Faheem Mitha wrote:
> 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
Or simulate from a gamma (using rgamma()) and multiply
> 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?
Both 1 and 2 are more expensive than rchisq() or rgamma(), at least if T
is large. A simple experiment shows that rnorm takes about 2-3 times as
long as rexp, but that 10^6 of them still only takes a few seconds.
In general the answer to the question "Which is quicker?" is "Try it and
see" (and the most common result is "It doesn't matter").
> 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
> 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?
You can obviously start by checking the mean, variance and perhaps a few
One more general possibility is to use the large deviation bound used in
checking the built-in generators, in tests/p-r-random-tests.R. The
empirical CDF F_n from a sample of size n and the true CDF F are related by
P(sup | F_n - F | > t) < 2 exp(-2nt^2)
If you get a large enough sample n this should detect any error in the
distribution (provided you have F computed correctly). However, AFAIK this
test has never actually detected an error in any of our random number
generators so I don't know how useful it is in practice.
If there are special cases of your simulation where you know the correct
answer this can also help with checking.
Thomas Lumley Asst. Professor, Biostatistics
tlumley at u.washington.edu University of Washington, Seattle
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
More information about the R-help