Simpson 2016 - Off-Grid Poisson-Form Likelihood

Summary

With a finite basis $Z(s)={\sum }_j{z}_j{\varphi }_j(s)$, the LGCP log-likelihood becomes $C-{\tilde{\alpha }}^{\top }\exp (A_{1}z)+1^{\top }{A}_{2}z$: the sum over points is exact and only the integral needs a $p$-node quadrature. This is the log-likelihood of $N+p$ independent Poisson variables with pseudo-observations, so it fits standard latent-Gaussian machinery (INLA).

Overview

The key idea is to separate (i) the Gaussian field model and (ii) the likelihood approximation.

Main Content

Finite-dimensional field (eq. 2)

$Z(s)={\sum }_{i=1}^n{z}_i{\varphi }_i(s)$, $z$ multivariate Gaussian, $\{{\varphi }_i\}$ linearly independent deterministic basis functions. Special cases: process convolution (Higdon 1998), predictive processes (Banerjee et al. 2008) and fixed-rank kriging (Cressie & Johannesson 2008) with correlated weights, and SPDE models (Lindgren et al. 2011) with compactly supported piecewise-linear ${\varphi }_i$.

Poisson-form likelihood (eq. 3-4)

For a quadrature rule ${\int }_{\Omega }f\approx {\sum }_{i=1}^p{\tilde{\alpha }}_{i}f({\tilde{s}}_i)$,

$$\log \pi (y\mid z)\approx C-{\tilde{\alpha }}^{\top }\exp (A_{1}z)+1^{\top }{A}_{2}z,[A_1{]}_{ij}={\varphi }_j({\tilde{s}}_i),[A_2{]}_{ij}={\varphi }_j(s_i).$$

Define $\log \eta =(z^{\top }{A}_1^{\top },z^{\top }{A}_2^{\top }{)}^{\top }$, $\alpha =({\tilde{\alpha }}^{\top },0_N^{\top }{)}^{\top }$ and pseudo-observations $y=(0_p^{\top },1_N^{\top }{)}^{\top }$. Then $\pi (y\mid z)\approx C{\prod }_{i=1}^{N+p}{\eta }_i^{y_i}{e}^{-{\alpha }_i{\eta }_i}$: $N+p$ conditionally independent Poisson variables with means ${\alpha }_i{\eta }_i$.

Why. The term ${\sum }_{i}Z(s_i)=1^{\top }{A}_{2}z$ is exact because $Z$ is continuous; the integral $\int e^Z$ is replaced by a weighted sum of $e^{Z({\tilde{s}}_i)}$. Weight-zero pseudo-points give ${\eta }_i^1$ factors (the $\prod \lambda (s_i)$ term). Related to Baddeley & Turner (2000), who used one point per cell for pseudo-likelihoods; here the integration scheme is designed carefully (Thm 1: error governed by quadrature order).

Integration scheme on a triangular mesh: dual-cell midpoint rule (nodes, weights ${\tilde{\alpha }}_i=|V_i|$ the dual cell volumes; second order on a regular grid, first order on irregular meshes), or Gaussian rules per triangle.

Examples

Numerical check of eq. 3 on a 1-D LGCP

Setup: domain $[0,1]$, piecewise-linear hat basis on $n$ nodes, quadrature = trapezoid with $p=200$ nodes. The code shows that the exact data term plus quadrature integral reproduces the Poisson-form likelihood.

import numpy as np
rng = np.random.default_rng(0)
n = 15; nodes = np.linspace(0, 1, n); z = rng.normal(0, 1, n).cumsum() * 0.3
phi = lambda s: np.maximum(0, 1 - np.abs(s[:, None] - nodes[None, :]) * (n - 1))   # hat functions
Z = lambda s: phi(s) @ z
pts = np.sort(rng.uniform(0, 1, 12))                 # observed points
sq = np.linspace(0, 1, 200); w = np.full(200, 1 / 199); w[[0, -1]] /= 2   # trapezoid
A1, A2 = phi(sq), phi(pts)
ll_direct = -(w * np.exp(Z(sq))).sum() + Z(pts).sum()
eta = np.concatenate([np.exp(A1 @ z), np.exp(A2 @ z)]); alpha = np.concatenate([w, np.zeros(len(pts))])
y = np.concatenate([np.zeros(len(sq)), np.ones(len(pts))])
ll_pois = (y * np.log(eta) - alpha * eta).sum()
print(ll_direct, ll_pois)   # identical

In PyMC: pass pseudo-data to pm.Poisson(mu=alpha*eta, observed=y), equivalent to pm.Potential with the explicit form.

Connections

• Simpson 2016 - Matern SPDE and GMRF Representation supplies ${\varphi }_i$, $Q$.
• Simpson 2016 - Convergence Theory gives Hellinger rates.
• The same trick (point evaluation matrix + quadrature matrix) is how HSGP LGCPs avoid gridding (Hilbert Space Gaussian Processes).

See Also

Simpson 2016 - Going Off Grid - Overview