[R] generating a gamma random variable

[Thomas Lumley]

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
> factor.

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
> 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?
>

You can obviously start by checking the mean, variance and perhaps a few
more moments.

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

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._