Theory Thursday- Simulation of a Poisson Process

You are using discrete-event simulation to analyze a process or system. Imagine that your arrivals occur according to a Poisson Process. Without using specialized simulation software, this post shows how to code up a Poisson arrival process. $Y_n$ will represent the n’th arrival time.

Stationary arrivals:
If your arrival rate does not vary over time, then your task is easy. You have two main options:
1. Generate exponential random variables $X_1$ , $X_2$ , … , representing interarrival times. Set $Y_1 = X_1$ , $Y_2 = Y_1+X_2$ , $Y_3=Y_2+X_3$ ,… Stop generating random variables when $Y_m>T$ for some m, where T is your time interval length that you want to simulate.
2. Generate a Poisson random variable $N(T)$ , representing the number of arrivals over a time interval of length T. Then generate $N(T)$ uniform [0,T] random variables, representing arrival times. Sort the arrival times in ascending order to obtain $Y_1$ , $Y_2$ , …

Non-stationary arrivals:
If your arrival rate varies over time, you’ll need an extra step or two of logic to generate your arrival times. Here are two options for simulation:
1. Find the highest arrival rate in your time interval, $\lambda_{max}$ . Generate a stationary Poisson process with rate $\lambda_{max}$ as in option 1 in the stationary arrival section above. For each arrival simulated, generate a uniform[0,1] random variable. If this uniform random variable is less than $\lambda(t)/\lambda_{max}$ , where $\lambda(t)$ is the arrival rate at the generated arrival time, accept the arrival as “real” and keep it. If not, discard the arrival. The “real” arrivals make up an accurate non-stationary Poisson arrival process.
2. Divide the time period [0,T] into small time increments. For each time increment, generate a uniform[0,1] random variable. For each time increment, if the uniform variable is less than or equal to $\lambda(t) dt$ , where $dt$ is the time increment size, then an arrival occurs during that increment. Assign the arrival time to be a random time in the increment. This method is only approximately correct, but it is good enough in most cases and may be faster to simulate in certain cases.