Tài liệu Random Numbers part 3 pptx

7.2 Transformation Method: Exponential and Normal Deviates 287

visit website http://www.nr.com or call 1-800-872-7423 (North America only),

or send email to [email protected] (outside North America).

readable files (including this one) to any server

computer, is strictly prohibited. To order Numerical Recipes books,

diskettes, or CDROMs

Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine￾Copyright (C) 1988-1992 by Cambridge University Press.

Programs Copyright (C) 1988-1992 by Numerical Recipes Software.

Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)

7.2 Transformation Method: Exponential and

Normal Deviates

In the previous section, we learned how to generate random deviates with

a uniform probability distribution, so that the probability of generating a number

between x and x + dx, denoted p(x)dx, is given by

p(x)dx =

n dx 0 <x< 1

0 otherwise (7.2.1)

The probability distribution p(x) is of course normalized, so that

Z ∞

−∞

p(x)dx =1 (7.2.2)

Now suppose that we generate a uniform deviate x and then take some prescribed

function of it, y(x). The probability distribution of y, denoted p(y)dy, is determined

by the fundamental transformation law of probabilities, which is simply

|p(y)dy| = |p(x)dx| (7.2.3)

or

p(y) = p(x)

dx

dy

(7.2.4)

Exponential Deviates

As an example, suppose that y(x) ≡ − ln(x), and that p(x) is as given by

equation (7.2.1) for a uniform deviate. Then

p(y)dy =

dx

dy

dy = e−ydy (7.2.5)

which is distributed exponentially. This exponential distribution occurs frequently

in real problems, usually as the distribution of waiting times between independent

Poisson-random events, for example the radioactive decay of nuclei. You can also

easily see (from 7.2.4) that the quantity y/λ has the probability distribution λe−λy.

So we have

#include <math.h>

float expdev(long *idum)

Returns an exponentially distributed, positive, random deviate of unit mean, using

ran1(idum) as the source of uniform deviates.

{

float ran1(long *idum);

float dum;

do

dum=ran1(idum);

while (dum == 0.0);

return -log(dum);

}