8 Random-Variate Generation

This chapter deals with procedures for sampling from a variety ofwidely used continuous and discrete distributions. Previous discussions and examples indicated the usefulness of probability distributions in modeling activities that are generally unpredictable or uncertain. For example, interarrival times and service times at queues and demands for a product are quite often unpredictable in nature, a least to a certain extent. Usually, such variables are modeled as random variables with some specified probability distribution, and standard statistical procedures exist for estimating the parameters of the hypothesized distribution and for testing the validity of the assumed statistical model.

In this chapter, it is assumed that a distribution has been completely specified, and ways are sough to generate samples from this distribution to be used as input to a simulation model. The purpose of the chapter is to explain and illustrate some widely used techniques for generating random variates, not to give a state-of-the-art survey of the most efficient techniques. In practice, most simulation modelers will use either existing routines available in programming libraries or the routines built into a simulation language. However, some programming languages do not have built-in routines for all of the regularly used distributions, and some computer installations do not have random-variategeneration libraries; in such cases the modeler must construct an acceptable routine. Even though the chance of this happening is small, it is nevertheless worthwhile to understand how random-variate generation occurs.

This chapter discusses the inverse-transform technique and, more briefly, the acceptance-rejection technique and special properties. Another technique, the composition method, is discussed by Devroye [1986], Dagpunar [1988], and Law [2007]. All the techniques in this chapter assume that a source of uniform [0,1] random numbers,$R_{1} , R_{2} , \ldots$. is readily available, where each $R_{i}$ has pdf

and cdf

Throughout this chapter R and $R_{1} , R_{2} , \ldots$. represent random numbers uniformly distributed on [0,1] and generated by random-number generation techniques or taken from a random number table, such as Table A.1 in the Appendix. Many of the algorithms in this chapter have been translated into Visual Basic for Applications and can be found in the spreadsheet RandomNumberTools.xls at www.bcnn.net.

1 Inverse-Transform Technique

The inverse-transform technique can be used to sample from the exponential, uniform, Weibull, and triangular distributions, and from empirical distributions. Additionally, it is the underlying principle for sampling from a wide variety of discrete distributions. The technique will be explained in detai for the exponential distribution and then applied to other distributions. Computationally, it is the most straightforward, but not always the most efficient, technique.

1.1 Exponential Distribution

The exponential distribution has the probability density function (pdf)

and the cumulative distribution function (cdf)

The parameter λ can be interpreted as the mean number of occurrences per time unit. For example, if interarrival times $X_{1} , X_{2} , X_{3} , . . .$. had an exponential distribution with rate λ, then λ could be interpreted as the mean number of arrivals per time unit, or the arrival rate. Notice that, for any i,

and so $1 / \lambda$ is the mean interarrival time. The goal here is to develop a procedure for generating values $X_{1} , X_{2} , X_{3} , . . .$. that have an exponential distribution.

The inverse-transform technique can be utilized, at least in principle, for any distribution, but it is most useful when the cdf,$F ( x )$), is of a form so simple that its inverse,$F^{- 1}$, can be computed easily. (The notation $F^{- 1}$ denotes the solution of the equation $r = F ( {\boldsymbol {x}} )$ for x; it does not denote $1 / F . )$ One step-by-step procedure for the inverse-transform technique, illustrated by the exponentia distribution, is shown here:

Step 1. Compute the cdf of the desired random variable $X$ For the exponential distribution, the cdf is $F ( x ) = 1 - e^{- \lambda x} , x \geq 0$ Step 2. Set $F ( X ) = R$ on the range of X.

For the exponential distribution, it becomes $1 - e^{- \lambda X} = R$ on the range $X \geq 0$ X is a random variable (with the exponential distribution in this case), so $1 - e^{- \lambda X}$ is also a random variable, here called R. As will be shown later, R has a uniform distribution over the interval [0, 1].

Step 3. Solve the equation $F ( X ) = R$ for X in terms of R.

For the exponential distribution, the solution proceeds as follows:

Equation (1) is called a random-variate generator for the exponential distribution. In general, Equation (1) is written as $X = F^{- 1} ( R )$. Generating a sequence of values is accomplished through Step 4.

Step 4. Generate (as needed) uniform random numbers $R_{1} , R_{2} , R_{3} , . . .$. and compute the desired random variates by

For the exponential case,$F^{- 1} ( R ) = ( - 1 / \lambda ) \ln ( 1 - R )$ by Equation (1), so

for $i = 1 , 2 , 3 , \dots$. One simplification that is usually employed in Equation (2) is to replace $1 - R_{i}$ by $R_{i}$ to yield

This alternative is justified by the fact that both $R_{i}$ and $1 - R_{i}$ are uniformly distributed on [0,1].

Example 1

Table 1 gives a sequence of random numbers from Table A.1 and the computed exponential variates,$X_{i} ,$, given by Equation (2) with the value $\lambda = 1$. Figure 1(a) is a histogram of 200 values,$R_{1} , R_{2} , \ldots , R_{200}$ from the uniform distribution, and Figure 1(b) is a histogram of the 200 values,$X_{1} , X_{2} , \ldots , X_{200}$, computed by Equation (2). Compare these empirical histograms with the theoretical density functions in Figure 1(c) and (d). As illustrated here, a histogram is an estimate of the underlying density function.

Table 1 Generation of Exponential Variates $X_{i}$ with Mean 1, given Random Numbers $R_{i}$| i | 1 | 2 | 3 | 4 | 5 |

Figure 2 gives a graphical interpretation of the inverse-transform technique. The cdf shown is $F ( x ) = 1 - e^{- x}$, an exponential distribution with rate $\lambda = 1$. To generate a value $X_{1}$ with cdf $F ( x )$ a random number $R_{1}$ between 0 and 1 is generated, then a horizontal line is drawn from $R_{1}$ to the graph of the cdf, then a vertical line is dropped to the x axis to obtain $X_{1}$, the desired result. Notice the inverse relation between $R_{1}$ and $X_{1}$, namely

and

In general, the relation is written as

and

Why does the random variable $X_{1}$ generated by this procedure have the desired distribution? Pick a value $x_{0}$ and compute the cumulative probability

To see the first equality in Equation (4), refer to Figure 2, where the fixed numbers $x_{0}$ and $F ( x_{0} )$ are drawn on their respective axes. It can be seen that $X_{1} \leq x_{0}$ when and only when $R_{1} \leq F ( x_{0} )$). Since $0 \leq F ( x_{0} ) \leq 1$, the second equality in Equation (4) follows immediately from the fact that $R_{1}$ is uniformly distributed on [0,1]. Equation (4) shows that the cdf of $X_{1}$ is $F ;$ hence,$X_{1}$ has the desired distribution.

Figure 1 (a) Empirical histogram of 200 uniform random numbers; (b) empirical histogram of 200 exponential variates; (c) theoretical uniform density on [0,1]; (d) theoretical exponentia density with mean 1.

1.2 Uniform Distribution

Consider a random variable X that is uniformly distributed on the interval [a, b]. A reasonable guess for generating X is given by

Figure 2 Graphical view of the inverse-transform technique.

(Recall that R is always a random number on [0, 1].) The pdf of X is given by

The derivation of Equation (5) follows Steps 1 through 3 of Section 1.1:

Step 1. The cdf is given by

Step 2. Set $F ( X ) = ( X - a ) / ( b - a ) = R$ Step 3. Solving for X in terms of R yields $X = a + ( b - a ) R ,$, which agrees with Equation (5).

1.3 Weibull Distribution

The Weibull distribution was introduced as a model for time to failure for machines or electronic components. When the location parameter ν is set to 0:

where $\alpha > 0$ and $\beta > 0$ are the scale and shape parameters of the distribution. To generate a Weibul variate, follow Steps 1 through 3 of Section 1.1:

Step 1. The cdf is given by $F ( X ) = 1 - e ^ {- ( x / \alpha ) ^ {\beta}} , x \geq 0$ Step 2. Set $F ( X ) = 1 - e ^ {- ( X / \alpha ) ^ {\beta}} = R .$ Step 3. Solving for X in terms of R yields

The derivation of Equation (6) is left as Exercise 10 for the reader. By comparing Equations (6) and (1), it can be seen that, if X is a Weibull variate, then $X^{\beta}$ is an exponential variate with mean $\alpha^{\beta}$ Conversely, if Y is an exponential variate with mean $\mu ,$, then $Y^{1 / \beta}$ is a Weibull variate with shape parameter $\beta$ and scale parameter $\alpha = \mu^{1 / \beta}$

1.4 Triangular Distribution

Consider a random variable X that has pdf

as shown in Figure 3. This distribution is called a triangular distribution with endpoints $( 0 , 2 )$ and mode at 1. Its cdf is given by

For $0 \leq X \leq 1$

Figure 3 Density function for a particular triangular distribution.

and for $1 \leq X \leq 2$

In Equation $( 7 ) , 0 \leq X \leq 1$ implies that $\begin{array} {r} {0 \leq R \leq \frac {1} {2}} \end{array}$, in which case $X = {\sqrt {2 R}}$. By Equation (8),$1 \leq X \leq 2$ implies that $\begin{array} {r} {\frac {1} {2} \leq R \leq 1} \end{array}$, in which case $X = 2 - \sqrt {2 ( 1 - R )}$. Thus, X is generated by

Exercises 2, 3, and 4 give the student practice in dealing with other triangular distributions. Notice that, if the pdf and cdf of the random variable X come in parts (i.e., require different formulas over different parts of the range of X), then the application of the inverse-transform technique for generating X will result in separate formulas over different parts of the range of R, as in Equation (9).

1.5 Empirical Continuous Distributions

If the modeler has been unable to find a theoretical distribution that provides a good model for the input data, then it may be necessary to use the empirical distribution of the data. One possibility is to simply resample the observed data itself. This is known as using the empirical distribution, and it makes particularly good sense when the input process is known to take on a finite number of values. See Section 1.7 for an example of this type of situation and for a method for generating random inputs.

On the other hand, if the data are drawn from what is believed to be a continuous-valued inpu process, then it makes sense to interpolate between the observed data points to fill in the gaps. This section describes a method for defining and generating data from a continuous empirical distribution.

Five observations of fire-crew response times (in minutes) to incoming alarms have been collected to be used in a simulation investigating possible alternative staffing and crew-scheduling policies. The data are

Before collecting more data, it is desired to develop a preliminary simulation model that uses a response-time distribution based on these five observations. Thus, a method for generating random variates from the response-time distribution is needed. Initially, it will be assumed that response times X have a range $0 \leq X \leq c .$, where c is unknown, but will be estimated by $\widehat {c} = \operatorname* {m a x} \{X_{i} : i =$$1 , \ldots , n \} = 2 . 76 \quad$, where $\{X_{i} , i = 1 , \dots , n \}$ are the raw data and $n = 5$ is the number of observations.

Arrange the data from smallest to largest and let $x_{( 1 )} \leq x_{( 2 )} \leq \cdot \cdot \cdot \leq x_{( n )}$ denote these sorted values. The smallest possible value is believed to be 0, so define $x_{( 0 )} = 0$. Assign the probability $1 / n = 1 / 5$ to each interval $x_{( i - 1 )} < x \leq x_{( i )}$, as shown in Table 2. The resulting empirical cdf,$\widehat F ( x )$ is illustrated in Figure 4. The slope of the ith line segment is given by

The inverse cdf is calculated by

when $( i - 1 ) / n < R \leq i / n .$ For example, if a random number $R_{1} = 0 . 71$ is generated, then $R_{1}$ is seen to lie in the fourth interval (between $3 / 5 = 0 . 60$ and $4 / 5 = 0 . 80 )$; so, by Equation (10),

The reader is referred to Figure 4 for a graphical view of the generation procedure.

Table 2 Summary of Fire-Crew Response-Time Data

Figure 4 Empirical cdf of fire-crew response times.

In Example 2, each data point was represented in the empirical cdf. If a large sample of data is available (and sample sizes from several hundred to tens of thousands are possible with automated data collection), then it might be more convenient (and computationally efficient) to first summarize the data into a frequency distribution with a much smaller number ofintervals and then fit a continuous empirical cdf to the frequency distribution. Only a slight generalization of Equation (10) is required to accomplish this. Now the slope of the ith line segment is given by:

where $c_{i}$ is the cumulative probability of the first i intervals of the frequency distribution and $x_{( i - 1 )} <$$x \leq x_{( i )}$ is the ith interval. The inverse cdf is calculated as

when $c_{i - 1} < R \leq c_{i}$ Suppose that 100 machine repair times have been collected. The data are summarized in Table 3 in terms of the number of observations in various intervals. For example, there were 31 observations between 0 and 0.5 hour, 10 between 0.5 and 1 hour, and so on. Suppose it is known that all repairs take at least 15 minutes, so that $X \geq 0 . 25$ hours always. Then we set $x_{( 0 )} = 0 . 25$, as shown in Table 3 and Figure 5.

Table 3 Summary of Repair-Time Data

Figure 5 Generating variates from the empirical distribution function for repair-time data $( X \geq 0 . 25 )$ For example, suppose the first random number generated is $R_{1} = 0 . 83$. Then, since $R_{1}$ is between $c_{3} = 0 . 66$ and $c_{4} = 1 . 00$

As another illustration, suppose that $R_{2} = 0 . 33$. Since $c_{1} = 0 . 31 < R_{2} \leq 0 . 41 = c_{2}$

The point $( R_{2} = 0 . 33 , X_{2} = 0 . 6 )$ is also shown in Figure 5.

Now reconsider the data of Table 3. The data are restricted in the range $0 . 25 \leq X \leq 2 . 0$ but the underlying distribution might have a wider range. This provides one important reason for attempting to find a theoretical probability distribution (such as the gamma or Weibull) for the data: that these distributions do allow a wider range—namely,$0 \leq X < \infty$. On the other hand, an empirica distribution adheres closely to what is present in the data itself, and the data are often the best source of information available.

When data are summarized in terms of frequency intervals, it is recommended that relatively short intervals be used, for doing so results in a more accurate portrayal of the underlying cdf. For example, for the repair-time data of Table 3, for which there were $n = 100$ observations, a much more accurate estimate could have been obtained by using 10 to 20 intervals, certainly not an excessive number, rather than the four fairly wide intervals actually used here for purposes of illustration.

Several comments are in order:

1. A computerized version of this procedure will become inefficient as the number of intervals n increases. A systematic computerized version is often called a table-lookup generation scheme, because, given a value of R, the computer program must search an array of $c_{i}$ values to find the interval i in which R lies, namely the interval i satisfying

The more intervals there are, the longer, on average, the search will take if it is implemented in the crude way described here. The analyst should consider this tradeoff between accuracy of the estimating cdf and computational efficiency when programming the procedure. If a large number of observations are available, the analyst may well decide to group the observations into intervals from 20 to 50, say, and then use the procedure of Example 3—or a more efficient table-lookup procedure could be used, such as the one described in Law [2007].

2. In Example 2, it was assumed that response times X satisfied $0 \leq X \leq 2 . 76$. This assumption led to the inclusion of the points $x_{( 0 )} = 0$ and $x_{( 5 )} = 2 . 76$ in Figure 4 and Table 2. If it is known a priori that X falls in some other range—for example, if it is known that response times are always between 15 seconds and 3 minutes, that is,

—then the points $x_{( 0 )} = 0 . 25$ and $x_{( 6 )} = 3 . 0$ would be used to estimate the empirical cdf of response times. Notice that, because of inclusion of the new point $x_{( 6 )}$, there are now six intervals instead of five, and each interval is assigned probability $1 / 6 = 0 . 167$. Exercise 12 illustrates the use of these additional assumptions.

1.6 Continuous Distributions without a Closed-Form Inverse

A number of useful continuous distributions do not have a closed form expression for their cdf or its inverse; examples include the normal, gamma, and beta distributions. For this reason, it is often stated that the inverse-transform technique for random-variate generation is not available for these distributions. It can, in effect, become available if we are willing to approximate the inverse cdf, or numerically integrate and search the cdf. Although this approach sounds inaccurate, notice that even a closed-form inverse requires approximation to evaluate it on a computer. For example, generating exponentially distributed random variates via the inverse cdf $X = F^{- 1} ( R ) = - \ln ( 1 - R ) / \lambda$ requires a numerical approximation for the logarithm function. Thus, there is no essential difference between using an approximate inverse cdf and approximately evaluating a closed-form inverse. The problem with using approximate inverse cdfs is that some of them are computationally slow to evaluate.

To illustrate the idea, consider a simple approximation to the inverse cdf of the standard normal distribution, proposed by Schmeiser [1979]:

This approximation gives at least one-decimal-place accuracy for $0 . 0013499 \leq R \leq 0 . 9986501$ Table 4 compares the approximation with exact values (to four decimal places) obtained by numerical integration for several values of R. Much more accurate approximations exist that are only slightly more complicated. A good source of these approximations for a number of distributions is Bratley, Fox, and Schrage [1996].

Table 4 Comparison of Approximate Inverse with Exact Values (to four decimal places) for the Standard Normal Distribution

1.7 Discrete Distributions

All discrete distributions can be generated via the inverse-transform technique, either numerically through a table-lookup procedure or, in some cases, algebraically, the final generation scheme being in terms of a formula. Other techniques are sometimes used for certain distributions, such as the convolution technique for the binomial distribution. Some of these methods are discussed in later sections. This subsection gives examples covering both empirical distributions and two of the standard discrete distributions, the (discrete) uniform and the geometric. Highly efficient table-lookup procedures for these and other distributions are found in Bratley, Fox, and Schrage [1996] and in Ripley [1987].

Example 4: An Empirical Discrete Distribution

At the end of any day, the number of shipments on the loading dock of the IHW Company is either 0, 1, or 2, with observed relative frequency of occurrence of 0.50, 0.30, and 0.20, respectively. Interna consultants have been asked to develop a model to improve the efficiency of the loading and hauling operations; as part of this model, they will need to be able to generate values X to represent the number of shipments on the loading dock at the end of each day. The consultants decide to model X as a discrete random variable with the distribution given in Table 5 and shown in Figure 6.

Table 5 Distribution of Number of Shipments X

The probability mass function (pmf), p(x), is given by

and the cdf,$F ( x ) = P ( X \leq x )$, is given by

Recall that the cdf of a discrete random variable always consists of horizontal line segments with jumps of size p(x) at those points, x, that the random variable can assume. For example, in Figure 6, there is a jump of size $p ( 0 ) = 0 . 5 \mathrm{a t} x = 0 .$, of size $p ( 1 ) = 0 . 3 \mathrm{a t} x = 1$, and of size $p ( 2 ) = 0 . 2$ at $x = 2$ Figure 6 Cdf of number of shipments X.

For generating discrete random variables, the inverse-transform technique becomes a tablelookup procedure, but, unlike in the case of continuous variables, interpolation is not required. To illustrate the procedure, suppose that $R_{1} = 0 . 73$ is generated. Graphically, as illustrated in Figure $^ {6 ,}$ first locate $R_{1} = 0 . 73$ on the vertical axis, next draw a horizontal line segment until it hits a $\mathrm {\bf {\ddot {j} u m p}} ^ {, \mathrm{*}}$ in the cdf, and then drop a perpendicular to the horizontal axis to get the generated variate. Here $R_{1} = 0 . 73$ is transformed to $X_{1} = 1$. This procedure is analogous to the procedure used for empirica continuous distributions in Section 1.5 and illustrated in Figure 4, except that the final step, linear interpolation, is eliminated.

The table-lookup procedure is facilitated by construction of a table such as Table 6. When $R_{1} = 0 . 73$ is generated, first find the interval in which $R_{1}$ lies. In general, for $R = R_{1}$, if

then set $X_{1} = x_{i}$. Here,$r_{0} = 0 ; x_{0} = - \infty ; x_{1} , x_{2} , \ldots , x_{n}$ are the possible values of the random variable; and $r_{k} = p ( x_{1} ) + \cdot \cdot \cdot + p ( x_{k} ) , k = 1 , 2 , \cdot \cdot \cdot , n .$ For this example,$n = 3 , x_{1} = 0 , x_{2} = 1$, and $x_{3} = 2 ;$ hence,$r_{1} = 0 . 5 , r_{2} = 0 . 8$, and $r_{3} = 1 . 0$. (Notice that $r_{n} = 1 . 0$ in all cases.)

Table 6 Table for Generating the Discrete Variate X

Since $r_{1} = 0 . 5 < R_{1} = 0 . 73 \leq r_{2} = 0 . 8$, set $X_{1} = x_{2} = 1$. The generation scheme is summarized as follows:

Example 4 illustrated the table-lookup procedure; the next example illustrates an algebraic approach that can be used for certain distributions.

Example 5: A Discrete Uniform Distribution

Consider the discrete uniform distribution on $\{1 , 2 , \ldots , k \}$ with pmf and cdf given by

and

Let $x_{i} = i$ and $r_{i} = p ( 1 ) + \cdot \cdot \cdot + p ( x_{i} ) = F ( x_{i} ) = i / k {\mathrm{~ f o r ~}} i = 1 , 2 , \dots , k$. Then, from Inequality (13), it can be seen that, if the generated random number R satisfies

then $X$ is generated by setting $X = i$. Now Inequality (14) can be solved for i:

Let y denote the smallest integer $\geq y$. For example,$\left\lceil 7 . 82 \right\rceil = 8 , \left\lceil 5 . 13 \right\rceil = 6$, and $\lceil - 1 . 32 \rceil = - 1$ For $y \ge 0 , \lceil y \rceil$ is a function that rounds up. This notation and Inequality (15) yield a formula for generating X, namely

For example, consider the generating of a random variate X that is uniformly distributed on $\{1 , 2 , \ldots , 10 \}$. The variate X might represent the number of pallets to be loaded onto a truck. Using Table A.1 as a source of random numbers R and using Equation (16) with $k = 10$ yields

The procedure discussed here can be modified to generate a discrete uniform random variate with any range consisting of consecutive integers. Exercise 13 asks the student to devise a procedure for one such case.

Example 6: The Geometric Distribution

Consider the geometric distribution with pmf

where $0 < p < 1$. Its cdf is given by

for $x = 0 , 1 , 2 , . . .$. Using the inverse-transform technique [i.e., Inequality (13)], recall that a geometric random variable X will assume the value x whenever

where R is a generated random number assumed $0 < R < 1$. Solving Inequality (17) for x proceeds as follows:

But $1 - p < 1$ implies that ln $\left( 1 - p \right) < 0$, so that

Thus,$X = x$ for that integer value of x satisfying Inequality (18). In brief, and using the round-up function , we have

Since $p$ is a fixed parameter, let $\beta = - 1 / \ln ( 1 - p )$. Then $\beta > 0$ and, by Equation (19),$X =$$\lceil - \beta \ln ( 1 - R ) - 1 \rceil$. By Equation $( 1 ) , - \beta \ln ( 1 - R )$ is an exponentially distributed random variable with mean $\beta ;$ so one way of generating a geometric variate with parameter $p$ is to generate (by any method) an exponential variate with parameter $\beta^{- 1} = - \ln ( 1 - p )$, subtract one, and round up.

Occasionally, there is needed a geometric variate X that can assume values q $, q + 1 , q + 2 , \ldots \}$ with pm ${\mathrm{\Sigma}} ^ {\prime} p ( x ) = p ( 1 - p ) ^ {x - q} ( x = q , q + 1 , \ldots )$. Such a variate X can be generated, via Equation (19), by

One of the most common cases is $q = 1$

Example 7

Generate three values from a geometric distribution on the range $\{X \ge 1 \}$ with mean 2. Such a geometric distribution has pmf $\begin{array} {r} {\boldsymbol {p} ( \boldsymbol {x} ) = p ( 1 - p ) ^ {x - 1} ( \boldsymbol {x} = 1 , 2 , \dots )} \end{array}$ with mean $1 / p = 2 , \mathrm{o r} p = 1 / 2$ Thus, X can be generated by Equation (20) with $q = 1 , p = 1 / 2$, and $1 / \ln ( 1 - p ) = - 1 . 443$. Using Table $\mathrm{A} . 1 , R_{1} = 0 . 932 , R_{2} = 0 . 105$, and $R_{3} = 0 . 687$ yields

Exercise 15 deals with an application of the geometric distribution.

2 Acceptance-Rejection Technique

Suppose that an analyst needed to devise a method for generating random variates, X, uniformly distributed between 1/4 and 1. One way to proceed would be to follow these steps:

Step 1. Generate a random number R.

Step 2a. If $R \ge 1 / 4_{\cdot}$, accept $X = R$, then go to Step 3.

Step 2b. I $\lceil R < 1 / 4$, reject R, and return to Step 1.

Step 3. If another uniform random variate on [1/4, 1] is needed, repeat the procedure beginning at Step 1. If not, stop.

Each time Step 1 is executed, a new random number R must be generated. Step 2a is an “acceptance” and Step 2b is a “rejection” in this acceptance-rejection technique. To summarize the technique, random variates R with some distribution (here uniform on [0, 1]) are generated until some condition (here $R > 1 / 4 )$ is satisfied. When the condition is finally satisfied, the desired random variate X (here uniform on $[ 1 / 4 , 1 ] )$ can be computed $( X = R )$. This procedure can be shown to be correct by recognizing that the accepted values of R are conditioned values; that is, R itself does not have the desired distribution, but R conditioned on the event $\{R \ge 1 / 4 \}$ does have the desired distribution. To show this, take $1 / 4 \leq a < b \leq 1 ;$; then

which is the correct probability for a uniform distribution on [1/4, 1]. Equation (21) says that the probability distribution of R, given that R is between 1/4 and 1 (all other values of R are thrown out), is the desired distribution. Therefore, if $1 / 4 \le R \le 1$, set $X = R$ The efficiency of an acceptance-rejection technique depends heavily on being able to minimize the number of rejections. In this example, the probability of a rejection is $P ( R < 1 / 4 ) = 1 / 4 ,$, so that the number of rejections is a geometrically distributed random variable with probability of “success” being $p = 3 / 4$ and mean number of rejections $( 1 / p - 1 ) = 4 / 3 - 1 = 1 / 3$. (Example 6 discussed the geometric distribution.) The mean number of random numbers R required to generate one variate X is one more than the number of rejections; hence, it is $4 / 3 = 1 . 33 .$. In other words, to generate 1000 values of X would require approximately 1333 random numbers R.

In the present situation, an alternative procedure exists for generating a uniform variate on [1/4, 1]—namely, Equation (5), which reduces to $X = 1 / 4 + ( 3 / 4 ) R$. Whether the acceptancerejection technique or an alternative procedure, such as the inverse-transform technique [Equation (5)], is the more efficient depends on several considerations. The computer being used, the skills of the programmer, and the relative inefficiency of generating the additional (rejected) random numbers needed by acceptance-rejection should be compared to the computations required by the alternative procedure. In practice, concern with generation efficiency is left to specialists who conduct extensive tests comparing alternative methods (i.e., until a simulation model begins to require excessive computer runtime due to the generator being used).

For the uniform distribution on [1/4, 1], the inverse-transform technique of Equation (5) is undoubtedly much easier to apply and more efficient than the acceptance-rejection technique. The main purpose of this example was to explain and motivate the basic concept of the acceptancerejection technique. However, for some important distributions, such as the normal, gamma and beta, the inverse cdf does not exist in closed form and therefore the inverse-transform technique is difficult. These more advanced techniques are summarized by Bratley, Fox, and Schrage [1996], and Law [2007].

In the following subsections, the acceptance-rejection technique is illustrated for the generation of random variates for the Poisson, nonstationary Poisson, and gamma distributions.

2.1 Poisson Distribution

A Poisson random variable N with mean $\alpha > 0$ has pmf

More important, however, is that N can be interpreted as the number of arrivals from a Poisson arrival process in one unit of time. The interarrival times,$A_{1} , A_{2} , \ldots$. of successive customers are exponentially distributed with rate α (i.e., α is the mean number of arrivals per unit time); in addition, an exponential variate can be generated by Equation (3). Thus, there is a relationship between the (discrete) Poisson distribution and the (continuous) exponential distribution:

if and only if

Equation (22),$N = n ,$, says there were exactly n arrivals during one unit of time; but Inequality (23) says that the nth arrival occurred before time 1 while the $( n + 1 ) \mathrm {{s t}}$ t arrival occurred after time 1. Clearly, these two statements are equivalent. Proceed now by generating exponential interarriva times until some arrival, say $n + 1$, occurs after time 1; then set $N = n$ For efficient generation purposes, Inequality (23) is usually simplified by first using Equation (3),$A_{i} = ( - 1 / \alpha )$ ln R , to obtain

Next multiply through $\boldsymbol {\mathsf {b y}} - \boldsymbol {\alpha}$, which reverses the sign of the inequality, and use the fact that a sum of logarithms is the logarithm of a product, to get

Finally, use the relation $e^{\ln x} = x$ for any number x to obtain

which is equivalent to Inequality (23). The procedure for generating a Poisson random variate, N, is given by the following steps:

Step 1. Set $n = 0 , P = 1$

Step 2. Generate a random number $R_{n + 1}$, and replace P by $P \cdot R_{n + 1}$ Step 3. If $\textit {P} < \textit {e} ^ {- \alpha}$, then accept $N = n$. Otherwise, reject the current n, increase n by one, and return to step 2.

Notice that, upon completion of Step 2, P is equal to the rightmost expression in Inequality (24). The basic idea of a rejection technique is again exhibited; if $P \ge e^{- \alpha}$ in Step 3, then n is rejected and the generation process must proceed through at least one more trial.

How many random numbers will be required, on the average, to generate one Poisson variate N? If N n, then n 1 random numbers are required, so the average number is given by

which is quite large if the mean α of the Poisson distribution is large.

Example 8

Generate three Poisson variates with mean $\alpha = 0 . 2$. First, compute $e^{- \alpha} = e^{- 0 . 2} = 0 . 8187$. Next, get a sequence of random numbers R from Table A.1 and follow the previously described Steps 1 to 3:

Step 1.$\mathrm{S e t} n = 0 , P = 1$ Step 2.$R_{1} = 0 . 4357 , P = 1 \cdot R_{1} = 0 . 4357$ Step 3. Since $P = 0 . 4357 < e^{- \alpha} = 0 . 8187$, accept N 0.

Step 1–3.$( R_{1} = 0 . 4146 \mathrm{l e a d s t o} N = 0 . )$ Step 1. Set $n = 0 , P = 1$ Step 2.$R_{1} = 0 . 8353 , P = 1 \cdot R_{1} = 0 . 8353 .$ Step 3. Since $P \ge e^{- \alpha}$, reject n 0 and return to Step 2 with n 1.

Step 2.$R_{2} = 0 . 9952 , P = R_{1} R_{2} = 0 . 8313 .$ Step 3. Since $P \ge e^{- \alpha}$, reject n 1 and return to Step 2 with n 2.

Step 2.$R_{3} = 0 . 8004 , P = R_{1} R_{2} R_{3} = 0 . 6654 .$ Step 3. Since $P < e^{- \alpha}$, accept N 2.

The calculations required for the generation of these three Poisson random variates are summarized as follows:

It took five random numbers R to generate three Poisson variates here $( N = 0 , N = 0$, and $N = 2 )$ but in the long run, to generate, say, 1000 Poisson variates with mean $\alpha = 0 . 2$, it would require approximately $1000 ( \alpha + 1 )$ or 1200 random numbers.

Example 9

Buses arrive at the bus stop at Peachtree and North Avenue according to a Poisson process with a mean of one bus per 15 minutes. Generate a random variate N which represents the number o arriving buses during a 1-hour time slot.$\operatorname{N o w} ,$ N is Poisson distributed with a mean of four buses per hour. First, compute $e^{- \alpha} = e^{- 4} = 0 . 0183$. Using a sequence of 12 random numbers from Table A.1 yields the following summarized results:

It is immediately seen that a larger value of α (here $\alpha = 4 )$ usually requires more random numbers; if 1000 Poisson variates were desired, approximately $1000 ( \alpha + 1 ) = 5000$ random numbers would be required.

When α is large, say,$\alpha \geq 15$, the rejection technique outlined here becomes quite expensive, but fortunately an approximate technique based on the normal distribution works quite well. When the mean α is large, then

is approximately normally distributed with mean zero and variance 1; this observation suggests an approximate technique. First, generate a standard normal variate Z, by Equation (28) in Section 3.1, then generate the desired Poisson variate N by

where is the round-up function described in Section 1.7.$\operatorname{I f} \alpha + {\sqrt {\alpha}} Z - 0 . 5 < 0$, then set $N = 0 . )$ The $" 0 . 5 '$ used in the formula makes the round-up function become a “round to the nearest integer” function. Equation (25) is not an acceptance-rejection technique, but, when used as an alternative to the acceptance-rejection method, it provides a fairly efficient and accurate method for generating Poisson variates with a large mean.

2.2 Nonstationary Poisson Process

Another type of acceptance-rejection method (which is called thinning) can be used to generate interarrival times from a nonstationary Poisson process (NSPP) with arrival rate $\lambda ( t ) , 0 \leq t \leq T$. An NSPP is an arrival process with an arrival rate that varies with time.

Consider, for instance, the arrival-rate function given in Table 7 that changes every hour. The idea behind thinning is to generate a stationary Poisson arrival process at the fastest rate $( 1 / 5$ customer per minute in the example), but “accept” or admit only a portion of the arrivals, thinning out jus enough to get the desired time-varying rate. Next we give the generic algorithm, which generates $\boldsymbol {\mathcal {T}} _ {i}$ as the time of the ith arrival. Remember that, in a stationary Poisson arrival process, the times between arrivals are exponentially distributed.

Step 1. Let $\begin{array} {l} {{\lambda^{*} = \operatorname* {m a x} _ {0 \leq t \leq T} \lambda ( t )}} \\ {{1 .}} \end{array}$ be the maximum of the arrival rate function and set $t = 0$ and $i = 1$ Step 2. Generate E from the exponential distribution with rate $\lambda^{*}$ and let $t = t + E$(this is the arrival time of the stationary Poisson process).

Step 3. Generate random number R from the $U ( 0 , 1 )$ distribution. If $R \leq \lambda ( t ) / \lambda^{*}$ then $\mathcal {T} _ {i} = t$ and $i = i + 1 .$ Step 4. Go to Step 2.

The thinning algorithm can be inefficient if there are large differences between the typical and the maximum arrival rate. However, thinning has the advantage that it works for any integrable arriva rate function, not just a piecewise-constant function, as in this example.

Table 7 Arrival Rate for NSPP Example

Example 10

For the arrival-rate function in Table 7, generate the first two arrival times.

Step 1.$\begin{array} {r} {\lambda^{*} = \operatorname* {m a x} _ {0 \leq t \leq T} \lambda ( t ) = 1 / 5 , t = 0 \mathrm{a n d} i = 1} \end{array}$ Step 2. For random number R 0.2130, E 5 ln(0.213) 13.13 and t 0 13.13 13.13.

Step 3. Generate $R = 0 . 8830$. Since $R = 0 . 8830 \not \equiv \lambda ( 13 . 13 ) / \lambda^{*} = ( 1 / 15 ) / ( 1 / 5 ) = 1 / 3$, do not generate the arrival.

Step 4. Go to Step 2.

Step 2. For random number $R = 0 . 5530 , E = - 5 \ln ( 0 . 553 ) = 2 . 96 \mathrm{a n d} t = 13 . 13 + 2 . 96 = 16 . 09$ Step 3. Generate $R = 0 . 0240 . \mathrm{~ S i n c e ~} R = 0 . 0240 \leq \lambda ( 16 . 09 ) / \lambda^{*} = ( 1 / 15 ) / ( 1 / 5 ) = 1 / 3 , \mathrm{~ s e t ~}$ t 16.09 and i i 12.

Step 4. Go to Step 2.

Step 2. For random number $R = 0 . 0001 , E = - 5 \ln ( 0 . 0001 ) = 46 . 05 \mathrm{a n d} t = 16 . 09 + 46 . 05 =$62.14.

Step 3. Generate $R = 0 . 1443$. Since $R = 0 . 1443 \leq \lambda ( 62 . 14 ) / \lambda^{*} = ( 1 / 12 ) / ( 1 / 5 ) = 5 / 12$, set t 62.14 and i i 13.

Step 4. Go to Step 2.

2.3 Gamma Distribution

Several acceptance-rejection techniques for generating gamma random variates have been developed. (See Bratley, Fox, and Schrage [1996] and Law [2007].) One ofthe more efficient is by Cheng [1977]; the mean number of trials is between 1.13 and 1.47 for any value of the shape parameter $\beta \geq 1$ If the shape parameter $\beta$ is an integer, say,$\beta = k ,$, one possibility is to use the convolution technique in Example 12, because the Erlang distribution is a special case of the more genera gamma distribution. On the other hand, the acceptance-rejection technique described here would be a highly efficient method for the Erlang distribution, especially if $\beta = k$ were large. The routine generates gamma random variates with scale parameter θ and shape parameter $\beta \mathrm{\cdot}$—that is, with mean $1 / \theta$ and variance $1 / \beta \theta^{2}$ The steps are as follows

Step 1. Compute $a = 1 / ( 2 \beta - 1 ) ^ {1 / 2} , b = \beta - \ln 4$ Step 2. Generate $R_{1}$ and $R_{2}$. Set $V = R_{1} / ( 1 - R_{1} )$ Step 3. Compute $X = \beta V^{a}$ Step 4a. If $\begin{array} {r} {X > b + ( \beta a + 1 ) \ln ( V ) - \ln ( R_{1} ^ {2} R_{2} )} \end{array}$), reject X and return to Step 2.

Step 4b. If $X \le b + ( \beta a + 1 ) \ln ( V ) - \ln ( R_{1} ^ {2} R_{2} )$, use X as the desired variate.

The generated variates from Step 4b will have mean and variance both equal to $\beta .$. If it is desired to have mean $1 / \theta$ and variance $1 / \beta \theta^{2}$, then include Step 5.

(Step 5. Replace X by $X / ( \beta \theta ) .$)

The basic idea of all acceptance-rejection methods is again illustrated here, but the proof of this example is beyond the scope of this text. In Step 3,$X = \beta V^{a} = \beta [ R_{1} / ( 1 - R_{1} ) ] ^ {a}$ is not gamma distributed, but rejection of certain values of X in Step 4a guarantees that the accepted values in Step 4b do have the gamma distribution.

Example 11

Downtimes for a high-production candy-making machine have been found to be gamma distributed with mean 2.2 minutes and variance 2.10 minutes^2. Thus,$1 / \theta = 2 . 2$ and $1 / \beta \theta^{2} = 2 . 10$, which together imply that β 2.30 and θ 0.4545.

Step 1. a 0.53, b 0.91.

Step 2. Generate R 0.832, R 0.021. Set V 0.832/(10.832) 4.952.

Step 3. Compute X 2.3(4.952)^0.53 5.37.

Step 4. X 5.37 > 0.91 [2.3(0.53) 1] ln(4.952) ln[(0.832)^20.021] 8.68, so reject X and return to Step 2.

Step 2. Generate $R_{1} = 0 . 434 , R_{2} = 0 . 716 . \mathrm {{S e t}} V = 0 . 434 / ( 1 - 0 . 434 ) = 0 . 767 .$ Step 3. Compute X 2.3(0.767)^0.53 2.00.

Step 4. Since X 2.000.91 [2.3(0.53) 1] ln(0.767) ln[(0.434)^20.716] 2.32, accept X.

Step 5. Divide X by βθ 1.045 to get X 1.91.

This example took two trials (i.e., one rejection) to generate an acceptable gamma-distributed random variate, but, on the average, to generate, say, 1000 gamma variates, the method will require between 1130 and 1470 trials, or equivalently, between 2260 and 2940 random numbers. The method is somewhat cumbersome for hand calculations, but is easy to program on the computer and is one of the most efficient gamma generators known.

3 Special Properties

“Special properties” are just as the name implies: They are variate-generation techniques that are based on features of a particular family of probability distributions, rather than being general-purpose techniques like the inverse-transform or acceptance-rejection techniques.

3.1 Direct Transformation for the Normal and Lognormal Distributions

Many methods have been developed for generating normally distributed random variates. The inversetransform technique cannot easily be applied, however, because the inverse cdf cannot be written in closed form. The standard normal cdf is given by

This section describes an intuitively appealing direct transformation that produces an independen pair of standard normal variates with mean zero and variance 1. The method is due to Box and Muller [1958]. Although not as efficient as many more modern techniques, it is easy to program in a scientific language, such as FORTRAN, C, C , Visual Basic, or Java. We then show how to transform a standard normal variate into a normal variate with mean $\mu$ and variance $\sigma^{2}$. Once we have a method (this or any other) for generating $X$ from an $N ( \mu , \sigma^{2} )$ distribution, then we can generate a lognorma random variate $Y$ with parameters $\mu$ and $\sigma^{2}$ by using the direct transformation $Y = e^{X}$. (Recall that $\mu$ and $\sigma^{2}$ are not the mean and variance of the lognormal.)

Consider two standard normal random variables,$Z_{1}$ and $Z_{2}$, plotted as a point in the plane as shown in Figure 7 and represented in polar coordinates as

Figure 7 Polar representation of a pair of standard normal variables.

It is known that $B^{2} = Z_{1} ^ {2} + Z_{2} ^ {2}$ has the chi-square distribution with 2 degrees of freedom, which is equivalent to an exponential distribution with mean 2. Thus, the radius B can be generated by use of Equation (3):

By the symmetry of the normal distribution, it seems reasonable to suppose, and indeed it is the case, that the angle is uniformly distributed between 0 and 2π radians. In addition, the radius B and the angle $\theta$ are mutually independent. Combining Equations (26) and (27) gives a direct method for generating two independent standard normal variates.$Z_{1}$ and $Z_{2} .$, from two independent random numbers,$R_{1}$ and $R_{2} {\mathrm{:}}$

To illustrate the generation scheme, consider Equation (28) with $R_{1} = 0 . 1758$ and $R_{2} = 0 . 1489$ Two standard normal random variates are generated as follows: