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
advance.
Sincerely, Faheem Mitha.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
generating a gamma random variable
2 messages · Faheem Mitha, 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._