# Methods of Generating Random Variates

Techniques for generating random variates are presented. Included among the techniques is the generation of non-homogeneous Poisson processes.

## Distribution Function Inversion

The cumulative distribution function, , for the random variable, , is the probability that the value of a random variable will be less than or equal to the function argument,

Discrete valued random variables can be generated by partitioning the interval between zero and one into sections whose widths are equal to the probability of each value. The value of the random variable is determined by the interval in which a pseudo-random number, U, falls. The probability of U falling in an interval is equal to the width of the interval, which, in turn, is the probability of the corresponding value occurring. This is equivalent to inverting the cumulative distribution function as illustrated in the following figure.

The same technique can be applied to continuous valued random variables. For example, the cumulative distribution for an exponential random variable is,

So, the inverse distribution function of a uniform random number will generate an exponential variate as

When the cumulative distribution function is easily inverted, this technique is recommended. Unfortunately, not very many of the more commonly used probability distributions have easily inverted cumulative distribution functions. Exceptions that are particularly useful are order statistics from a sample of uniform random variables.

Order Statistics: Order statistics are sorted samples; the minimum order statistic is the smallest value in a sample. Suppose we want to know the first time that one of several identical independent components will fail. We could generate the lifetimes of each component and sort them to find the shortest. If there are many components, it would be easier to generate the minimum order statistic from the lifetime distribution.

Order statistics can be generated by evaluating an inverse cumulative distribution function at the corresponding uniform order statistic. The ith smallest of K uniform random variables has a beta distribution with parameters i and K-i+1. A common special case is where we want the smallest of K independent values of a random variable. The cumulative distribution function of the smallest of K independent uniform random variables is,

Therefore,

.

If the lifetime of each of K independent components has an exponential distribution, the distribution of the time until the first component fails is equal to

## Other Methods

Other methods for generating random variates include acceptance/rejection, composition, and special relationships. Like inversion, these other methods are not always possible to use in their pure form, and special algorithms that combine these methods have been invented. The reference by Devroye (1986) contains many of these algorithms.

Acceptance/rejection involves bounding the probability density and generating a point uniformly within this bounding area. If the generated point falls below the density function, the horizontal value of the generated point is used as the random variate. If the point falls outside the density function, the generated point is rejected and a new point is tried. This is continued until an accepted point that falls within the region under the density function is found. This is analogous to throwing uniform points onto the bounding region and accepting those points that fall below the probability density. The algorithm for generating beta variates, beta_fn, in the SIGMA C library uses an acceptance/rejection algorithm due to Cheng (1978).

Composition involves breaking the probability distribution into regions from which it is relatively easy to generate variates. One of these regions is selected with the probability in that region and the variate is generated from that region.

Special properties exploit relationships between different types of random variables. Some commonly used relationships include generating an Erlang variate as a sum of exponentials, generating a geometric variate as the integer part of an exponential, generating a chi-square variate as the sum of squared normal variates, and generating a Poisson as the count of exponential variates whose sum falls within an interval. (When the Poisson rate is large, there are better methods for generating Poisson variates given in the references).

## Generating Non-Homogeneous Poisson Processes

Random arrivals to a service system can often be modeled using a Poisson process. This process has proven to be a rather good model for many processes that occur naturally. It is also used for other types of exogenous events that drive a simulation model, such as equipment failures or flaws occurring in a piece of fabric. The parameter for a Poisson process is its rate, λ, expressed in the number of events per time unit (e.g., customer-arrivals/hour or flaws/square-foot).

There are numerous methods for generating Poisson processes which exploit different properties of these processes. The fact that the times between Poisson events have exponential distributions can be used to simply make the delay time for a self- scheduling edge have an exponential distribution with a mean equal to the inverse of the Poisson rate. For example, if the arrivals to a queue are going to be modeled as a Poisson process with an rate equal to RATE, then the edge delay time for the self-scheduling edge that models sequential arrivals would have a delay time of

```(1.0/RATE)*ERL{1} (recall that Erlang-1 and exponential distributions are same).
```

If we know the total number of Poisson events that occur in an interval of time, we can better match the process by conditioning on this knowledge. We do this using the fact that the distribution of the times of K Poisson events in an interval have the same distribution as K uniform random numbers on the same interval. Say we know that K Poisson events occurred in an interval between 0 and T. We can generate the minimum order statistic of K uniforms over the interval for our first event time, T1. The second event time will be the minimum of the remaining K-1 uniforms distributed over the remaining interval between T1 and T. The rest of the K events in the original interval can be generated in this manner.

Unfortunately, processes in the real world do not often occur with a constant rate. The arrival rate may be relatively high during a rush hour and slack (or zero) late at night. A Poisson process with a changing rate is called a non-homogeneous Poisson process. An easy way to model a non-homogeneous Poisson process is by a technique called "thinning" (Lewis and Shedler, 1979). Here we simply generate Poisson events at the maximum rate and keep them with a probability proportional to the current rate. To illustrate: assume we know that during an 24-hour day, customers will arrive at our carwash with the following hourly rates: R, R,..., R (the rate can be zero if the facility is closed). Let RMAX be the maximum of these rates. The event graph to generate Poisson arrival events according to this daily demand cycle is called NONPOIS.MOD (it is shown in the following figure.