<!doctype html public "-//w3c//dtd html 4.0 transitional//en">
<html>
<body bgcolor="#FFFFFF">
Marco Taboga wrote:
<blockquote TYPE=CITE><style></style>
<font face="Arial"><font size=-1>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)</font></font>
<font face="Arial"><font size=-1>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).</font></font>
<font face="Arial"><font size=-1>Sincerely,Marco.</font></font></blockquote>
<p><br>(sorry for my english which is french translation).
<br>That method (above) is naturally the primary method for stochastic
simulation. But it's depending of many conditions and it is computationally
delicate.
<br>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.
<br>1) First you simulate a random variable following gamma law with parameters
a and 1 with 0<a<1.
<p>For this you use the fact that the density of law gamma(a,1),
g(a,1)(x) verifies the inequality :
<br>g(a,1)(x) <= ((a+e)/(ae*gamma(a))g(x) where g(x) is the function
<br>g(x)=(ae/(a+e))(x^(a-1)) 1{ 0<=x<1}+e^(-x) 1{x>=1}
<br> 1_1 You simulate simply a random variable Y with density
g by leading Y=G^(-1)(U), U uniform on [0,1] and
<br> 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)}
<br> 1-2 Simulation of random uniform variable V on [0,1]
independant of U
<br> 1-3 if V <= gamma(a,1)(Y)/[((a+e)/(a*e*gamma(a)))*g(Y)]
more precisely
<br>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.
<p>2) You simulate a random variable X following gamma(a,b) a >0, b>0 as
folow
<br>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
<br>above.
<br>2-1 We simulate a random variable J of law gamma([a],1) as sum of [a]
exponential variables
<br> J=-log[U1 x ...........U[a]] (product in log of
U from 1 to [a]) where the U[i] are uniform, independant
<br> on [0,1]
<br>2-3 lead X=b(J+H)
<br>
<p>Aboubakar Maitournam.
</body>
</html>