#+BEGIN_COMMENT
.. title: Simulating Arrival Times for a Non-homogeneous Poisson Process
.. slug: nhpp
.. date: 2024-01-01 12:09:31 UTC+01:00
.. tags:
.. category:
.. link:
.. description:
.. type: text
.. has_math: true
#+END_COMMENT
#+title: Simulating Arrival Times for a Non-homogeneous Poisson Process
#+author: Nicky
#+setupfile: ../preamble.org
#+PROPERTY: header-args:python :session :eval no-export :exports both
#+STARTUP: indent
#+STARTUP: showall
* COMMENT Test the page
#+begin_src shell
/usr/bin/emacs --batch -l ../plugins/orgmode/init.el --eval '(nikola-html-export "nhpp.org" "../output/output.html")'
#+end_src
* Intro
We provide an algorithm to generate arrival times of a non-homogeneous Poisson process (nhpp) when the rate function $\lambda(s)$ is piecewise constant.
Specifically, we assume given two sequences $\{\lambda_{j}\}_{j=0}^{\infty}$ and $\{s_{j}\}_{j=0}^{\infty}$, with $s_{0}=0$, the $s_{j}$ strictly increasing, and $s_{k}\to \infty$ as $k\to\infty$, and then we take $\lambda(s) = \lambda_{j}$ for $s\in [s_{j}, s_{j+1}]$.
An example would be a daily scheme of arrival rates of patients per hour at a hospital, where $\lambda_{0} = 1$ from $[0, 8)$, $\lambda_{1} = 10$ from $[8, 12)$, $\lambda_{2} = 8$ from $[12, 18)$ and $\lambda_{4}=3$ from $[18, 24)]$, update $\{s_{j}\}$ accordingly).
* The Maths
One way to simulate a non-homogenous Poisson process with piece-wise constant rates is as follows.
1. For each interval $[s_j, s_{j+1}]$ compute the duration $t_j = s_{j+1}-s_j$.
2. Generate a random deviate $R_j \sim \Pois{\lambda_j t_j}$.
3. Generate $R_j$ uniform deviates $U_i \sim \Unif{[s_j, s_{j+1}]}$.
4. The uniform deviates form the arrivals times on $[s_j, s_{j+1}]$.
This procedure is easy, but the computation of a Poisson distributed random deviate is somewhat costly.
Perhaps the next procedure is more efficient--I have not tested it, so I don't know.
The paper [[https://web.ics.purdue.edu/~pasupath/PAPERS/2011pasB.pdf][Generating Nonhomogeneous Poisson Processes]] by Raghu Pasupathy explains that we need to find the smallest $t$ such that \(\Lambda (s, t) = \int_{s}^{t} \lambda(y)\d y =x \), where $x = -\log u$ and $u\sim\Unif{[0.1]}$.
** A constant rate throughout
If the rate $\lambda(x) \equiv \lambda$ is contant throughout, finding such \(t\) is simple because $\Lambda(s, t) = \lambda (t-s)$ so that \(t= s + s/\lambda\).
This gives the following algorithm:
1. $s = 0$.
2. \(\mineq \log(u)/\lambda, \quad u \sim \Unif{[0,1]}\).
3. Yield $s$ and return to step 2.
** Piecewise Constant
Let us generalize this proceduce so that we can deal with piecewise constant $\lambda(\cdot)$.
Noting that the sequence $\{s_{j}\}$ places a grid on $\R^{+}$, we define
\begin{align*}
j &= \max\{i: s_{i} \leq s\}, & k &=\min\left\{i : \Lambda(s, s_{i}) >x\right\},
\end{align*}
i.e., $s_{j}$ is the point on the grid just below $s$, and $s_{k}$ is the smallest point on the grid such that the integral $\Lambda(s, s_{k})$ exceeds $x$.
For the algorithm, if $s+x/\lambda_{j}< s_{j+1}$, then we can set $s += x/\lambda_{j}$, and return $s$ as the next arrival time.
However, if if $s+x/\lambda_{j}\geq s_{j+1}$, we use the following idea:
\begin{equation*}
x = \Lambda(s, t) = \Lambda(s_{j}, s_{k}) - \Lambda(s_{j}, s) - \Lambda(t, s_{k}).
\end{equation*}
Since $\lambda(\cdot)$ is $\lambda_{j}$ on $[s_{j}, s]$ and $\lambda(\cdot) = \lambda_{k-1}$ on $[t, s_{k})$,
\begin{align*}
\Lambda(s_{j}, s) &= \lambda_{j}(s-s_{k}), & \Lambda(t, s_{k}) = \lambda_{k-1} (s_{k}-t).
\end{align*}
Thus, if we have $s_{j}$ and $s_{k}$, the next arrival time $t$ after $s$ is easily solved from the equation
\begin{equation*}
\lambda_{k-1}(s_{k} -t ) = x + \Lambda(s_{j}, s) - \Lambda(s_{j}, s_{k}).
\end{equation*}
That is,
\begin{equation}\label{eq:nh8}
t = s_{k } - (\Lambda(s_{j}, s_{k}) - x - \Lambda(s_{j}, s) )/\lambda_{k-1}.
\end{equation}
* The Algorithm
With this idea, we obtain the following algorithm. Assume that $\lambda_{0} > 0$ (Otherwise search for the first $\lambda_{j}>0$, and start from there).
1. $s=0$, $j=0$
2. $s \mineq \log(u)/\lambda_{j}$
3. If $s < s_{j+1}$, yield $s$ and return to step 2, else
4. $T = \lambda_{j}(s_{j+1}-s)$.
5. While $T<0$: $j \pluseq 1$, $T \pluseq \lambda_{j}(s_{j+1}-s_{j})$.
6. $s = s_{j+1} - T/\lambda_{j}$.
7. Yield $s$ and return to step 2.
* The Python ```
Iterator[float]:
rate = next(rates)
down, up = next(times), next(times)
while np.isclose(rate, 0):
rate = next(rates)
down, up = up, next(times)
now = max(start, down)
while True:
now += gen.exponential() / rate
T = rate * (up - now)
while T < 0:
rate = next(rates)
down, up = up, next(times)
T += rate * (up - down)
now = up - T / rate
if now > end:
break
yield now
#+end_src
* The Tests
#+begin_src python
from datetime import datetime, timedelta
from itertools import cycle, accumulate
from collections import Counter
import numpy as np
from nhpp import nhpp_range
def test_floats():
gen = np.random.default_rng(3)
start = 0
end = 40
periods = [1, 1, 3]
rates = [0, 1, 10]
times = []
counter = Counter()
counter[0] = 0
for time in nhpp_range(
start, end, accumulate(cycle(periods), initial=start), cycle(rates)
):
times.append(time)
counter[int(time) % 5] += 1
print(counter)
A = np.array(times)
X = A[1:] - A[:-1]
print(X.mean(), X.std())
theoretical_mean_rate = sum(p * r for r, p in zip(rates, periods)) / sum(
periods
)
print(1 / theoretical_mean_rate)
def test_datetimes():
start = datetime(2023, 8, 9, 10, 0)
end = datetime(2023, 8, 10, 18, 0)
rates = list(range(13))
periods = [2] * len(rates)
time = start
for delta in nhpp_range(
0,
(end - start).total_seconds() / 3600,
accumulate(cycle(periods), initial=0),
cycle(rates),
):
time += timedelta(hours=delta)
print(time)
print(start, time)
if __name__ == "__main__":
test_floats()
test_datetimes()
#+end_src
```