Likelihood Inference for Temporal Poisson Intensities

From Event Times to a Likelihood

A temporal Poisson model describes both how many events occur and when they occur. If the rate changes over time, the event times themselves contain information about the shape of that change. This tutorial derives the likelihood for those times and uses it to estimate two common intensity models. The observation window is fixed at \([0,T]\), and events after \(T\) are not observed.

Let \((\Omega,\mathcal{F},\mathbb{P})\) be a probability space, and let \((N_t)_{t\ge0}\) be a nonhomogeneous Poisson process with state space \(S=\{0,1,2,\ldots\}\) and continuous time \(t\ge0\). Its deterministic intensity is \(\lambda(t;\theta)\ge0\), where \(\theta\) is an unknown parameter. Assume the intensity is integrable on \([0,T]\), and define the cumulative intensity \(\Lambda_\theta(T)=\int_0^T\lambda(u;\theta)\,du\). The count \(N_T\) is Poisson with mean \(\Lambda_\theta(T)\). When \(\Lambda_\theta(T)>0\), conditional on \(N_T=n\), the \(n\) event times are distributed as the order statistics of \(n\) independent draws with density \(\lambda(t;\theta)/\Lambda_\theta(T)\).

The likelihood is a density for the observed event times, not a probability of observing those exact real values. To derive it when \(n\ge1\) and \(\Lambda_\theta(T)>0\), multiply the probability of the count by the conditional density of the ordered times. The count probability is \(e^{-\Lambda_\theta(T)}\Lambda_\theta(T)^n/n!\). The joint density of the ordered times conditional on \(N_T=n\) is \(n!\prod_{k=1}^n[\lambda(t_k;\theta)/\Lambda_\theta(T)]\) on the ordered region. Multiplication cancels both \(n!\) and \(\Lambda_\theta(T)^n\), giving the stated likelihood. For \(n=0\), the count probability itself is \(e^{-\Lambda_\theta(T)}\), as required. If \(\Lambda_\theta(T)=0\), the process has no events in the window with probability one; an observed positive count then has likelihood zero.

Taking logarithms gives the log-likelihood \(\ell(\theta)=\sum_{k=1}^n\log\lambda(t_k;\theta)-\int_0^T\lambda(u;\theta)\,du\), whenever the intensity is positive at the observed event times. The first term rewards intensity where events occurred. The integral is an exposure penalty: a model assigning a large rate throughout the observation window predicts many events, including at times when none were observed. For a differentiable model with an interior parameter value, its score vector is \(U(\theta)=\nabla\ell(\theta)=\sum_k\nabla\log\lambda(t_k;\theta)-\int_0^T\nabla\lambda(u;\theta)\,du\). An interior maximum must have score zero. That equation is necessary, but by itself it does not establish that a solution is a maximum; curvature and boundary behavior must also be checked.

Constant Rate: The Exposure Penalty in Action

Exponential-Trend Intensity

A simple model for a rate that rises or falls smoothly is \(\lambda(t;a,b)=a e^{bt}\), with \(a>0\) and \(b\in\mathbb{R}\). Here \(a\) is the rate at time zero and \(b\) controls the temporal trend: positive \(b\) means an increasing rate, while negative \(b\) means a decreasing rate. Write \(A(b)=\int_0^T e^{bu}\,du\). The log-likelihood is \(\ell(a,b)=n\log a+b\sum_k t_k-aA(b)\). The score equations are \(n/a-A(b)=0\) and \(\sum_k t_k-aA'(b)=0\). The first equation gives \(\widehat a(b)=n/A(b)\) when \(n>0\). Substitution into the second yields a one-dimensional equation for \(b\).

Proof. From \(\widehat a(b)=n/A(b)\), the profiled log-likelihood is \(\ell_{\mathrm{prof}}(b)=n\log n-n\log A(b)+b\sum_k t_k-n\). Differentiating gives \(\ell_{\mathrm{prof}}'(b)=\sum_k t_k-nA'(b)/A(b)\). Define \(f_b(u)=e^{bu}/A(b)\) on \([0,T]\); its integral is \(1\), so it is a probability density. Differentiating \(A\) under the integral sign gives \(A'(b)=\int_0^T ue^{bu}\,du\) and \(A''(b)=\int_0^T u^2e^{bu}\,du\). Therefore \(m(b)=\int_0^T uf_b(u)\,du\), and \(m'(b)=A''(b)/A(b)-(A'(b)/A(b))^2=\int_0^T u^2f_b(u)\,du-(\int_0^T uf_b(u)\,du)^2\), which is the variance under \(f_b\). Since \(T>0\) and \(f_b\) is positive throughout the interval, this distribution is not concentrated at one value; its variance is strictly positive. Consequently \(\ell_{\mathrm{prof}}''(b)=-nm'(b)<0\), so the profile score is strictly decreasing and has at most one root. \(\square\)

Piecewise-Constant Intensities

Suppose the window is divided into disjoint intervals \(I_j\) of lengths \(\Delta_j>0\), and the rate is a separate constant \(\lambda_j\ge0\) on each interval. Let \(n_j\) be the number of observed events in \(I_j\). Because the integrated intensity is \(\sum_j\lambda_j\Delta_j\) and each event contributes the rate of its interval, the likelihood factors as \(\prod_j\lambda_j^{n_j}e^{-\lambda_j\Delta_j}\). For \(n_j>0\), the log-likelihood contribution \(n_j\log\lambda_j-\lambda_j\Delta_j\) has score \(n_j/\lambda_j-\Delta_j\), so its maximum is \(\widehat\lambda_j=n_j/\Delta_j\). Its second derivative is \(-n_j/\lambda_j^2<0\), and the contribution tends to \(-\infty\) at both ends of the positive parameter range. If \(n_j=0\), the likelihood contribution is \(e^{-\lambda_j\Delta_j}\), maximized at \(\lambda_j=0\) when zero is allowed. If the model insists on \(\lambda_j>0\), there is no attained maximum in that interval; the likelihood supremum is approached as \(\lambda_j\) decreases to zero.

Uncertainty, Interpretation, and a Common Error

The observed information is \(J(\widehat\theta)=-\nabla^2\ell(\widehat\theta)\). Under the usual regularity conditions for likelihood asymptotics—an identifiable model, a true interior parameter, sufficient smoothness, positive nonsingular information, and an appropriate large-sample regime—\(J(\widehat\theta)^{-1}\) is used to approximate the covariance matrix of the estimator. This is the standard likelihood-based uncertainty approximation; it is not a universal exact formula. In a temporal model, a large-sample regime might involve repeated independent observation windows or increasing observation time, but the relevant conditions still have to hold. Boundary estimates, such as a zero rate in an interval with no events, do not satisfy the interior-parameter condition, so ordinary symmetric standard-error intervals may be inappropriate.

For a piecewise-constant model, the counts in disjoint intervals are independent Poisson variables with means \(\lambda_j\Delta_j\). Thus \(n_j/\Delta_j\) is unbiased for \(\lambda_j\), and its exact variance is \(\lambda_j/\Delta_j\): a Poisson variable with mean \(\lambda_j\Delta_j\) has variance \(\lambda_j\Delta_j\), and dividing by \(\Delta_j^2\) gives the stated variance. Substituting the fitted rate yields the familiar information-based estimate \(n_j/\Delta_j^2\) when \(n_j>0\). This exact variance calculation for the piecewise model does not make the exponential-trend standard errors exact.

Common pitfall: omitting the integrated intensity. With three events in a window of length \(2\) under a constant rate, using only the event factors gives the incorrect likelihood \(L_{\mathrm{wrong}}(\theta)=\theta^3\). It grows without bound as \(\theta\) increases, so it has no finite maximum. The correct likelihood includes the exposure term: \(L(\theta)=\theta^3e^{-2\theta}\). Its score \(3/\theta-2\) vanishes at \(3/2\), and its negative second derivative there is \(4/3>0\), confirming a finite maximum. The missing exponential term is exactly what penalizes rates that predict excessive events over the observation window.

The key distinction is between the count-only likelihood and the event-time likelihood. In a constant-rate model, event locations do not change the estimate once the count and window length are known. In a time-varying model, the event locations enter through \(\prod_k\lambda(t_k;\theta)\), while the integrated intensity still accounts for the full exposure. Keeping both pieces is essential for estimating temporal patterns and interpreting uncertainty.