Suppose we want to simulate a large number n of independent exponentials with rate 1—call them X1,X2, . . . , Xn. If we were to employ the inverse transform technique we would require one logarithmic computation for each exponential generated. One way to avoid this is to first simulate Sn, a gamma random variable with parameters (n, 1) (say, by the method of Section 11.3.3).

Now interpret Sn as the time of the nth event of a Poisson process with rate 1 and use the result that given Sn the set of the first n−1 event times is distributed as the set of n−1 independent uniform (0,Sn) random variables. Based on this, explain why the following algorithm simulates n independent exponentials:

Step 1: Generate Sn, a gamma random variable with parameters (n, 1).

Step 2: Generate n −1 random numbers U1,U2, . . . , Un−1.

Step 3: Order the Ui, i = 1, . . . , n−1 to obtain U(1)

Step 4: Let U(0) = 0,U(n) = 1, and set Xi = Sn(U(i) −U(i−1)), i = 1, . . . , n.

When the ordering (step 3) is performed according to the algorithm described in Section 11.5, the preceding is an efficient method for simulating n exponentials when all n are simultaneously required. If memory space is limited, however, and the exponentials can be employed sequentially, discarding each exponential from memory once it has been used, then the preceding may not be appropriate.