Simpson 2016 - LGCP Likelihood and the Lattice Approximation

Summary

An inhomogeneous Poisson process has likelihood ; an LGCP puts , a Gaussian random field. The classical fit bins points into lattice cells with independent Poisson counts, which is a valid GLMM but wasteful because the grid must be fine enough to resolve point locations, not just the field.

Overview

Point data arise in ecology, biology and economics; MCMC for Cox processes is hard to tune, so INLA on a lattice (Illian et al. 2012) became popular. Simpson et al. argue the lattice plays two roles and the second is unnatural.

Main Content

Inhomogeneous Poisson process likelihood (eq. 1)

On bounded with intensity surface , , , and for a pattern

Log-Gaussian Cox process (LGCP)

A Cox (doubly stochastic Poisson) process with , a Gaussian random field; conditional on it is an inhomogeneous Poisson process. It models aggregation caused by observed or unobserved environmental variation and is a latent Gaussian model, so INLA applies. Log-likelihood given : .

Lattice (counts-in-cells) approximation

Partition into cells ; counts independent given , with . This is a GLMM and converges as cell size (Waagepetersen 2004; Corollary A1).

Why wasteful: (a) a dense covariance for limits lattice size; block-Toeplitz tricks (Moller et al. 1998) break under second-order approximations (manifold MCMC, INLA); (b) CAR models are fast but need regular lattices and are hard to make resolution-consistent on irregular ones; (c) binning is the dominant error (Corollary A1 gives only ), so the grid must be much finer than the field needs; (d) no local refinement: unsampled regions (Sec. 7.3) cost as much as sampled ones.

Examples

Gridded LGCP likelihood as used in PyMC practice

The gridded likelihood is exactly what E07/E86-style PyMC models use with a Hilbert-space GP: counts per cell .

Untested PyMC 5 sketch (2-D HSGP via pm.gp.HSGP):

import numpy as np, pymc as pm
# y: counts per grid cell (n_cells,), X: cell centres (n_cells, 2), A: cell areas
with pm.Model():
    mu = pm.Normal("mu", 0, 2)
    ell = pm.Gamma("ell", 3, 3 / 0.2); eta = pm.Exponential("eta", 1)
    cov = eta**2 * pm.gp.cov.Matern32(2, ls=ell)
    gp = pm.gp.HSGP(m=[20, 20], c=1.5, cov_func=cov)
    f = gp.prior("f", X=X)
    pm.Poisson("y", mu=A * pm.math.exp(mu + f), observed=y)   # counts-in-cells

Pitfall (E07): exp-link makes tails heavy; use tight priors on (“taming priors behind an exp link”). Read Hilbert Space Gaussian Processes for the basis details.

Connections

Simpson 2016 - Off-Grid Poisson-Form Likelihood removes the binning; Simpson 2016 - Convergence Theory quantifies the error.

See Also

Simpson 2016 - Going Off Grid - Overview