generating a gamma random variable
Marco Taboga wrote:
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.
(sorry for my english which is french translation).
That method (above) is naturally the primary method for stochastic
simulation. But it's depending of many conditions and it is
computationally delicate.
There is a way to simulate a gamma law of parameters a and b, based on
the standard method above (simplified) and the variation of the method
of rejection.
1) First you simulate a random variable following gamma law with
parameters a and 1 with 0<a<1.
For this you use the fact that the density of law gamma(a,1),
g(a,1)(x) verifies the inequality :
g(a,1)(x) <= ((a+e)/(ae*gamma(a))g(x) where g(x) is the function
g(x)=(ae/(a+e))(x^(a-1)) 1{ 0<=x<1}+e^(-x) 1{x>=1}
1_1 You simulate simply a random variable Y with density g by
leading Y=G^(-1)(U), U uniform on [0,1] and
G^{-1}(z)=[((a+e)/e)z]^(1/a) 1{z<(e/(a+e))}-log((1-z)((a+e)/ae)) 1
{z>=e/(a+e)}
1-2 Simulation of random uniform variable V on [0,1] independant of
U
1-3 if V <= gamma(a,1)(Y)/[((a+e)/(a*e*gamma(a)))*g(Y)] more
precisely
if V <= e^(-Y)1{0<=y<1} + Y^(a-1)1{1<=y} we lead X=Y or if V >....... we
return to 1_1.
2) You simulate a random variable X following gamma(a,b) a >0, b>0 as
folow
2-1 consider a-[a] the fractionnary part of (a-[a] <1), we simulate a
random variable H following law gamma(a-[a],1) according to the
algorithm 1
above.
2-1 We simulate a random variable J of law gamma([a],1) as sum of [a]
exponential variables
J=-log[U1 x ...........U[a]] (product in log of U from 1 to [a])
where the U[i] are uniform, independant
on [0,1]
2-3 lead X=b(J+H)
Aboubakar Maitournam.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://stat.ethz.ch/pipermail/r-help/attachments/20011024/8bfa61ad/attachment.html