generating a gamma random variable

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

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.

And there are source code and references for how R does it, yet another
different way, in src/nmath/rgamma.c

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