Robert_V._Hogg,_Joseph_W._McKean,_Allen_T._Craig

4.8. The Method of Monte Carlo 295

Then the corresponding stream of exponential observations is

0 .641, 1.95, 0.696, 1.13, 0. 274.

As the next example shows, we can generate Poisson realizations using this expo-

nential generation.

Example 4.8.2(Simulating Poisson Processes).LetXbe the number of occur-

rences of an event over a unit of time and assume that it has a Poisson distribution

with meanλ, (3.2.1). LetT 1 ,T 2 ,T 3 ,...be the interarrival times of the occurrences.

Recall from Remark 3.3.1 thatT 1 ,T 2 ,T 3 ,...are iid with the common Γ(1, 1 /λ)-

distribution. Note thatX=kif and only if

∑k

j=1Tj≤1and

∑k+1

j=1Tj>1. Using

this fact and the generation of Γ(1, 1 /λ) variates discussed above, the following

algorithm generates a realization ofX (assume that the uniforms generated are

independent of one another).

1. SetX=0andT=0.


2. GenerateUuniform (0,1) and letY=−(1/λ) log(1−U).


3. SetT=T+Y.


4. IfT>1, outputX;
   else setX=X+1andgotostep2.
   The R functionpoisrandprovides an implementation of this algorithm, generating
   nsimulations of a Poisson distribution with parameterλ. As an illustration, we
   obtained 1000 realizations from a Poisson distribution withλ= 5 by running R
   with the R codetemp = poisrand(1000,5), which stores the realizations in the
   vectortemp. The sample average of these realizations is computed by the command
   mean(temp). In the situation that we ran, the realized mean was 4.895.
   Example 4.8.3(Monte Carlo Integration).Suppose we want to obtain the integral
   ∫b
   ag(x)dxfor a continuous functiongover the closed and bounded interval [a, b].
   If the antiderivative ofgdoes not exist, then numerical integration is in order. A
   simple numerical technique is the method of Monte Carlo. We can write the integral
   as ∫
   b

a

g(x)dx=(b−a)

∫b

a

g(x)

1

b−a

dx=(b−a)E[g(X)],

whereXhas the uniform (a, b) distribution. The Monte Carlo technique is then to
generate a random sampleX 1 ,...,Xnof sizenfrom the uniform (a, b) distribution

and computeYi=(b−a)g(Xi). ThenYis an unbiased estimator of

∫b
ag(x)dx.
Example 4.8.4(Estimation ofπby Monte Carlo Integration).For a numerical
example, reconsider the estimation ofπ. Instead of the experiment described in
Example 4.8.1, we use the method of Monte Carlo integration. Letg(x)=4

√

1 −x^2

for 0<x<1. Then

π=

∫ 1

0

g(x)dx=E[g(X)],

whereXhas the uniform (0,1) distribution. Hence we need to generate a random
sampleX 1 ,...,Xnfrom the uniform (0,1) distribution and formYi=4

√

1 −Xi^2.