SIMULATION NON-STATIONARY ARRIVAL PROCESS

16th, April, 2019.

Input themselves are generally not of interest and the outputs that the inputs imply are usually of greater interest. It is apparent that inaccurate judgement about input modelling can lead to more error. Since fitting distributions to available data using statistical framework as maximum likelihood estimation is required to understand their inference, our post will focus more on explaining how to avoid errors in input modelling in some certain cases.
Let see one example that we usually face when analysing the queueing in the hospital reception system. There is an electronic server with a touch screen replacing the human receptionist. In practical, there are some kinds of customers who are less comfortable using the electronic receptionist, then the service times for each inviduals may be more variable or slower. The aim is to evaluate the additional delay that this situations might cause. Hospital record data are available to calculate the overall arrival rate. Since customers can make independent decision when they arrive and the exponentially distributed interarrival times are available for some simulation languages, the arrival process of a large number of potential customers are usually assumed by Poisson Process. This is the input model and to fit it, the arrival rate was usually estimated by the total number of arrivals over a certain time interval over the length of this interval.
Note that this arrival process has been treated as stationary Poisson process, i.e. constant arrival rate, which is really strong assumption. There are some obvious reasons that we can reject this assumptions.

  • The distribution of the data sometimes can be different from Poisson. For instance, when a reception only serves scheduled patients and a fixed number of appointments are planned each day, then the arrival process are nearly deterministic and moreover scheduled patients typically arrive not randomly, therefore not Poisson process.
  • The distribution of data can be different from day to day. What happens when the arrival rate changes over time? For example, Saturdays and Sunday are usually busier than weekdays, there could be nonstationary across weeks.
This is a very simple example to emphasize that the input modelling should be examined more carefully, since the aim of input is to get the relevant inference about the outputs. I am not trying to convince you to search for “the exact distribution” to fit the data, but at least, we should approximate the input processes by the available data as best as we can.
In this post, we will present two algorithm to simulate the non-stationary arrival process that can overcome some limitations that we mentioned above. Examples of nonstationary arrival processes that we meet every day are the arrival of fans to football stadium, the arrivals of e-mails to a mail serve, etc. Even though we can define the joint distribution of all the arrival times, it is more natural and convenient for simulation discrete interarrival times or the times between one arrival to the next one.
The idea of methods for simulation is to transform a basic and easily controllable input process to the input process we want. In more details, in order to obtain the nonstationary arrival process, we transform the independent and identically distributed interarrival times to the arrival rate we want while retaining a measure of variability. This makes us think about renewal process, which is the basic form of random process. Let \(N(t)\) be the arrival counting process (number of arrivals by time \( t\)) of a nonstationary arrival process. We define the integrated rate function be the expected number of arrivals by time \(t\) $$ \Lambda(t) = E(N(t))$$ and the arrival rate is define by differentiating the integrated rate function with respect to \(t\), which can be differentiable $$ \lambda(t) = \frac{d}{dt} \Lambda(t)$$ And for the equilibrium renewal process, there is a relationship between them $$ \Lambda(t) = \lambda(t) t$$
Inversion algorithm using \(\Lambda(t)\)

In Inversion a stationary Poisson process with rate one is used as a base, which is still a random process but with a mean interarrival times of \(1\). In order to simulate the sample of interarrival times:
  • Consider the time of the \(n\)-th arrival in the stationary process and evaluate this at \( \Lambda^{-1}(t) \).This is the time of arrival in the non-stationary process.
  • The previous arrival time is subtracted from this new arrival time, generating the \(n\)-th interarrival time.
  • This is done iteratively to simulate a non-stationary arrival process from a stationary one.

centered image

Fig.1 - llustration of inversion method when \(\Lambda(t)=t^2\).
Figure \(1\) illustrates the inversion algorithm for an arrival process with arrival rate \(\lambda(t)=2t\) and therefore \(\Lambda(t)=t^2\). Since the rate is an increasing function of time, arrivals are expected to be more and more dense as time increases. In the vertical axis, the circle points are the arrival times in a rate one base process. We transformed these times by \(\Lambda^{-1}(t)=\sqrt{t}\) into arrival times in the nonstationary process on the horizontal axis.
In order to see why \(\Lambda^{-1}(t)\) indeed maps the arrival times of an equilibrium renewal process with rate \(1\) into the desired nonstationary arrival process with rate \( \lambda(t)\), please find the very short proof in the references.

Thinning algorithm using \(\lambda(t)\)

In Thinning a sample of interarrival times can be generated by:
  • Let \(\tilde{\lambda} = \max_t \lambda(t)\) and generate a random number \(U \sim U(0,1)\).
  • Evaluate the rate function at the time of the \(n^{\text{th}}\) arrival in the maximum rate stationary process, and then calculate the ratio \(\lambda(t)/\tilde{\lambda}\).
  • If \(U \leq \lambda(t)/\tilde{\lambda}\) then we add this interarrival time, else the timer continues.

centered image

Fig.2. - Illustration of the thinning method when \(\lambda(t) = 6 + 4 sin(t)\).
In order to see why thinning version indeed maps the arrival times of an equilibrium renewal process with rate \( \max_{t} \lambda(t)\) into the desired nonstationary arrival process with rate \( \lambda(t)\), please find the proof in the references.
Even though both algorithm lead to arrival processes with the same desired arrival rate which changes over time, they generally do not give exactly the same processes in term of probability such as variance of the number of arrivals by time \(t\), except for the Poisson arrival process that are probabilistically the same.Here is the comparison of two algorithm by highlighting the advantage and drawback of each methods:

centered image

Table1. –Comparison of Inversion and Thinning method.
Indeed we can simulate the desired nonstationary arrival process from the arrival rate or the integrated rate function, which are unknown, but they can be estimated from the data. For example, the integrated rate function which is the expected number of arrivals by time t can be calculated by taking the average number of arrivals by time \(t\) over \(k\) realizations. Note that the integrated rate function are now piecewise constant function which is not convenient for the inversion algorithm. Fortunately, this issue can be tackled by using linearly interpolate between observed arrival times to get strictly increasing integrated rate function. For the arrival rate, it can be very straightforward by computing the average number of arrivals during a sufficiently small interval over the k realizations, then divide by the length of this interval.
I hope this post can help you understand the importance of input modelling as well as how to simulate the nonstationary arrival process with the desired time-varying arrival rate.

Reference:
1. Barry Nelson, Foundations and Methods of Stochastic Simulation: A first course, 2013.


Comments Please: