nstat.extras.spatial.lgcp

Log-Gaussian Cox process rate maps by the Laplace approximation.

Pure NumPy/SciPy. Turns a spatial point pattern into a posterior rate map with calibrated credible bands — the uncertainty a bare kernel-density rate map cannot provide.

The model:

1. Bin the events onto a grid; cell counts \(y_m\) are exact integers (Møller, Syversveen & Waagepetersen 1998).

2. Place a Gaussian-process prior on the log-intensity field \(f = \log\Lambda\), \(f \sim \mathcal N(m_0, K)\) with a Matern-5/2 covariance.

3. Find the posterior mode by Newton / IRLS (Rasmussen & Williams, Algorithm 3.1) with the Poisson observation model \(y_m \sim \text{Poisson}(a\,e^{f_m})\) (a = cell area):

   \[W = \operatorname{diag}(a\,e^{f}), \qquad f \leftarrow m_0 + (K^{-1} + W)^{-1} \bigl[W(f - m_0) + (y - a\,e^{f})\bigr].\]

4. The posterior covariance at the mode is \(\Sigma = (K^{-1} + \hat W)^{-1}\), with per-cell log-rate variance \(v_m = \Sigma_{mm}\).

The rate map is then read off as a log-normal band (the calibrated payoff): for a credible level with two-sided z-score \(z\),

\[\text{rate\_mean} = e^{\hat f + v/2}, \qquad \text{rate\_lo/hi} = e^{\hat f \mp z\sqrt{v}} .\]

The band is wider in data-sparse cells: where there are no events, \(\hat W \to 0\), so \(v_m \to K_{mm}\) (the GP prior variance) — this is the property the companion test asserts.

An optional heavier GP path may bridge gpflow (install pip install nstat-toolbox[spatial-gp]); the default is dependency-free pure NumPy/SciPy.

References

• Rasmussen CE, Williams CKI (2006). Gaussian Processes for Machine Learning, Algorithm 3.1 (Laplace approximation for GP classification), adapted here to the Poisson likelihood.

• Diggle PJ et al. (2013). Spatial and spatio-temporal log-Gaussian Cox processes. Statistical Science 28(4):542.

Functions

lgcp_fit(points, domain, *[, grid, kernel, ...])	Fit an LGCP rate map by the Laplace approximation.
lgcp_fit_glm(points, domain, basis, prior, *)	Basis-projected LGCP fit by penalized Poisson IRLS.

Classes

LGCPResult(grid, counts, f_mode, f_var, ...)	Fitted LGCP rate map (plain NumPy).
MaternPrior(nu, length_scale[, ...])	Matern GP prior on the coefficient vector of a 2-D B-spline basis.

class nstat.extras.spatial.lgcp.LGCPResult(grid, counts, f_mode, f_var, cell_area, edges, n_iter, converged)

Bases: object

Fitted LGCP rate map (plain NumPy).

Variables:

• grid (numpy.ndarray) – (M, d) cell-centre coordinates.

• counts (numpy.ndarray) – (M,) exact integer cell counts.

• f_mode (numpy.ndarray) – (M,) posterior-mode log-intensity \(\hat f\).

• f_var (numpy.ndarray) – (M,) posterior log-intensity variance \(v_m\).

• cell_area (float) – Area (or volume) of one grid cell.

• edges (list) – Per-dimension bin edges (length G+1 each).

• n_iter (int) – Newton/IRLS iterations to convergence.

• converged (bool) – Whether the mode-finding hit the tolerance before max_iter.

Parameters:

• grid (ndarray)

• counts (ndarray)

• f_mode (ndarray)

• f_var (ndarray)

• cell_area (float)

• edges (list)

• n_iter (int)

• converged (bool)

grid: ndarray

counts: ndarray

f_mode: ndarray

f_var: ndarray

cell_area: float

edges: list

n_iter: int

converged: bool

rate_map(level=0.9)

Posterior rate map with a log-normal credible band.

Parameters:

level (float) – Two-sided credible level (0.90 → \(z = 1.645\)).

Returns:

mean, lo, hi – Each (M,): the log-normal posterior mean rate \(e^{\hat f + v/2}\) and the lower/upper band \(e^{\hat f \mp z\sqrt v}\). Rates are per unit area (multiply by cell_area for expected counts).

Return type:

tuple[np.ndarray, np.ndarray, np.ndarray]

rate_mean()

Convenience: just the log-normal posterior-mean rate.

Return type:

ndarray

intensity_fn()

Return a callable X -> rate (nearest-cell lookup of the mean).

Useful as the lambda_hat argument to the nstat.extras.spatial.spatial_gof estimators — though note the plug-in caveat (use a held-out intensity for honest coverage).

nstat.extras.spatial.lgcp.lgcp_fit(points, domain, *, grid=20, kernel='matern52', length_scale=0.12, variance=1.0, prior_mean=None, max_iter=50, tol=1e-08, jitter=1e-06, backend='numpy')

Fit an LGCP rate map by the Laplace approximation.

Parameters:

• points (ndarray) – (n, d) event coordinates.

• domain – ((lo, hi), ...) rectangular analysis window per dimension.

• grid (int) – Number of cells per dimension (G). The analysis grid is G**d cells.

• kernel (str) – GP covariance. "matern52" (default, the standard smooth-but-not-too-smooth LGCP prior), "matern32", or "exponential" (Matern-1/2).

• length_scale (float) – GP range parameter \(\ell\).

• variance (float) – GP marginal variance \(\sigma^2\).

• prior_mean (float | None) – Constant prior mean \(m_0\) for the log-rate. None (default) uses \(\log(n/|D|) - \sigma^2/2\) so the prior-mean rate matches the empirical event density.

• max_iter (int) – Newton/IRLS stopping controls.

• tol (float) – Newton/IRLS stopping controls.

• jitter (float) – Diagonal jitter on the covariance for numerical stability.

• backend (str) – "numpy" (default, dependency-free) or "gpflow" (requires pip install nstat-toolbox[spatial-gp]). The gpflow path is a thin convenience that fits a variational GP and reads back the same posterior-mode/variance interface; it is optional and the NumPy path is the reference implementation.

Returns:

Call LGCPResult.rate_map() for the credible-band rate map.

Return type:

LGCPResult

Notes

Confidence: high — this is the standard Laplace GP (Rasmussen & Williams Alg. 3.1) with the Poisson link; the recovered mean tracks the true field and the band widens in data-sparse cells (both asserted in the companion tests).

nstat.extras.spatial.lgcp.lgcp_fit_glm(points, domain, basis, prior, *, grid=32, prior_mean=None, max_iter=50, tol=1e-06)

Basis-projected LGCP fit by penalized Poisson IRLS.

Bins events into a grid x grid cell grid, projects the log-intensity onto a 2-D B-spline basis (BSplineBasis2D), and finds the posterior mode of the basis coefficients under a MaternPrior evaluated at the Greville abscissae of the basis. The resulting per-cell log-rate posterior \((\hat f, v)\) is returned as the same LGCPResult that lgcp_fit() produces — so LGCPResult.rate_map() is available unchanged.

The cubic cost of the dominant linear solve scales with the basis dimension \(K\), not the grid size \(M\), so this is the routine to reach for on grids where the dense-GP lgcp_fit() would be slow (e.g. \(G \ge 50\)).

Parameters:

• points (ndarray) – (n, 2) event coordinates.

• domain – ((xlo, xhi), (ylo, yhi)) rectangular window.

• basis – A BSplineBasis2D whose grid_x / grid_y are the cell-centre axes of the analysis grid (length grid each). Use BSplineBasis2D.from_grid() to construct.

• prior (MaternPrior) – MaternPrior evaluated at the basis’ Greville abscissae (de Boor 1978) — the prior on the coefficient vector.

• grid (int) – Cells per axis (default 32). Must equal len(basis.grid_x) / len(basis.grid_y).

• prior_mean (float | ndarray | None) – Constant prior mean \(m_0\) for the log-rate (added on the coefficient axis; scalar broadcasts). None (default) sets \(m_0 = \log(\max(n,1)/|D|) - \sigma^2/2\) (matching lgcp_fit()’s default and reflecting the log-normal mean correction of Møller et al. 1998).

• max_iter (int) – Newton/IRLS stopping controls.

• tol (float) – Newton/IRLS stopping controls.

Returns:

With the same field semantics as lgcp_fit(): f_mode and f_var are per-cell log-rates (so rate_map() works without modification).

Return type:

LGCPResult

Notes

Confidence: high — this is the standard penalized GLM/GAM Laplace approximation (Wood 2017, Ch. 3; Diggle et al. 2013) with a Matern prior on B-spline coefficients. The basis-projection device is Diggle, Moraga, Rowlingson & Taylor 2013, Statistical Science.

References

• Møller J, Syversveen AR, Waagepetersen RP (1998). Log Gaussian Cox processes. Scand. J. Statistics 25(3):451-482.

• Rasmussen CE, Williams CKI (2006). Gaussian Processes for Machine Learning.

• Diggle PJ, Moraga P, Rowlingson B, Taylor BM (2013). Spatial and spatio-temporal log-Gaussian Cox processes. Statistical Science 28(4):542.

• Wood SN (2017). Generalized Additive Models: An Introduction with R.

See also

lgcp_fit

dense per-cell GP fit; the reference implementation (cheaper at small G, expensive at large G).

class nstat.extras.spatial.lgcp.MaternPrior(nu, length_scale, marginal_var=1.0, jitter=1e-06)

Bases: object

Matern GP prior on the coefficient vector of a 2-D B-spline basis.

Used by lgcp_fit_glm() to place a smooth GP prior directly on the \(K\) B-spline coefficients (evaluated at their Greville abscissae, de Boor 1978) rather than on the \(M = G^2\) cell log-rates. When \(K \ll M\) this is the difference between an \(O(K^3)\) IRLS step and an \(O(M^3)\) one — for a 64x64 grid (\(M=4096\)) with a typical \(K=64\)-coefficient basis the speedup is ~64x in the dominant cubic cost.

Parameters:

• nu (float) – Matern smoothness; one of {0.5, 1.5, 2.5} (closed-form Matern).

• length_scale (float) – Range parameter \(\ell > 0\).

• marginal_var (float) – Marginal variance \(\sigma^2 > 0\) (default 1.0).

• jitter (float) – Diagonal stabilizer added to \(K\) (default 1e-6); on a numpy.linalg.LinAlgError from Cholesky the jitter is retried once at 10 * jitter.

Notes

Derived quantities (each a method taking coords): K is the dense (K, K) covariance; cholesky the lower Cholesky factor \(L\) with \(K = LL^\top\); K_inv the inverse \(K^{-1}\) via scipy.linalg.cho_solve(); log_det \(\log\det K = 2\sum_i \log L_{ii}\).

Each derived quantity is computed once per coords and cached on the instance, keyed by the array’s shape and bytes.

References

• Rasmussen CE, Williams CKI (2006). Gaussian Processes for Machine Learning, Ch. 4.

• Møller J, Syversveen AR, Waagepetersen RP (1998). Log Gaussian Cox processes. Scand. J. Statistics 25(3):451-482.

nu: float

length_scale: float

marginal_var: float = 1.0

jitter: float = 1e-06

K(coords)

Dense Matern covariance (K, K) on coords.

Parameters:

coords (ndarray)

Return type:

ndarray

cholesky(coords)

Lower Cholesky factor of K().

Parameters:

coords (ndarray)

Return type:

ndarray

K_inv(coords)

Inverse covariance (computed via the cached factor).

Parameters:

coords (ndarray)

Return type:

ndarray

log_det(coords)

Log-determinant of K().

Parameters:

coords (ndarray)

Return type:

float