Thinning Algorithms for Nonhomogeneous Poisson Processes

Why Thin a Poisson Process?

When an event rate changes over time, drawing successive waiting times from one exponential distribution does not in general simulate the desired process. Thinning provides another exact method: first generate candidate events from a homogeneous Poisson process with a constant rate that is at least as large as the target rate, then retain each candidate with a time-dependent probability. The retained events occur more often where the target intensity is high and less often where it is low.

Work on a fixed interval \([0,T]\), where \(T>0\), on a probability space \((\Omega,\mathcal{F},\mathbb{P})\). Let \(\lambda:[0,T]\to[0,\infty)\) be a bounded measurable function. A counting process \((N_t)_{0\le t\le T}\), with state space \(\{0,1,2,\ldots\}\) and continuous time index, is a nonhomogeneous Poisson process with intensity \(\lambda\) if \(N_0=0\), counts on disjoint time intervals are independent, and for \(0\le a<b\le T\), \(N_b-N_a\) has a Poisson distribution with mean \(\int_a^b\lambda(s)\,ds\). Restricting the times and intervals to \([0,T]\) ensures that every intensity value and integral in this definition is defined.

Definition (Thinning algorithm): Choose a finite constant \(M>0\) such that \(\lambda(t)\le M\) for every \(t\in[0,T]\). Generate candidate times from a homogeneous Poisson process of rate \(M\) on \([0,T]\). Independently generate a Uniform\([0,1]\) mark \(U\) for each candidate at time \(t\). Retain that candidate if \(U\le \lambda(t)/M\), and otherwise discard it. If \(M=0\) is a valid bound, then \(\lambda\) is zero throughout the interval and the algorithm returns no events.

The bound must hold at every time where the algorithm evaluates the intensity. A convenient choice is any known upper bound, not necessarily the smallest one. The method requires that the intensity be evaluable and that a candidate process of rate \(M\) be simulatable. The probability of acceptance is \(\lambda(t)/M\), which is between zero and one because of the bound. Using \(\le\) rather than \(<\) makes no difference to the distribution: a uniform mark equals a fixed threshold with probability zero.

To generate the candidates, draw independent exponential waiting times with rate \(M\), accumulating them until the next candidate would exceed \(T\). Attach independent uniform marks to the candidates. This produces the homogeneous Poisson process and independent marks used below. The exactness argument concerns the distribution of the resulting process, not a guarantee that any one simulated path will resemble a typical path.

Why the Accepted Events Have the Right Law

Theorem (Exactness of thinning): Suppose \(T>0\), \(\lambda:[0,T]\to[0,\infty)\) is bounded and measurable, and \(0<M<\infty\) satisfies \(\lambda(t)\le M\) for every \(t\in[0,T]\). The thinning algorithm produces a nonhomogeneous Poisson process with intensity \(\lambda\) on \([0,T]\).

Proof. Fix an interval \(I=(a,b]\subseteq[0,T]\), and write \(|I|=b-a>0\). Let \(K_I\) be the number of candidate events in \(I\). Since candidates form a homogeneous Poisson process of rate \(M\), \(K_I\) has a Poisson distribution with mean \(M|I|\). Conditional on \(K_I=k\), the \(k\) candidate locations, when considered without ordering, are independent Uniform\((a,b]\) random variables. This is the conditional-uniform property of Poisson arrival times. The independent uniform marks attached to them make the location-mark pairs independent conditional on \(K_I=k\).

For any one such pair \((V,U)\), the probability of acceptance is

p_I=\mathbb{P}\!\left(U\le\frac{\lambda(V)}{M}\right) =\frac{1}{|I|}\int_a^b\frac{\lambda(s)}{M}\,ds =\frac{1}{M|I|}\int_a^b\lambda(s)\,ds.

Thus, conditional on \(K_I=k\), the accepted count \(A_I=N_b-N_a\) has a binomial distribution with \(k\) trials and success probability \(p_I\). For \(0\le z\le1\), its probability generating function is

\begin{aligned} \mathbb{E}[z^{A_I}] &=\mathbb{E}\!\left[(1-p_I+p_Iz)^{K_I}\right]\\ &=\exp\!\left(M|I|\big((1-p_I+p_Iz)-1\big)\right)\\ &=\exp\!\left(\left(\int_a^b\lambda(s)\,ds\right)(z-1)\right). \end{aligned}

The second equality uses the probability generating function \(\mathbb{E}[w^{K_I}]=\exp(M|I|(w-1))\) of a Poisson random variable, with \(w=1-p_I+p_Iz\). The final expression is the generating function of a Poisson random variable with mean \(\int_a^b\lambda(s)\,ds\). Hence the accepted count on \(I\) has the required Poisson law.

For disjoint intervals, the candidate processes restricted to those intervals are independent, by independent increments of the homogeneous Poisson process. The marks are also independent across candidates. Therefore the accepted counts on disjoint intervals are independent. There are almost surely finitely many candidates on \([0,T]\), so the construction gives a well-defined counting process with \(N_0=0\). These facts verify the defining properties of a nonhomogeneous Poisson process on the stated time interval. \(\square\)

The proof also explains the mechanism: on a short region where the intensity is a fraction \(r\) of the proposal rate, approximately that fraction of candidates is retained. The argument establishes the exact count law over every interval, including when \(\lambda\) varies continuously rather than being piecewise constant.

Worked Examples

Worked Example: One Simulated Candidate Path

On \([0,2]\), let \(\lambda(t)=1\) for \(0\le t<1\) and \(\lambda(t)=2\) for \(1\le t\le2\). Use \(M=2\). Suppose the generated candidate times are \(0.2,0.7,1.1,1.6\), with respective independent marks \(0.4,0.7,0.9,0.3\). Determine the retained events.

The proposed rate is valid: in the first part of the interval, \(\lambda(t)/M=1/2\), and in the second part, \(\lambda(t)/M=2/2=1\). At time \(0.2\), \(0.4\le1/2\), so retain the event. At time \(0.7\), \(0.7>1/2\), so discard it. At time \(1.1\), \(0.9\le1\), so retain it. At time \(1.6\), \(0.3\le1\), so retain it. The retained event times are therefore \(0.2,1.1,1.6\).

Check: Each candidate in \([0,1)\) was compared with \(1/2\), and the candidate times at least \(1\) were compared with \(1\). The path has three retained events; that is a single random outcome, not a claim that the expected count is three. The mean count on \([0,2]\), computed from the intensity, is \(\int_0^1 1\,dt+\int_1^2 2\,dt=1+2=3\).

Worked Example: Counts on Two Intervals

Let \(\lambda(t)=2\) on \([0,\tfrac12]\) and \(\lambda(t)=1\) on \((\tfrac12,\tfrac32]\), with \(M=2\). Find the distributions of the accepted counts on each interval, and the probability of no accepted events over the whole time range.

On the first interval the length is \(\tfrac12\), so the candidate mean is \(2\cdot\tfrac12=1\). The acceptance probability at each candidate time there is \(2/2=1\). Thus every candidate is retained, and the accepted count is Poisson with mean \(1\).

On the second interval the length is \(1\), so the candidate mean is \(2\cdot1=2\). The acceptance probability is \(1/2\). The exactness result gives an accepted count that is Poisson with mean \(1\cdot1=1\). The counts in the two intervals are independent. Consequently, the probability of no accepted events in either interval is the product of the two zero-count probabilities, \(e^{-1}e^{-1}=e^{-2}\).

Check: Direct integration gives total intensity \(2\cdot\tfrac12+1\cdot1=1+1=2\), so the full-interval count is Poisson with mean \(2\), whose probability of zero is \(e^{-2}\), as computed using the independent interval counts.

Worked Example: A Linear Intensity and Proposal Efficiency

On \([0,2]\), take \(\lambda(t)=1+t\) and \(M=3\). Check the bound, calculate the expected numbers of proposals and accepted events, and find the acceptance probability for a candidate at \(t=\tfrac32\).

For \(0\le t\le2\), \(1\le1+t\le3\), so \(M=3\) is a valid bound. The expected number of proposals is \(MT=3\cdot2=6\). The target integrated intensity is

\int_0^2(1+t)\,dt =\left[t+\frac{t^2}{2}\right]_0^2 =(2+2)-(0+0)=4.

Therefore the expected number of accepted events is \(4\), and the expected number discarded is \(6-4=2\). At \(t=\tfrac32\), the intensity is \(1+\tfrac32=\tfrac52\), so the acceptance probability is \((\tfrac52)/3=\tfrac56\).

Check: The expected proposal count splits into \(4\) accepted and \(2\) rejected, summing to \(6\). The integrated target mean is \(4\), agreeing with the expected accepted count required by the nonhomogeneous Poisson law.

Proposal Cost and Choosing the Bound

Proposition (Expected proposals and rejections): Under the assumptions of the exactness theorem, the expected number of candidates on \([0,T]\) is \(MT\), the expected number retained is \(\int_0^T\lambda(s)\,ds\), and the expected number rejected is \(MT-\int_0^T\lambda(s)\,ds\).

Proof. The candidate count has a Poisson distribution with mean \(MT\), so its expectation is \(MT\). Conditional on \(K=k\) candidates on \([0,T]\), each candidate has a uniform location and an independent mark. The probability that a candidate is accepted, averaged over its location and mark, is \(\frac{1}{MT}\int_0^T\lambda(s)\,ds\). Hence the conditional expected accepted count is \(k\frac{1}{MT}\int_0^T\lambda(s)\,ds\). Taking expectations and using \(\mathbb{E}K=MT\) gives \(\int_0^T\lambda(s)\,ds\). The number rejected equals the number proposed minus the number accepted, so its expectation is \(MT-\int_0^T\lambda(s)\,ds\). \(\square\)

A bound close to the largest intensity usually reduces wasted proposals. For a nonconstant intensity, choosing \(M\) far above its maximum still gives the correct accepted process, but increases the expected number of candidates. If the intensity is known only through a valid upper bound, correctness takes priority over efficiency: a loose valid bound is safe, while an invalid one is not.

Common Mistakes and a Practical Checklist

Common mistake: using an envelope that is too small. For \(\lambda(t)=1+t\) on \([0,2]\), taking \(M=1\) would give a supposed acceptance probability \(\lambda(\tfrac32)/M=(\tfrac52)/1=\tfrac52\) at \(t=\tfrac32\). This is not a probability. Even if one replaced values above \(1\) by \(1\), the resulting accepted rate at that time could be no greater than the proposal rate \(1\), whereas the target rate is \(\tfrac52\). The correction is to use a valid bound, such as \(M=3\), for which the acceptance probability at \(t=\tfrac32\) is \(\tfrac56\).

Common mistake: accepting with probability \(M/\lambda(t)\). That reverses the thinning ratio and can exceed one when the target rate is below the proposal rate. The correct rule is \(\lambda(t)/M\): a candidate is kept with probability equal to target intensity divided by proposal intensity.

Before running the algorithm, check the following:

Specify the finite time horizon and ensure the intensity is defined throughout it.
Choose \(M\) with \(0\le\lambda(t)\le M<\infty\) throughout the interval.
Generate candidate events at rate \(M\), not at the target rate.
Use a fresh independent Uniform\([0,1]\) mark for each candidate and retain it when \(U\le\lambda(t)/M\).
Interpret the output as the events on the specified time interval; the proof guarantees the Poisson count laws and independent increments there.

Check Your Understanding

For each question, state any bound you use and verify it against the intensity.

On \([0,1]\), let \(\lambda(t)=2-t\). Is \(M=2\) a valid bound? What is the acceptance probability at \(t=\tfrac12\)?
If a candidate process has rate \(M=4\) and the target intensity is \(1\) throughout an interval of length \(3\), what are the expected candidate and accepted counts on that interval?
For \(\lambda(t)=t\) on \([0,3]\), choose a valid positive constant bound and calculate the probability of accepting a candidate at \(t=1\).
Why are accepted counts on disjoint intervals independent in the thinning construction?
What goes wrong if a proposed bound \(M\) is smaller than the target intensity at some time?