[R] generating a gamma random variable

[Aboubakar Maitournam]

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