1.3 Fighting the Exponential Distribution is Futile

Perhaps the mathematical arguments in the preceding section have not fully convinced you that exponential inter-arrival times and Poisson arrival counts arise naturally in many systems. We therefore use a simulation to show how exponential inter-arrival times emerge when many nearly periodic patient streams are combined.

The simulation is based on an internship of one of my master’s students with the task of evaluating the scheduling system for patients with inflammatory bowel disease (IBD) at a hospital. Patients with IBD can be categorized as regular and urgent. Regular patients required check-ups approximately every six months, whereas urgent patients needed to be seen within one week, that is, five working days.¹¹This procedure has been changed recently; however, the ideas of this case apply much more generally. Each consultation takes 10 minutes. Immediately after a consultation, regular patients scheduled their next visit for approximately six months later. There are about 1000 patients in the system and a working year consists of 200 working days. Given that each patient is seen twice per year, a physician needed to consult with an average of 10 patients daily.²²Two visits per year per patient times 1000 patients distributed over 200 days per year is 10 patients per day. The problem was to design a planning procedure that allowed enough flexibility to meet the time constraints for the urgent patients without leading to idle times for the physicians. Of course, new patient intakes occurred, but for simplicity we assume that the patient population remains constant.

Planning regular patients six months in advance seems to lead to stable arrival times, but this is not true, as there are many practical disturbances. A significant number of patients rescheduled their appointment date, for example, due to altered vacation plans. There were no-shows and yet other patients relocated to other cities; some died. However, even without such practical problems and even when each patient individually shows relatively predictable behavior, the simulation results below show that at the population level the inter-arrival times are well modeled as exponentially distributed rvs.

In conclusion, whether the planner plans six months in advance or does not plan at all, there will not be much difference in the variability of the number of patients arriving per day. In other words, the idea that making appointments half a year in advance will reduce variability by a large amount is simply not true. The verdict must then be clear: When rules do not have an effect, it is best to drop them altogether. It is simply better to send patients a reminder for a new appointment three weeks before a new visit is due, rather than planning six months in advance. Far fewer replanning actions are needed, and the probability for no-shows decreases.

We now turn to developing the simulation to show how the exponential distribution arises ’out of thin air’. We focus only on regular patients and neglect urgent patients as this is a relatively small group. Assume that the hospital has just \(50\) regular patients. For ease, scale the time so that each patient plans the next appointment uniformly in the interval \([0.9 , 1.1]\), independent of previous appointments and other patients.³³Like this, the inter-arrival times between two consecutive appointments correspond to half a year (100 working days) plus or minus 10 days (observe that \([0.9, 1.1] \times 100 = [90, 110]\)).

The question of interest is to estimate the distribution of inter-arrival times between two patients as seen by the planner of the hospital. Of course, with only one patient this is straightforward, the planner sees \(\sim \Unif{[0.9, 1.1]}\). But what happens for many patients?

Fig. 1.1 provides the answer to this question, but we first discuss the algorithms used to produce the figure, as they clarify how to interpret it.

The simulation uses \(N=50\) patients and \(T=10\) half-year periods. For each patient \(k\), the inter-arrival times \(X_i^k\sim\Unif{[0.9,1.1]}\) are independent, and the arrival times satisfy \(A_i^k=A_{i-1}^k+X_i^k\), as in (1.2.2). Alg. 1.3.1 applies these recursions to all patients.

\begin{algorithm}
\caption{Generate the arrival times of all patients.}
\begin{algorithmic}[1]
\State \textbf{Input:} number of patients \(N\), number of visits \(T\)
\State \textbf{Output:} arrival times \(\{A_i^k: 1\leq k\leq N,\ 0\leq i\leq T\}\)
\Procedure{GeneratePatientArrivals}{$N,T$}
    \For{\(1\leq k\leq N\)}
        \State \(A_0^k \gets 0\)
        \For{\(1\leq i\leq T\)}
            \State Draw \(X_i^k\sim\Unif{[0.9,1.1]}\)
            \State \(A_i^k\gets A_{i-1}^k+X_i^k\)
        \EndFor
    \EndFor
    \Return \(\{A_i^k: 1\leq k\leq N,\ 0\leq i\leq T\}\)
\EndProcedure
\end{algorithmic}

\end{algorithm}

Next, we build the upper left panel of Fig. 1.1. Each row contains a number of dots, and each dot corresponds to an arrival of a patient.⁴⁴The \(x\) coordinate of the \(i\)th dot of patient \(k\) is \(A_{i}^k\) and its \(y\) coordinate is \(k\). As we have \(50\) patients, we have \(50\) such rows of dots.

We plot the arrivals as seen by the planner as a jitter plot at \(y=0\). A jitter plot adds small random noise in the \(y\) direction to the data points. By spreading out the points slightly, jittering allows for a clearer visualization of the distribution of individual points. Without the jitter, we would just see small line segments instead of somewhat rectangular shapes.

The hospital combines the arrivals of all patients into one ordered sequence. Alg. 1.3.2 sorts the individual arrival times and computes the time between each pair of consecutive arrivals as seen by the hospital.

Using Alg. 1.3.2, we compute the inter-arrival times seen by the hospital and assemble them in 50 bins. On average, we have \(\lambda=50\) arrivals per unit time; thus, the mean inter-arrival time is \(1/\lambda = 1/50\). Next, recall that the pdf of the exponential distribution is \(\lambda e^{-\lambda x}\). To plot the pdf, we take as support for the plot \(5\) times the mean time between arrivals, that is, \(5/\lambda = 0.1\). We then plot the histogram of the inter-arrival times and the pdf in the upper right panel.

\begin{algorithm}
\caption{Merge patient arrivals and compute hospital inter-arrival times.}
\begin{algorithmic}[1]
\State \textbf{Input:} patient arrival times \(\{A_i^k:1\leq k\leq N,\ 1\leq i\leq T\}\)
\State \textbf{Output:} inter-arrival times \((X_2,\ldots,X_{NT})\)
\Procedure{MergePatientArrivals}{$\{A_i^k\}$}
    \State \((A_1,\ldots,A_{NT})\gets\operatorname{sort}\{A_i^k:1\leq k\leq N,\ 1\leq i\leq T\}\)
    \For{\(2\leq j\leq NT\)}
        \State \(X_j\gets A_j-A_{j-1}\)
    \EndFor
    \Return \((X_2,\ldots,X_{NT})\)
\EndProcedure
\end{algorithmic}

\end{algorithm}

What can we see in the upper panels of Fig. 1.1? In the upper left panel, the arrival times of each patient are at first tightly concentrated around \(1\), but as time progresses, the arrival times become more spread out. The jitter plot around \(y=0\) also demonstrates this. The rectangles become wider until they almost overlap. In the right upper panel, we see that in comparison to the exponential density, there are relatively many short inter-arrival times. This is expected because all patients start at the same time, so arrival times will be concentrated for the first couple of visits.⁵⁵But even with these unrealistic starting times, the exponential distribution is not a bad fit. Therefore, letting all patients start at the same time is not a good model for a hospital that has been serving patients for years.

To compensate for these unnatural start times, we set the initial arrival time of each patient as uniformly distributed on \([0,1]\). With this change, we run the above code again and plot the results in the lower left and right panels of Fig. 1.1. Now we see that the jitter rectangles, that is, the arrivals as seen by the hospital, overlap right from the start, and the exponential distribution is an excellent fit.

To conclude, if many patients arrive at the hospital and each contributes some small variation to the inter-arrival times, the hospital sees exponential inter-arrival times for the population. Clearly, this is not specific to hospitals: It applies to any system that serves a large pool of customers, each arriving with some variation. To avoid any confusion, here we model the inter-arrival times for one patient as uniformly distributed, but the same phenomenon occurs when the inter-arrival times for one patient have other types of distribution.

In the above we sorted the arrivals of the patients. Alg. 1.3.3 shows a recursive⁶⁶And very efficient. algorithm to carry out the sorting.⁷⁷Think deeply about this algorithm. It is very useful. We use the following sorting algorithm below. Is uses the procedure \(\pr{ChoosePivot}\) to select an element uniformly from its non-empty input list.

\begin{algorithm}
\caption{Sort a list \(\Sset\) in increasing order.}
\begin{algorithmic}[1]
\State \textbf{Input:} a list $\Sset$.
\State \textbf{Output:} the elements of $\Sset$ in increasing order.
\Procedure{QuickSort}{$\Sset$}
    \If{$|\Sset| \leq 1$}
        \Return $\Sset$
    \EndIf
    \State $p \gets \pr{ChoosePivot}(\Sset)$
    \State $\Lset \gets \qb{x \in \Sset : x < p}$
    \State $\Mset \gets \qb{x \in \Sset : x = p}$
    \State $\Rset \gets \qb{x \in \Sset : x > p}$
    \Return $\pr{QuickSort}(\Lset) + \Mset + \pr{QuickSort}(\Rset)$
\EndProcedure
\end{algorithmic}

\end{algorithm}

\begin{algorithm}
\begin{algorithmic}[1]
\State \textbf{Input:} arrival times \(\{A_i^k:1\leq k\leq N,\ 1\leq i\leq T\}\)
\State \textbf{Output:} inter-arrival times \((X_2,\ldots,X_{NT})\)
\Procedure{MergeArrivals}{$\{A_i^k\}$}
    \State \((A_1,\ldots,A_{NT})\gets\operatorname{sort}\{A_i^k:1\leq k\leq N,\ 1\leq i\leq T\}\)
    \For{\(2\leq j\leq NT\)}
        \State \(X_j\gets A_j-A_{j-1}\)
    \EndFor
    \Return \(\operatorname{sort}(X_2,\ldots,X_{NT})\)
\EndProcedure
\end{algorithmic}

\end{algorithm}

\begin{algorithm}
\begin{algorithmic}[1]
\State \textbf{Input:} a list $\Sset$.
\State \textbf{Output:} the elements of $\Sset$ in increasing order.
\Procedure{QuickSort}{$\Sset$}
    \If{$|\Sset| \leq 1$}
        \Return $\Sset$
    \EndIf
    \State $p \gets \pr{ChoosePivot}(\Sset)$
    \State $\Lset \gets \text{???}$
    \State $\Mset \gets \qb{x \in \Sset : x = p}$
    \State $\Rset \gets \qb{x \in \Sset : x > p}$
    \Return $\pr{QuickSort}(\Lset) + \Mset + \pr{QuickSort}(\Rset)$
\EndProcedure
\end{algorithmic}

\end{algorithm}

\begin{algorithm}
\begin{algorithmic}[1]
\State \textbf{Input:} a list $\Sset$.
\State \textbf{Output:} the elements of $\Sset$ in increasing order.
\Procedure{QuickSort}{$\Sset$}
    \If{$|\Sset| \leq 1$}
        \Return $\Sset$
    \EndIf
    \State $p \gets \pr{ChoosePivot}(\Sset)$
    \State $\Lset \gets \qb{x \in \Sset : x < p}$
    \State $\Mset \gets \qb{x \in \Sset : x = p}$
    \State $\Rset \gets \qb{x \in \Sset : x > p}$
    \State \textbf{return} ???
\EndProcedure
\end{algorithmic}

\end{algorithm}

The next two questions provide some further practice with the keys and values of maps.

\begin{algorithm}
\caption{Store a probability mass function in a map.}
\begin{algorithmic}[1]
\State \(p\gets\pr{Map}()\)
\State \(p(0)\gets1/4\)
\State \(p(1)\gets1/2\)
\State \(p(2)\gets1/4\)
\end{algorithmic}

\end{algorithm}

\begin{algorithm}
\caption{Update values stored in a map.}
\begin{algorithmic}[1]
\State \(M\gets\pr{Map}()\)
\State \(M(2)\gets5\)
\State \(M(5)\gets M(5)+2\)
\State \(M(2)\gets M(2)+1\)
\end{algorithmic}

\end{algorithm}

It remains to formalize the observations from the above plots.⁸⁸If you are not interested in probability theory, you can skip the rest of this section. In fact, we should distinguish between two limiting regimes. Patients recurrently plan appointments, separated around some mean period, that is, roughly every half a year. The upper left panel in Fig. 1.1 suggests that if the number of appointments increases, the appointments for each patient eventually become uniformly distributed on the period length (here half a year). The lower left panel suggests that if the appointments are uniformly distributed to begin with, they remain uniformly distributed. So, our first step is to prove the claim that in the limit \(t\to \infty\), stochastically recurring inter-arrival times become uniformly distributed over an interval. The second claim is that, in the limit of many patients, each arriving at a uniformly distributed time in an interval, the inter-arrival times of all patients together become exponentially distributed, in a scaled sense.

The theorems and proofs are elegant and combine order statistics, the beta distribution, first-step analysis, and some infinitesimal calculus to give the reasoning a heuristic flavor. Finally, we discuss an important smoothing property of taking expectations. This idea is also relevant in data science when moving averages are used. Let us first show that periodically planned arrivals, with some stochasticity, tend to become uniformly distributed over an interval. To formalize this claim, we consider a discrete-time model. There are on average \(C\) days between two successive visits of a patient (for the moment, the time on a day is of no importance). For instance, let the inter-arrival times \(X_{i}\) be such that \(\P{X_{i}=C-1} = \P{X_{i}=C+1}=1/2\). In this case, if the process starts on day \(0\) and \(X_{1} = C-1\), the next visit will be on day \(C-1\), but if \(X_{1}=C+1\), the next visit will be on day \(1\) of the next cycle. Looking at the arrival process like this, the patient moves on a circle with states \(\{0, 1, \ldots, C-1\}\), and with every visit the patient moves with probability \(1/2\) one state to the left or to the right. Thus, when \(X=C-1\), we can just as well model this as \(X=-1\).⁹⁹There is a small subtlety when the cycle length is even and \(\P{X=C-1} = \P{X=C+1} = 1/2\). As an even number is a multiple of \(2\), the odd (even) days can only be visited after an odd (even) number of visits. The problem is easily resolved by assuming, in addition, that \(C\) is odd or \(\P{X = 0} > 0\).

Applied to the hospital, after many visits of a patient—each roughly half a year apart— the visits will appear to be uniformly distributed on a `circle’ of half-year length.

There is one point that remains about the lower left panel of Fig. 1.1; recall that for this panel we started the simulation with uniformly distributed arrival times. What is the distribution of the arrival times on the circle if we start like this? Let us analyze this for the case in which the inter-arrival times are uniform on \([0,1]\) and we take the circle as \(C=[0, 1)\). Observe that when \(C=[0, 1)\), \(U+V \bmod 1\), that is, the sum of two iid uniform rvs on the circle, is just the fractional part of \(U + V\).

We now deal with the case with many patients, that is, when the number \(n\) of patients goes to infinity. By the above results, we model the arrival times as iid uniform rvs \(\{U_k\}_{k=1}^{n}\) such that \(U_{k}\sim \Unif{[0,1)}\). Recall that in the simulation we sorted the arrival times of all patients together to obtain the ordered sequence of arrivals as seen by the hospital. In probabilistic terms, the hospital sees the order statistic of \(\{U_k\}_{i=1}^{n}\). Write this as \(0=A^{n}_{0} <A^{n}_1 < A^n_2 < \cdots < A_n^{n}< A_{n+1}^{n} = 1\).¹³¹³We neglect the probability of simultaneous arrivals; these have probability zero. Define the size of the inter-arrivals as seen by the hospital as \(X_{k}^{n} = A_k^n - A_{k-1}^{n}\), \(k=1, \ldots, n+1\).

As the population size \(n\) increases, the size of the first gap, i.e., the properly scaled arrival time of the first patient, becomes increasingly like an exponentially distributed random variable.¹⁴¹⁴We need to scale with a factor \(n\) because the number of patients increases.

Next, we show that all the hospital inter-arrival times have the same density.

Combining the above two lemmas leads straightaway to the next theorem.