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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Sun, 21 Oct 2001, Bob O'Hara wrote:> 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. > > > Is there a reason why you don't want to use either rgamma or rchisq?Come to think of it, there is no reason why I should not use rchisq. This would depend on T, but that is Ok, since T is a constant within the calculation, and the method I'm using above depends on T being fixed anyway. This would probably be the most efficient way, since to get a gamma distribution would then just involve multiplying by a constant. Background: I'm using this simulation to run coupled Markov chains. T is a contant within the calculation.> > There seem to be at least two different ways to simulate from a > > chi-squared distribution with n degrees of freedom. > > > There are much quicker ways than you suggested. A rejection algorithm > is one choice. I'll save him the embarrassment of having to indulge in > self promotion by suggesting that the classic text is Brian Ripley's > 1987 book 'Stochastic Simulation'. There is also an algorithm in Press > et al.'s 'Numerical Recipes in C'.Yes, actually Stochastic Simulation is on my desk as I write, and it suggests the methods I outlined above for simulating from chisq/gamma. I'm still wondering about how to check (in general) whether a bunch of rvs are distributed as gamma or not. How do people check such things in practice? Thanks for your reply. 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
One method to generate a value x taken by a random variable X with a strictly
increasing distribution function F(x) is the following:
- Generate a value y taken by a random variable Y distributed as uniform on
[0,1] (every programming language has routines to do this)
- Evaluate the inverse of the distribution function F at the point y you
just generated
- The number you get x=F^(-1)(y) is the value x you wanted to generate
(drawn from the distribution you chose)
This method is nice, but it can be computationally very expensive, since, in
most cases (as in the case of a gamma), to compute the inverse of F you have to
solve numerically the equation F(x)-y=0 and every iteration of the method you
use to solve the equation requires the numerical computation of a definite
integral to find F(x).
The reward you get for this great amount of computation is that your random
generator for X is as good as the random generator for Y (the uniform random
variable).
Sincerely,
Marco.
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
https://stat.ethz.ch/pipermail/r-help/attachments/20011022/ed049ac9/attachment.html