\[ \newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right)} \newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]} \newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]} \newcommand{\SampleVar}[1]{\widehat{\mathrm{Var}}\left[ #1 \right]} \newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]} \newcommand{\TrueRegFunc}{\mu} \newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}} \DeclareMathOperator{\tr}{tr} \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator*{\argmax}{argmax} \DeclareMathOperator{\det}{det} \DeclareMathOperator{\dim}{dim} \newcommand{\TrueNoise}{\epsilon} \newcommand{\EstNoise}{\widehat{\TrueNoise}} \newcommand{\Signal}{S} \newcommand{\SignalNoise}{N} \newcommand{\AutoCov}{\gamma} \newcommand{\NoiseAutoCov}{\xi} \newcommand{\LinearExpect}[1]{\mathbb{L}\left[ #1 \right]} \newcommand{\Periodicity}{\tau} \newcommand{\Loadings}{\mathbf{b}} \newcommand{\FactorNoise}{\epsilon} \newcommand{\Uniqueness}{\mathbf{\psi}} \newcommand{\AltLoadings}{\mathbf{c}} \newcommand{\AltFactorNoise}{\zeta} \newcommand{\AltUniqueness}{\mathbf{\rho}} \newcommand{\DynMatrix}{\mathbf{a}} \newcommand{\DynNoise}{\eta} \newcommand{\DynVar}{\xi} \newcommand{\KalmanGain}{\mathbf{K}} \newcommand{\indep}{\perp} \newcommand{\Neighborhood}{\partial} \]

Previously

Data consist of scalars / vectors / categoricals at particular coordinates, \(X(r,t)\)
We have been interested in relating \(X(r,t)\) to \(X(q,s)\) (or to \(Y(r,t)\) or to \(Y(q,s)\), etc.)
Today: the data is just the coordinates: “something happened there” or “something happened then”

Examples of Point Data

… with some selected citations

“Spikes” from neurons in neuroscience (Rieke et al. 1997; Brown, Kass, and Mitra 2004)
Locations of plants in forestry / ecology (O’Sullivan and Perry 2013, sec. 2.2)
- Ignoring distinctions among plants of the same species
Locations of stars or galaxies in astronomy (Neyman and Scott 1952)
- Ignoring distinctions among types of stars
Locations of disease cases (Diggle 2014, ch. 9)
- Or locations and times of infections (Diggle 2014, ch. 10)
Locations and times of crimes (Reinhart and Greenhouse 2018)

Point process

A random process whose value is a set of points, also called events
- Contrast to an ordinary stochastic process, whose value is a function of time and/or space
Equivalently: a random measure giving the number of events in every set
- This measure is exclusively atomic \(=\) a sum of delta functions
- This is not a probability measure (unless there’s only one event)

Point process as a limit of binary-valued stochastic processes

Take time first for simplicity, say the interval \([0,T]\)
Divide time up into bins of width \(h\)
\(X(i) = 1\) if there was an event in bin \(i\)
- \(i \in 1: (T/h)\); don’t worry if \(T/h\) is not an integer
Now imagine shrinking the bins, \(h \rightarrow 0\)

Point process vs. counting process vs. time points

The counting process \(N(t) \equiv\) number of events in the interval \([0, t)\) \[\begin{eqnarray} N(t) & = & \sum_{i=0}^{t/h}{X(i)}\\ X(i) & = & N(ih) - N((i-1)h) \end{eqnarray}\]
The time coordinates \(\tau_1, \tau_2, \ldots \tau_{N(T)}\) tell us when events happened and let us recover \(X\) or \(N\), and vice versa
- Inter-event intervals: \(\tau_j - \tau_{j-1}\)

The most basic model

Events happen with a constant probability, independently across bins: \[ X(i) \sim_{IID} \mathrm{Bern}(p) \]
Counting process is binomial: \[ N(t) \sim \mathrm{Binom}(t/h, p) \]
Counting increments are independent and binomial: \[ N(t+s) - N(t) = \sum_{(t/h)+1}^{(t+s)/h}{X(i)} = \mathrm{Binom}(s/h, p) \]
Inter-event times are geometric: \[ \tau_j - \tau_{j-1} \sim \mathrm{Geom}(p) \]
- To wait \(k\) time-steps for the next event, we need \(k-1\) failures followed by 1 success

Shrink the bins

Let \(h \rightarrow 0\) but force \(p \propto h\), or at least \(p/h \rightarrow \lambda\)
Limits follow from basic prob: \[\begin{eqnarray} N(t) & \sim & \mathrm{Pois}(\lambda t)\\ N(t+s) - N(t) & \sim & \mathrm{Pois}(\lambda s)\\ \tau_{i+1} - \tau_{i} & \sim & \mathrm{Exp}(\lambda) \end{eqnarray}\]
Counting increments and inter-event intervals are still independent
\(X\) is now a Poisson process
\(\lambda\) is the intensity of \(X\)

Estimating the homogeneous Poisson process

Use the inter-event-time representation, \(\tau_1, \tau_2, \ldots \tau_{N_T}\)
- For convenience, set \(\tau_0 = 0\) and \(\tau_{N_T+1} = T\)
The log likelihood is \[\begin{eqnarray} \mathcal{L}(\lambda) & = & \log{\lambda^{N_T}} + \sum_{i=1}^{N_T+1}{\log{e^{-\lambda(\tau_i - \tau_{i-1})}}}\\ & = & N_T\log{\lambda} - \lambda \sum_{i=1}^{N_T+1}{\left(\tau_i - \tau_{i-1}\right)}\\ \label{eqn:homogenous-poisson-likelihood} & = & N_T\log{\lambda} - \lambda T \end{eqnarray}\]
Maximizing this gives \[ \hat{\lambda} = N_T/T \]
Since \(N_T \sim \mathrm{Pois}(\lambda T)\), this is unbiased and consistent
- \(\Expect{\hat{\lambda}} = T^{-1} (\lambda T) = \lambda\), \(\Var{\hat{\lambda}} = T^{-2} (\lambda T) = \lambda/T\)
EXERCISE: Can you get this log-likelihood as the limit of the log-likelihood of the Bernoulli process? Can you get this MLE as the limit of the Bernoulli MLEs?

Inhomogenous Poisson Processes

Now let \(\lambda(t)\) be a function of time, the intensity function
Mythology: we start at zero and put down points independently, with the probability per unit time, at time \(t\), being \(\lambda_t\)
We can still equally describe the process by the points or the counting process or the inter-event intervals
The counting process is still Poisson-distributed, with independent increments \[\begin{eqnarray} N(t) & \sim & \mathrm{Pois}\left(\int_{0}^{t}{\lambda(r) dr}\right) \\ N(t+s) - N(t) & \sim & \mathrm{Pois}\left(\int_{t}^{t+s}{\lambda(r) dr}\right) \end{eqnarray}\]

Time-rescaling inhomogeneous Poisson processes

Define re-scaled time \[ \tilde{\tau}(t) = \int_{0}^{t}{\lambda(r) dr} \] and re-scaled event times \(\tilde{\tau}_i\) for each of the event times \(\tau_i\)
- Stretch out physical time where the intensity is high, compress physical time where the intensity is low
- EXERCISE (easy): Re-express the counting process in terms of the re-scaled time
Re-scaled inter-event intervals are independently and exponentially distributed: \[ \tilde{\tau}_i - \tilde{\tau}_{i-1} \sim \mathrm{Exp}(1) \]
This gives us an easy way to test goodness-of-fit to a proposed \(\lambda(t)\): use it to re-scale the observed event times, and then test for fit to an exponential distribution (Brown et al. 2002)

More General Point Processes

Inhomogeneous Poisson process: \(\lambda(t)\) is a deterministic function of time
In general, the intensity can be a function of the past of the process and/or covariates
These will be random, so the intensity will be random so write \(\Lambda(t)\)
Typically \(\Lambda(t) = \lambda(t|\mathcal{H}(t))\), for some variable \(\mathcal{H}(t)\) recording the relevant history up to and including \(t\)
- Really \(\mathcal{H}(t)\) is a filtration
Intensities are still related to counts: \[ \Expect{N(t) - N(s)|\mathcal{H}(s)} = \int_{s}^{t}{\Expect{\lambda(r|\mathcal{H}(r))|\mathcal{H}(s)} dr} \]
- Seems to require averaging over all paths of the point process from time \(s\) to time \(t\)
  - Can sometimes short-cut this
Conversely: \[ \frac{d\Expect{N(t)}}{dt} = \Lambda_t \]
- Intensity \(=\) expected rate of increase of the count
Time-rescaling still works, and still gives a nice goodness-of-fit test (Brown et al. 2002)

Self-Exciting Processes

A self-exciting or Hawkes process (Hawkes 1971): \[ \lambda(t|\mathcal{H}_t) = \lambda_0 + \sum_{i: \tau_i < t}{m(t-\tau_i)} \] for some memory kernel or memory functuon \(m\)
Can be tricky to allow for inhibition as well as excitation
- Could always make the intensity multiplicative instead

The Likelihood for Point Processes

With event times \(\tau_1, \tau_2, \ldots \tau_{N_T}\) and intensity function \(\lambda(t)\), the log-likelihood becomes \[ \mathcal{L}(\lambda) = \sum_{i=1}^{N_T}{\log{\lambda(\tau_i)}} -\int_{0}^{T}{\lambda(t) dt} \]
- Go back to the discrete-time-bins picture and shrink the bins
- The sum favors putting high intensity around the observed data
- The integral penalizes just increasing the intensity everywhere
The fully nonparametric MLE is \[ \sum_{i=1}^{N_T}{\delta(t-\tau_i)} \] but this is not (usually) useful

Parametric Likelihood Estimates

Suppose \[ \lambda(t) = \lambda_0 + \beta z(t) \] for given covariate function \(z(t)\)
Maximizing the log-likelihood gives (EXERCISE) \[\begin{eqnarray} \sum_{i=1}^{n_T}{\frac{1}{\hat{\lambda}_0 + \hat{\beta}z(\tau_i)}} & = & T\\ \sum_{i=1}^{n_T}{\frac{z(\tau_i)}{\hat{\lambda}_0 + \hat{\beta}z(\tau_i)}} & = & T \int_{0}^{T}{z(t) dt} \end{eqnarray}\]
For a self-exciting Hawkes process, the estimating equations are not quite so nice, but one can still optimize

Non-Parametric Estimation (I)

Replace points with smooth smears

Start with the fully non-parametric MLE \[ \sum_{i=1}^{N_T}{\delta(t-\tau_i)} \]
Now replace \(\delta\) functions with, e.g., Gaussians or other smooth pdfs
- Basically the same idea as kernel density estimation
- Which pdf?
- How much bandwidth? (Use cross-validation)

Non-Parametric Estimation (II)

Smooth the counting process

\[ \frac{d\Expect{N_t}}{dt} = \Lambda_t \]

Get the realized counting process \(N\) over \([0,T]\)
Use nonparametric regression to get \(\widehat{\Expect{N(t)}}\)
- Use a method like splines that’s consistent for the first derivative
Differentiate to get the intensity: \[ \widehat{\Lambda(t)} = \frac{d}{dt}\widehat{\Expect{N(t)}} \]

Spatial and Spatio-Temporal Processes

The counting process becomes a measure, \(N(A) =\) number of events in set \(A\)
- Knowing \(N(A)\) is still equivalent to knowing all the event locations
There’s still an intensity
The expected count of events of a region is the integral of the intensity \[ \Expect{N(A)} = \int_{A}{\lambda(r) dr} \]
\(\lambda(r)\) is the expected growth rate (per unit area or volume) in the number of events as we expand a ball around \(r\) \[ \lim_{h\downarrow 0}{\frac{\Expect{N(B(r,h))}}{|B(r,h)|}} = \lambda(r) \] where \(B(r,h) =\) ball of radius \(h\) centered at \(r\)
- Similarly for \(\lambda(r,t)\)
Knowing the distances between all the events is no longer equivalent to knowing the location of all the events

Spatial and Spatio-Temporal Inference

Log-likelihood still takes the form \[ \mathcal{L}(\lambda) = \sum_{i}{\log{\lambda(r_i)}} -\int_{\mathbb{D}}{\lambda(r) dr} \] where \(\mathbb{D}\) is the region of observation
- Again, go back to the discrete-bins picture and take the limit
BUT lots of good models involve lots of latent variables making likelihood computationally painful

Example: Neyman-Scott Model (Neyman and Scott 1952)

Create centers \(Z_1, \ldots Z_N\) according to a homogeneous Poisson process
For each center \(i \in 1:N\)
- Draw \(K_i \sim g(k)\), a pmf on integers
- Draw points \(X_1, \ldots X_{K_i} \sim f(x-Z_i)\), \(f\) a pdf on the space
  - e.g., uniform on a ball of some radius
Observe the \(X\)s, not the \(Z\)s
Model is easy to simulate, but not easy to deal with analytically
- EXERCISE: write the simulator
Calculating the intensity without the \(Z\)s is hard!
- Still hard if you observe the \(Z\)s but can’t tell them from the \(X\)s

Summing Up

Point processes are, equivalently, random sets of isolated points, or random atomic measures
The most basic situation is the homogeneous Poisson process, and other processes get built up from there
The equivalent of the trend is the intensity function, which is \(\approx\) the density of the expected random measure
Likelihood-based inference is possible but quickly gets difficult as models grow even a little bit complicated
We will look at simulation-based inference for point processes in 741

Backup: Topics deliberately omitted, to save room

In no particular order:

The full taxonomy of point processes with stochastic intensities (Cox processes vs. doubly stochastic Poisson processes vs. …)
- Relatedly, the full (semi-) martingale theory of intensities as stochastic processes
Marked point processes (where there is a finite set of types of events, i.e., each point gets one of a finite number of “marks”)
Second-order properties, essentially “how do you define correlation within a point process” (e.g., Ripley’s \(K\) function) or across processes
Intensity functions conditional on the presence or absence of events, conditional distributions more generally (“Palm measure”)
Fourier transforms for point processes
Details on relating point processes to the history of (vector-valued) covariates; “decoding” the covariates from the point process (Rieke et al. 1997)

References

Brown, Emery N., Riccardo Barbieri, Valérie Ventura, Robert E. Kass, and Loren M. Frank. 2002. “The Time-Rescaling Theorem and Its Applications to Neural Spike Train Data Analysis.” Neural Computation 14. https://doi.org/https://doi.org/10.1162/08997660252741149.

Brown, Emery N., Robert E. Kass, and Partha P. Mitra. 2004. “Multiple Neural Spike Train Data Analysis: State-of-the-Art and Future Challanges.” Nature Neuroscience 7:456–61. http://www.stat.cmu.edu/~kass/papers/bkm.pdf.

Diggle, Peter J. 2014. Statistical Analysis of Spatial and Spatio-Temporal Point Patterns. Third. Boca Raton, Florida: CRC Press.

Hawkes, Alan G. 1971. “Spectra of Some Self-Exciting and Mutually Exciting Point Processes.” Biometrika 58:83–90. https://doi.org/10.2307/2334319.

Neyman, J., and E. L. Scott. 1952. “A Theory of the Spatial Distribution of Galaxies.” Astrophysical Journal 116:144–63. https://doi.org/10.1086/145599.

O’Sullivan, David, and George L. W. Perry. 2013. Spatial Simulation: Exploring Pattern and Process. Chichester, England: Wiley.

Reinhart, Alex, and Joel Greenhouse. 2018. “Self-Exciting Point Processes with Spatial Covariates: Modeling the Dynamics of Crime.” Journal of the Royal Statistical Society C 67:1305–29. https://doi.org/10.1111/rssc.12277.

Rieke, Fred, David Warland, Rob de Ruyter van Steveninck, and William Bialek. 1997. Spikes: Exploring the Neural Code. Cambridge, Massachusetts: MIT Press.