2. SYSTEMS HAVING POISSON ARRIVALS

Consider a system in which arrivals occur according to a Poisson process and suppose that we are interested in using simulation to compute E [D], where the value of D depends on the arrival process only through those arrivals before time t. For instance, D might be the sum of the delays of all arrivals by time t in a parallel multiserver queuing system. We suggest the following approach to using simulation to estimate E [D]. First, with N(t) equal to the number of arrivals by time t, note that

Each run of our suggested simulation procedure generates an independent estimate of E [D]. At each stage of a run, we will have a set S whose elements are arranged in increasing value and which represents the set of arrival times. For simplicity, we will present our approach under the assumption that E [D|N(t) = 0] can be easily computed and also that D can be determined by knowing the arrival times along with the service time of each arrival. A run is as follows:

1. Let N = 1. Generate a random number U₁ and let S = {tU₁}.
2. Suppose N(t) = 1, with the arrival occurring at time tU₁. Generate the service time of this arrival and compute the resulting value of D. Call this value D₁.
3. Let N = N + 1.
4. Generate a random number U_N and add tU_N in its appropriate place to the set S so that the elements is S are in increasing order.
5. Suppose N(t) = N, with S specifying the N arrival times; generate the service time of the arrival at time tU_N, and using the previously generated service times of the other arrivals, compute the resulting value of D. Call this value D_N.
6. If N < m, return to Step 3. If N = m, use the inverse transform method to generate the value of N(t) conditional on it exceeding m. If the generated value is m + k, generate k additional random numbers, multiply each by t, and add these k numbers to the set S. Generate the service times of these k arrivals and, using the previously generated service times, compute D. Call this value D_>m.

With D₀ = E [D|N(t) = 0], the estimate from this run is

Using the fact that the set of unordered arrival times, given that N(t) = j, is distributed as a set of j independent uniform (0,t) random variables, it follows that the preceding is an unbiased estimator of E [D]. Generating multiple runs and taking the average value of the resulting estimates yields the final simulation estimator.

We now show that EST has a smaller variance than does the raw simulation estimator D.

Theorem 1:

Proof: The quantity D can be simulated as follows:

It is now easy to check, by conditioning on N(t), that

The result now follows from the conditional variance formula. █

Remarks:

1. The use of Eq. (1) is the use of stratified sampling (see [3]), which has been previously suggested when analyzing Poisson arrival systems (see, for instance, [1, p. 228]). However, previous suggestions were to first estimate E [D|N(t) = 1] by a sequence of independent runs; then to estimate E [D|N(t) = 2] by a new sequence of independent runs (that are also independent of the runs used to compute E [D|N(t) = 1]), and so on. This differs from our idea of estimating all of the quantities E [D|N(t) = j],j = 1,…,m, and E [D|N(t) > m] in a single run, sequentially making use of the set of generated data values used to estimate E [D|N(t) = j] to speed the estimation of E [D|N(t) = j + 1], and so on.
2. It should be noted that the variance of our estimator
   
   is, because of the positive correlations introduced by reusing the same data, larger than it would be if the D_j were independent estimators. However, we believe that the increased speed of the simulation will usually more than make up for this increased variance.
3. When computing D_j+1, we can make use of quantities used in computing D_j. For instance, suppose D_i,j was the delay of arrival i when N(t) = j. Then if the new arrival time tU_j+1 is the kth smallest of the new set S, then D_i,j+1 = D_i,j for i < k.
4. Other variance reduction ideas can be used in conjunction with our approach. For instance, we can improve the estimator by using a linear combination of the service times as a control variable.

It remains to determine an appropriate value of m. A reasonable approach might be to choose m to make

sufficiently small. Because Var(N(t)) = λt, a reasonable choice would be of the form

for some positive integer k. One can often bound E [D|N(t) > m] and use this bound to determine the appropriate value of m. For instance, suppose D is the sum of the delays of all arrivals by time t in a single server system with mean service time 1. Then because this quantity will be maximized when all arrivals come simultaneously, we see that

Because the conditional distribution of N(t) given that it exceeds m will, when m is at least five standard deviations greater than E [N(t)], put most of its weight near m + 1, we see from the preceding equation that one can reasonably assume that

Using that for a standard normal random variable Z (see [2]),

we see, upon using the normal approximation to the Poisson, that for

, we can reasonably assume that

For instance, with λt = 10³ and k = 6, the preceding upper bound is about 0.0008.

3. COMPUTING THE EXPECTED VALUE OF AN INCREASING FUNCTION

Suppose that we want to use simulation to compute E [g(U₁,…,U_n)], where U₁,…,U_n are independent uniform(0,1) random variables and g is nondecreasing in each of its coordinates. Because of the monotonicity of g,

will often be a good predictor of g in the sense that

will be relatively small. Because of this, we suggest generating U₁,…,U_n by first generating the value of

and then generating U₁,…,U_n conditional on the generated value of the product. In generating the value of

on different simulation runs, we suggest using stratified sampling.

To implement the preceding idea, we need to first show how to generate

and how to generate U₁,…,U_n conditional on

. Note the following:

Now the generation can proceed as follows:

Let G_n denote the distribution function of a gamma random variable with parameters (n,1). Suppose that we are planning to do m simulation runs. Then on the kth simulation run, a random number U should be generated and T should be taken to equal G_n⁻¹((U + k − 1)/m) (i.e., we use stratified sampling when generating the successive values of T). We should then follow Steps 2 and 3 of the preceding and then for the values of U_i obtained in Step 3, compute g(U₁,…,U_n). The average of the values of g(U₁,…,U_n) obtained in the m runs is then taken as the estimator of E [g(U₁,…,U_n)].

Remarks:

3. Another way to do the stratification is to choose r values 0 = a₀ < a₁ < ··· < a_r < a_r+1 = ∞. Determine p_i = P{a_i−1 < T < a_i}, i = 1,…,r + 1, and then perform mp_i of the m planned simulation runs with T simulated conditional on its being in the interval (a_i−1,a_i). This conditional value can be simulated by using the rejection procedure based on an exponential random variable that is conditioned to lie in the same interval. Indeed, even better is to do a small preliminary simulation to estimate the quantities Var(g(U₁,…,U_n)|a_i−1 < T < a_i). If s_i²,i = 1,…,r + 1, are the estimates, then do a total of m(p_i s_i²/[sum ]_j p_j s_j²) of the runs conditional on a_i−1 < T < a_i, either by using the inverse transform or the rejection method. (If the inverse transform is used, then you should stratify within the subintervals.) The final estimate is

, where X_j is the average of the simulation runs done conditional on a_i−1 < T < a_i.

4. If g(u₁,…,u_n) is only increasing in k of its variables, say in u₁,…,u_k, then we can generate U₁,…,U_n on successive runs by using stratified sampling on

, generating U₁,…,U_k conditional on the value of this product and generating U_i,i > k, as independent uniforms.