nstat.extras.spatial.spatial_gof

Inhomogeneous second-order spatial goodness-of-fit (pure NumPy/SciPy).

A Python-only inhomogeneous second-order goodness-of-fit suite for spatial point processes. It provides the intensity-reweighted second-order summary statistics that test what a fitted intensity \(\hat\lambda(\mathbf{x})\) leaves over — the diagnostics a homogeneous \(K\)-function cannot give for a non-stationary neural field:

  • pair_correlation() — the SOIRS-reweighted pair correlation \(g(r)\); g(r) > 1 indicates clustering, g(r) < 1 repulsion, g(r) = 1 the (inhomogeneous) Poisson null.

  • k_inhom() — the inhomogeneous \(K\)-function of Baddeley, Moller & Waagepetersen (2000); for an inhomogeneous Poisson process K_inhom(r) = pi r^2 in 2-D.

  • l_function() — the variance-stabilized \(L(r) = \sqrt{K(r)/\pi}\), so L(r) - r = 0 under the null.

  • nearest_neighbour_FGJ() — the empty-space \(F\), nearest- neighbour \(G\), and \(J = (1-G)/(1-F)\) functions.

  • global_envelope() — a Monte-Carlo global-rank envelope test (Myllymaki et al. 2017) built by thinning the fitted inhomogeneous Poisson intensity.

Warning

Plug-in bias. Reweighting by an intensity \(\hat\lambda\) estimated from the same pattern deflates the variance of the estimator and shrinks the envelope below nominal coverage (the fitted field has already absorbed some of the clustering). Pass a held-out intensity for lambda_hat (e.g. a smoother fit to a disjoint fold, or a global-reweighting estimator; Shaw, Moller & Waagepetersen 2021) or treat the resulting coverage as optimistic. Every estimator below carries this caveat.

This is the static-spatial analogue of the time-rescaling KS test in nstat.extras.spatial.marked_gof; the SOIRS precondition (second- order intensity-reweighted stationarity — \(\lambda\) bounded away from zero, \(g\) a function of lag only) is the spatial mirror of the time-rescaling null.

References

  • Baddeley AJ, Moller J, Waagepetersen R (2000). Non- and semi-parametric estimation of interaction in inhomogeneous point patterns. Statistica Neerlandica 54(3):329.

  • Myllymaki M, Mrkvicka T, Grabarnik P, Seijo H, Hahn U (2017). Global envelope tests for spatial processes. JRSS-B 79(2):381.

  • Shaw T, Moller J, Waagepetersen R (2021). Globally intensity-reweighted estimators for K- and pair correlation functions. ANZJS 63(1):93.

Functions

cross_k_inhom(points_A, points_B, lambda_A, ...)

Inhomogeneous cross \(K\)-function (Baddeley-Moller-Waagepetersen 2000).

cross_pair_correlation(points_A, points_B, ...)

Cross pair correlation \(g_{AB}(r)\) for two label classes.

global_envelope(points, lambda_hat, r_grid, *)

Monte-Carlo global-rank envelope test against the inhomogeneous null.

k_inhom(points, lambda_hat, r_grid, *[, ...])

Inhomogeneous \(K\)-function (Baddeley-Moller-Waagepetersen 2000).

l_function(points, lambda_hat, r_grid, *[, ...])

Variance-stabilized inhomogeneous \(L\)-function.

nearest_neighbour_FGJ(points, r_grid, *[, ...])

Empty-space \(F\), nearest-neighbour \(G\), and \(J\).

pair_correlation(points, lambda_hat, r_grid, *)

SOIRS-reweighted inhomogeneous pair correlation \(g(r)\).
nstat.extras.spatial.spatial_gof.pair_correlation(points, lambda_hat, r_grid, *, bw=None, domain=None, edge_correction='epanechnikov')[source]

SOIRS-reweighted inhomogeneous pair correlation \(g(r)\).

Kernel estimator with an Epanechnikov edge kernel and the SOIRS reweighting \(1/(\hat\lambda_i\,\hat\lambda_j)\):

\[\hat g(r) = \frac{1}{2\pi r\,|D|} \sum_{i\neq j} \frac{k_{bw}(r - \|x_i - x_j\|)} {\hat\lambda(x_i)\,\hat\lambda(x_j)} .\]

Under the inhomogeneous Poisson null, \(g(r) = 1\) at all lags; g(r) > 1 is clustering and g(r) < 1 repulsion.

Parameters:

  • points (ndarray) – (n, d) event coordinates (d = 2 for the planar estimator).

  • lambda_hat – The reweighting intensity — a callable X -> lambda(X) or a per-point array. Use a held-out estimate (see module warning).

  • r_grid (ndarray) – Lags at which to evaluate \(g\).

  • bw (float | None) – Kernel bandwidth. Defaults to a Stoyan-style rule of thumb based on the lag spacing.

  • domain – ((lo, hi), ...) per dim. Defaults to the bounding box of points — pass the true analysis window when known.

  • edge_correction (str) – One of "epanechnikov" (default — the original SOIRS estimator), "isotropic" (Ripley 1976/1977, symmetric i<->j average), "translation" (Ohser 1983), or "border" (Baddeley-Rubak- Turner 2015). The three named modes require a rectangular domain; "border" returns NaN at radii where no event has boundary-distance >= r.

Returns:

\(\hat g(r)\) at each lag in r_grid.

Return type:

np.ndarray

Notes

Confidence: high on the homogeneous limit (g → 1 for Poisson) and the clustering/repulsion sign; the absolute scale is sensitive to the bandwidth and the plug-in caveat in the module docstring.

nstat.extras.spatial.spatial_gof.k_inhom(points, lambda_hat, r_grid, *, domain=None, edge_correction='epanechnikov')[source]

Inhomogeneous \(K\)-function (Baddeley-Moller-Waagepetersen 2000).

\[\hat K_{\text{inhom}}(r) = \frac{1}{|D|} \sum_{i\neq j} \frac{\mathbf{1}[\|x_i - x_j\| \le r]} {\hat\lambda(x_i)\,\hat\lambda(x_j)} ,\]

the standard estimator that reduces to the ordinary \(K\) when \(\hat\lambda\) is constant. For an inhomogeneous Poisson process, K_inhom(r) = pi r^2 in 2-D (the CSR benchmark).

Parameters:
Returns:

\(\hat K_{\text{inhom}}(r)\) at each lag.

Return type:

np.ndarray

nstat.extras.spatial.spatial_gof.l_function(points, lambda_hat, r_grid, *, domain=None, edge_correction='epanechnikov')[source]

Variance-stabilized inhomogeneous \(L\)-function.

\(L(r) = \sqrt{K_{\text{inhom}}(r)/\pi}\) (2-D), so the centred statistic \(L(r) - r\) is flat at zero under the inhomogeneous Poisson null. Returns L(r) (subtract r_grid for the centred form). The edge_correction kwarg passes through to k_inhom(); "border" may yield NaN rows where no event is usable as a focal point.

Parameters:
Return type:

ndarray

nstat.extras.spatial.spatial_gof.nearest_neighbour_FGJ(points, r_grid, *, domain=None, n_dummy=None, rng=None)[source]

Empty-space \(F\), nearest-neighbour \(G\), and \(J\).

  • \(G(r)\) — CDF of the nearest-neighbour distance from each event to the nearest other event.

  • \(F(r)\) — empty-space function: CDF of the distance from a grid of reference (dummy) locations to the nearest event.

  • \(J(r) = (1 - G(r))/(1 - F(r))\) — equals 1 under CSR; J < 1 clustering, J > 1 regularity.

These are homogeneous (unreweighted) summaries; they characterize inter-event spacing and are most useful as a complement to the reweighted \(g\)/\(K\) on roughly-stationary patterns.

Returns:

F, G, J – Each aligned with r_grid.

Return type:

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

Parameters:
nstat.extras.spatial.spatial_gof.global_envelope(points, lambda_hat, r_grid, *, n_sim=199, domain=None, statistic='pcf', bw=None, alpha=0.05, lambda_grid=None, grid_points=None, edge_correction='epanechnikov', rng=None)[source]

Monte-Carlo global-rank envelope test against the inhomogeneous null.

Simulates n_sim realizations of the fitted inhomogeneous Poisson process (by thinning a dominating homogeneous process at the maximum intensity over a reference grid), computes the chosen summary statistic on each, and builds the global-rank envelope of Myllymaki et al. (2017). The observed pattern’s statistic is tested for global containment.

Parameters:

  • points (ndarray) – See pair_correlation(). lambda_hat (held-out) is used both to reweight the statistic and as the simulation intensity.

  • lambda_hat – See pair_correlation(). lambda_hat (held-out) is used both to reweight the statistic and as the simulation intensity.

  • r_grid (ndarray) – See pair_correlation(). lambda_hat (held-out) is used both to reweight the statistic and as the simulation intensity.

  • domain – See pair_correlation(). lambda_hat (held-out) is used both to reweight the statistic and as the simulation intensity.

  • bw (float | None) – See pair_correlation(). lambda_hat (held-out) is used both to reweight the statistic and as the simulation intensity.

  • n_sim (int) – Number of Monte-Carlo null simulations (199 → 95% global level with alpha = 0.05).

  • statistic (str) – "pcf" for \(g(r)\) (default), "kinhom" for \(K_{\text{inhom}}(r)\), or "linhom" for \(L(r)\).

  • alpha (float) – Global type-I error level.

  • lambda_grid – Optional precomputed intensity values and their coordinates, used to find the dominating intensity for thinning. If omitted, a default reference grid over domain is built and lambda_hat evaluated on it.

  • grid_points (ndarray | None) – Optional precomputed intensity values and their coordinates, used to find the dominating intensity for thinning. If omitted, a default reference grid over domain is built and lambda_hat evaluated on it.

  • edge_correction (str) – Edge correction forwarded to the per-curve summary statistic ("epanechnikov" default, "isotropic", "translation", or "border"). The default branch is bit-identical to the pre-kwarg behaviour.

  • rng (Generator | None) – NumPy random generator for the simulations.

Returns:

Carries observed, lo, hi, inside, p_interval.

Return type:

EnvelopeResult

Notes

Honours the plug-in caveat: if lambda_hat was fit to the same pattern, the envelope coverage is optimistic. Prefer a held-out intensity. Confidence: high on the mechanics; coverage validity is conditional on the held-out reweighting.

nstat.extras.spatial.spatial_gof.cross_k_inhom(points_A, points_B, lambda_A, lambda_B, r_grid, *, domain=None, edge_correction='epanechnikov')[source]

Inhomogeneous cross \(K\)-function (Baddeley-Moller-Waagepetersen 2000).

For two disjoint label classes \(A\) and \(B\) with intensities \(\hat\lambda_A\) and \(\hat\lambda_B\),

\[\hat K_{\text{inhom}}^{AB}(r) = \frac{1}{|D|} \sum_{i \in A}\sum_{j \in B} \frac{\mathbf{1}[\|x_i - x_j\| \le r]} {\hat\lambda_A(x_i)\,\hat\lambda_B(x_j)} .\]

Under independent inhomogeneous Poisson labels (the bivariate Poisson null), \(K_{\text{inhom}}^{AB}(r) = \pi r^2\) in 2-D — the CSR-pair benchmark.

The edge-correction modes mirror the single-type k_inhom(): "epanechnikov" (default, unweighted), and the rectangular-domain modes "isotropic" (symmetric Ripley average over the two events of each pair), "translation" (Ohser 1983), and "border" (focal-event boundary distance >= r).

Parameters:

  • points_A (ndarray) – Event coordinates of each label class. No self-exclusion is applied: the disjoint label hypothesis is assumed.

  • points_B (ndarray) – Event coordinates of each label class. No self-exclusion is applied: the disjoint label hypothesis is assumed.

  • lambda_A – Reweighting intensities — callables X -> lambda(X) or per-point arrays aligned with points_A / points_B.

  • lambda_B – Reweighting intensities — callables X -> lambda(X) or per-point arrays aligned with points_A / points_B.

  • r_grid (ndarray) – Upper radii at which to evaluate the cumulative statistic.

  • domain – ((lo, hi), ...) per dim. Defaults to the bounding box of np.vstack([points_A, points_B]).

  • edge_correction (str) – See k_inhom().

Returns:

\(\hat K_{\text{inhom}}^{AB}(r)\) at each lag.

Return type:

np.ndarray

Notes

Confidence: high on the independent-pair limit (\(K^{AB} \to \pi r^2\) for independent inhomogeneous Poisson label classes); same plug-in caveat as the single-type estimator.

nstat.extras.spatial.spatial_gof.cross_pair_correlation(points_A, points_B, lambda_A, lambda_B, r_grid, *, bw=None, domain=None, edge_correction='epanechnikov')[source]

Cross pair correlation \(g_{AB}(r)\) for two label classes.

Kernel estimator with the Epanechnikov edge kernel and the cross-SOIRS reweighting \(1 / (\hat\lambda_A(x_i)\, \hat\lambda_B(x_j))\):

\[\hat g_{AB}(r) = \frac{1}{2\pi r\,|D|} \sum_{i \in A}\sum_{j \in B} \frac{k_{bw}(r - \|x_i - x_j\|)} {\hat\lambda_A(x_i)\,\hat\lambda_B(x_j)} .\]

Under independent inhomogeneous Poisson labels, \(g_{AB}(r) = 1\) at every lag; \(> 1\) indicates cross-attraction, < 1 cross-repulsion.

Parameters:
Returns:

\(\hat g_{AB}(r)\) at each lag in r_grid.

Return type:

np.ndarray

Notes

Confidence: high on the bivariate-Poisson limit (\(g_{AB} \to 1\) for independent label classes). Same plug-in caveat as the single-type estimator — prefer held-out \(\hat\lambda_A, \hat\lambda_B\).

class nstat.extras.spatial.spatial_gof.EnvelopeResult(r_grid, observed, lo, hi, inside, p_interval, n_sim)[source]

Bases: object

Result of a global-rank envelope test.

Variables:

  • r_grid (numpy.ndarray) – The lag values at which the summary statistic was evaluated.

  • observed (numpy.ndarray) – The observed summary curve, shape (len(r_grid),).

  • lo (numpy.ndarray) – Lower / upper global-rank envelope at the requested level, each shape (len(r_grid),).

  • hi (numpy.ndarray) – See lo.

  • inside (bool) – True iff the observed curve lies inside [lo, hi] at every lag — i.e. the global test does NOT reject the null.

  • p_interval (tuple[float, float]) – Conservative / liberal global-rank p-value interval (p_lo, p_hi) (Myllymaki et al. 2017, eq. for the rank envelope test).

  • n_sim (int) – Number of Monte-Carlo simulations used.

Parameters:
r_grid: ndarray
observed: ndarray
lo: ndarray
hi: ndarray
inside: bool
p_interval: tuple[float, float]
n_sim: int