r"""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
:math:`\hat\lambda(\mathbf{x})` leaves over — the diagnostics a
homogeneous :math:`K`-function cannot give for a non-stationary neural
field:
- :func:`pair_correlation` — the SOIRS-reweighted pair correlation
:math:`g(r)`; ``g(r) > 1`` indicates clustering, ``g(r) < 1`` repulsion,
``g(r) = 1`` the (inhomogeneous) Poisson null.
- :func:`k_inhom` — the inhomogeneous :math:`K`-function of
Baddeley, Moller & Waagepetersen (2000); for an inhomogeneous Poisson
process ``K_inhom(r) = pi r^2`` in 2-D.
- :func:`l_function` — the variance-stabilized
:math:`L(r) = \sqrt{K(r)/\pi}`, so ``L(r) - r = 0`` under the null.
- :func:`nearest_neighbour_FGJ` — the empty-space :math:`F`, nearest-
neighbour :math:`G`, and :math:`J = (1-G)/(1-F)` functions.
- :func:`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 :math:`\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
:mod:`nstat.extras.spatial.marked_gof`; the SOIRS precondition (second-
order intensity-reweighted stationarity — :math:`\lambda` bounded away
from zero, :math:`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.
"""
from __future__ import annotations
import numpy as np
from scipy.spatial.distance import cdist, pdist
from nstat.extras.spatial._edge import (
border_usable_mask,
frac_disc_in_rect,
frac_translation_rect,
)
from nstat.extras.spatial._envelopes import EnvelopeResult, global_rank_envelope
from nstat.extras.spatial._kernels import epanechnikov
_VALID_EDGE_CORRECTIONS = ("epanechnikov", "isotropic", "translation", "border")
def _validate_edge_correction(name: str) -> None:
if name not in _VALID_EDGE_CORRECTIONS:
raise ValueError(
"edge_correction must be one of "
"{'epanechnikov','isotropic','translation','border'}; "
f"got {name!r}"
)
def _ensure_rectangular_domain(domain, mode: str) -> None:
"""Check that ``domain`` is a 2-tuple of 2-tuples for a rectangle."""
ok = (
isinstance(domain, tuple)
and len(domain) == 2
and all(isinstance(d, tuple) and len(d) == 2 for d in domain)
)
if not ok:
raise NotImplementedError(
"only rectangular domain is supported for non-epanechnikov edge "
"corrections in v0.6.0"
)
# ----------------------------------------------------------------------
# Internal helpers
# ----------------------------------------------------------------------
def _as_points(points: np.ndarray) -> np.ndarray:
points = np.atleast_2d(np.asarray(points, dtype=float))
if points.ndim != 2:
raise ValueError(f"points must be (n, d); got shape {points.shape}")
return points
def _intensity_at(lambda_hat, points: np.ndarray) -> np.ndarray:
"""Evaluate the reweighting intensity at the given points.
``lambda_hat`` may be (a) a callable ``X -> lambda(X)``, or (b) a
1-D array aligned with ``points`` giving lambda at each point.
"""
if callable(lambda_hat):
lam = np.asarray(lambda_hat(points), dtype=float).ravel()
else:
lam = np.asarray(lambda_hat, dtype=float).ravel()
if lam.shape[0] != points.shape[0]:
raise ValueError(
"lambda_hat array must align with points; got "
f"{lam.shape[0]} vs {points.shape[0]}"
)
if np.any(lam <= 0):
raise ValueError(
"lambda_hat must be strictly positive (SOIRS requires lambda "
"bounded away from zero); got a non-positive value."
)
return lam
def _domain_measure(domain) -> float:
return float(np.prod([hi - lo for lo, hi in domain]))
# ----------------------------------------------------------------------
# Pair correlation g(r)
# ----------------------------------------------------------------------
[docs]
def pair_correlation(
points: np.ndarray,
lambda_hat,
r_grid: np.ndarray,
*,
bw: float | None = None,
domain=None,
edge_correction: str = "epanechnikov",
) -> np.ndarray:
r"""SOIRS-reweighted inhomogeneous pair correlation :math:`g(r)`.
Kernel estimator with an Epanechnikov edge kernel and the SOIRS
reweighting :math:`1/(\hat\lambda_i\,\hat\lambda_j)`:
.. math::
\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, :math:`g(r) = 1` at all lags;
``g(r) > 1`` is clustering and ``g(r) < 1`` repulsion.
Parameters
----------
points
``(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
Lags at which to evaluate :math:`g`.
bw
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
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
-------
np.ndarray
:math:`\hat g(r)` at each lag in ``r_grid``.
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.
"""
_validate_edge_correction(edge_correction)
pts = _as_points(points)
r_grid = np.asarray(r_grid, dtype=float)
if pts.shape[0] < 2:
return np.zeros_like(r_grid)
if pts.shape[1] != 2:
raise ValueError("pair_correlation is implemented for d=2 (planar).")
if domain is None:
domain = tuple((pts[:, k].min(), pts[:, k].max()) for k in range(pts.shape[1]))
area = _domain_measure(domain)
if bw is None:
bw = float(0.1 / np.sqrt(pts.shape[0] / area)) if area > 0 else 0.05
bw = max(bw, float(np.mean(np.diff(r_grid))) if len(r_grid) > 1 else 0.05)
lam = _intensity_at(lambda_hat, pts)
iu = np.triu_indices(pts.shape[0], k=1)
d = cdist(pts, pts)[iu]
wgt = 1.0 / (lam[iu[0]] * lam[iu[1]])
if edge_correction == "epanechnikov":
g = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
ring = 2.0 * np.pi * r * area
if ring <= 0:
g[k] = 0.0
continue
g[k] = 2.0 * np.sum(epanechnikov((d - r) / bw) / bw * wgt) / ring
return g
# Named-correction branches share the rectangular-window precondition.
_ensure_rectangular_domain(domain, edge_correction)
pi_i, pi_j = pts[iu[0]], pts[iu[1]]
if edge_correction == "isotropic":
g = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
ring = 2.0 * np.pi * r * area
if ring <= 0:
g[k] = 0.0
continue
ker = epanechnikov((d - r) / bw) / bw
# Per-pair Ripley weight, symmetric in i<->j at the pair distance.
w_iso = np.empty_like(d)
for m in range(d.size):
if ker[m] == 0.0:
w_iso[m] = 0.0
continue
fi = frac_disc_in_rect(pi_i[m], d[m], domain)
fj = frac_disc_in_rect(pi_j[m], d[m], domain)
w_iso[m] = 0.5 * (1.0 / fi + 1.0 / fj)
g[k] = 2.0 * np.sum(ker * wgt * w_iso) / ring
return g
if edge_correction == "translation":
# Translation weight depends only on the inter-event offset.
offs = pi_i - pi_j
w_tr = np.array(
[frac_translation_rect(o[0], o[1], domain) for o in offs], dtype=float
)
g = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
ring = 2.0 * np.pi * r * area
if ring <= 0:
g[k] = 0.0
continue
ker = epanechnikov((d - r) / bw) / bw
g[k] = 2.0 * np.sum(ker * wgt * w_tr) / ring
return g
# edge_correction == "border": restrict the focal points per radius.
g = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
ring = 2.0 * np.pi * r * area
if ring <= 0:
g[k] = 0.0
continue
usable = border_usable_mask(pts, domain, r)
if not usable.any():
g[k] = np.nan
continue
# Keep pairs whose first index is a usable focal point.
keep = usable[iu[0]]
if not keep.any():
g[k] = np.nan
continue
ker = epanechnikov((d[keep] - r) / bw) / bw
# Reduced "effective" window is the area of the eroded rectangle.
(xlo, xhi), (ylo, yhi) = domain
eff_area = max(0.0, (xhi - xlo) - 2 * r) * max(0.0, (yhi - ylo) - 2 * r)
if eff_area <= 0:
g[k] = np.nan
continue
ring_eff = 2.0 * np.pi * r * eff_area
g[k] = 2.0 * np.sum(ker * wgt[keep]) / ring_eff
return g
# ----------------------------------------------------------------------
# Inhomogeneous K and L
# ----------------------------------------------------------------------
[docs]
def k_inhom(
points: np.ndarray,
lambda_hat,
r_grid: np.ndarray,
*,
domain=None,
edge_correction: str = "epanechnikov",
) -> np.ndarray:
r"""Inhomogeneous :math:`K`-function (Baddeley-Moller-Waagepetersen 2000).
.. math::
\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 :math:`K` when
:math:`\hat\lambda` is constant. For an inhomogeneous Poisson
process, ``K_inhom(r) = pi r^2`` in 2-D (the CSR benchmark).
Parameters
----------
points, lambda_hat, domain
See :func:`pair_correlation`. **Held-out** ``lambda_hat``.
r_grid
Upper radii at which to evaluate the cumulative statistic.
edge_correction
See :func:`pair_correlation` — one of ``"epanechnikov"`` (default,
identical to the previous implementation), ``"isotropic"``,
``"translation"``, or ``"border"``.
Returns
-------
np.ndarray
:math:`\hat K_{\text{inhom}}(r)` at each lag.
"""
_validate_edge_correction(edge_correction)
pts = _as_points(points)
r_grid = np.asarray(r_grid, dtype=float)
if pts.shape[0] < 2:
return np.zeros_like(r_grid)
if domain is None:
domain = tuple((pts[:, k].min(), pts[:, k].max()) for k in range(pts.shape[1]))
area = _domain_measure(domain)
lam = _intensity_at(lambda_hat, pts)
iu = np.triu_indices(pts.shape[0], k=1)
d = pdist(pts)
wgt = 1.0 / (lam[iu[0]] * lam[iu[1]])
if edge_correction == "epanechnikov":
K = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
# factor 2: the (i, j) and (j, i) ordered pairs both contribute.
K[k] = 2.0 * np.sum(wgt[d <= r]) / area
return K
_ensure_rectangular_domain(domain, edge_correction)
pi_i, pi_j = pts[iu[0]], pts[iu[1]]
if edge_correction == "isotropic":
# Per-pair isotropic weight is evaluated at the inter-event distance.
fi = np.array([frac_disc_in_rect(pi_i[m], d[m], domain) for m in range(d.size)])
fj = np.array([frac_disc_in_rect(pi_j[m], d[m], domain) for m in range(d.size)])
w_iso = 0.5 * (1.0 / fi + 1.0 / fj)
K = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
mask = d <= r
K[k] = 2.0 * np.sum(wgt[mask] * w_iso[mask]) / area
return K
if edge_correction == "translation":
offs = pi_i - pi_j
w_tr = np.array(
[frac_translation_rect(o[0], o[1], domain) for o in offs], dtype=float
)
K = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
mask = d <= r
K[k] = 2.0 * np.sum(wgt[mask] * w_tr[mask]) / area
return K
# border
K = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
usable = border_usable_mask(pts, domain, r)
if not usable.any():
K[k] = np.nan
continue
keep = usable[iu[0]] & (d <= r)
if not keep.any():
K[k] = 0.0
continue
(xlo, xhi), (ylo, yhi) = domain
eff_area = max(0.0, (xhi - xlo) - 2 * r) * max(0.0, (yhi - ylo) - 2 * r)
if eff_area <= 0:
K[k] = np.nan
continue
K[k] = 2.0 * np.sum(wgt[keep]) / eff_area
return K
[docs]
def l_function(
points: np.ndarray,
lambda_hat,
r_grid: np.ndarray,
*,
domain=None,
edge_correction: str = "epanechnikov",
) -> np.ndarray:
r"""Variance-stabilized inhomogeneous :math:`L`-function.
:math:`L(r) = \sqrt{K_{\text{inhom}}(r)/\pi}` (2-D), so the centred
statistic :math:`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
:func:`k_inhom`; ``"border"`` may yield ``NaN`` rows where no event is
usable as a focal point.
"""
K = k_inhom(points, lambda_hat, r_grid, domain=domain, edge_correction=edge_correction)
# NaN propagates through sqrt and 0.0 floor.
out = np.where(np.isnan(K), np.nan, np.sqrt(np.maximum(K, 0.0) / np.pi))
return out
# ----------------------------------------------------------------------
# Nearest-neighbour F / G / J
# ----------------------------------------------------------------------
[docs]
def nearest_neighbour_FGJ(
points: np.ndarray,
r_grid: np.ndarray,
*,
domain=None,
n_dummy: int | None = None,
rng: np.random.Generator | None = None,
):
r"""Empty-space :math:`F`, nearest-neighbour :math:`G`, and :math:`J`.
- :math:`G(r)` — CDF of the nearest-neighbour distance from each event
to the nearest other event.
- :math:`F(r)` — empty-space function: CDF of the distance from a grid
of reference (dummy) locations to the nearest event.
- :math:`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 :math:`g`/:math:`K` on roughly-stationary patterns.
Returns
-------
F, G, J : tuple[np.ndarray, np.ndarray, np.ndarray]
Each aligned with ``r_grid``.
"""
pts = _as_points(points)
r_grid = np.asarray(r_grid, dtype=float)
rng = np.random.default_rng() if rng is None else rng
if domain is None:
domain = tuple((pts[:, k].min(), pts[:, k].max()) for k in range(pts.shape[1]))
n = pts.shape[0]
if n < 2:
z = np.zeros_like(r_grid)
return z, z, np.ones_like(r_grid)
# G: nearest-neighbour distances among events.
dmat = cdist(pts, pts)
np.fill_diagonal(dmat, np.inf)
nn = dmat.min(axis=1)
G = np.array([np.mean(nn <= r) for r in r_grid])
# F: empty-space distances from dummy points to nearest event.
if n_dummy is None:
n_dummy = max(100, 4 * n)
lows = np.array([lo for lo, _ in domain])
highs = np.array([hi for _, hi in domain])
dummies = lows + (highs - lows) * rng.uniform(size=(n_dummy, pts.shape[1]))
de = cdist(dummies, pts).min(axis=1)
F = np.array([np.mean(de <= r) for r in r_grid])
with np.errstate(divide="ignore", invalid="ignore"):
J = np.where((1.0 - F) > 1e-12, (1.0 - G) / (1.0 - F), np.nan)
return F, G, J
# ----------------------------------------------------------------------
# Global-rank envelope (Monte-Carlo)
# ----------------------------------------------------------------------
[docs]
def global_envelope(
points: np.ndarray,
lambda_hat,
r_grid: np.ndarray,
*,
n_sim: int = 199,
domain=None,
statistic: str = "pcf",
bw: float | None = None,
alpha: float = 0.05,
lambda_grid=None,
grid_points: np.ndarray | None = None,
edge_correction: str = "epanechnikov",
rng: np.random.Generator | None = None,
) -> EnvelopeResult:
r"""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, lambda_hat, r_grid, domain, bw
See :func:`pair_correlation`. ``lambda_hat`` (held-out) is used
both to reweight the statistic and as the simulation intensity.
n_sim
Number of Monte-Carlo null simulations (199 → 95% global level
with ``alpha = 0.05``).
statistic
``"pcf"`` for :math:`g(r)` (default), ``"kinhom"`` for
:math:`K_{\text{inhom}}(r)`, or ``"linhom"`` for :math:`L(r)`.
alpha
Global type-I error level.
lambda_grid, grid_points
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
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
NumPy random generator for the simulations.
Returns
-------
EnvelopeResult
Carries ``observed``, ``lo``, ``hi``, ``inside``, ``p_interval``.
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.*
"""
_validate_edge_correction(edge_correction)
pts = _as_points(points)
r_grid = np.asarray(r_grid, dtype=float)
rng = np.random.default_rng() if rng is None else rng
if domain is None:
domain = tuple((pts[:, k].min(), pts[:, k].max()) for k in range(pts.shape[1]))
area = _domain_measure(domain)
stat_fns = {
"pcf": lambda X: pair_correlation(
X, lambda_hat, r_grid, bw=bw, domain=domain,
edge_correction=edge_correction,
),
"kinhom": lambda X: k_inhom(
X, lambda_hat, r_grid, domain=domain,
edge_correction=edge_correction,
),
"linhom": lambda X: l_function(
X, lambda_hat, r_grid, domain=domain,
edge_correction=edge_correction,
),
}
if statistic not in stat_fns:
raise ValueError(f"statistic must be one of {list(stat_fns)}; got {statistic!r}")
stat_fn = stat_fns[statistic]
# Dominating intensity for the thinning simulation.
if grid_points is None or lambda_grid is None:
n_ref = 40
lows = np.array([lo for lo, _ in domain])
highs = np.array([hi for _, hi in domain])
axes = [np.linspace(lo, hi, n_ref) for lo, hi in domain]
mesh = np.meshgrid(*axes, indexing="xy")
grid_points = np.column_stack([m.ravel() for m in mesh])
lambda_grid = _intensity_at(lambda_hat, grid_points)
else:
grid_points = _as_points(grid_points)
lambda_grid = np.asarray(lambda_grid, dtype=float).ravel()
lam_max = float(np.max(lambda_grid))
observed = stat_fn(pts)
sims = np.empty((n_sim, len(r_grid)))
lows = np.array([lo for lo, _ in domain])
highs = np.array([hi for _, hi in domain])
for s in range(n_sim):
n_prop = rng.poisson(lam_max * area)
cand = lows + (highs - lows) * rng.uniform(size=(n_prop, pts.shape[1]))
lam_cand = _intensity_at(lambda_hat, cand) if n_prop > 0 else np.empty(0)
keep = rng.uniform(size=n_prop) < (lam_cand / lam_max) if n_prop > 0 else np.array([], bool)
Xs = cand[keep]
sims[s] = stat_fn(Xs)
return global_rank_envelope(observed, sims, r_grid, alpha=alpha)
# ----------------------------------------------------------------------
# Cross-type (bivariate) inhomogeneous second-order statistics
# ----------------------------------------------------------------------
def _cross_default_domain(pts_A: np.ndarray, pts_B: np.ndarray):
stacked = np.vstack([pts_A, pts_B])
return tuple((stacked[:, k].min(), stacked[:, k].max()) for k in range(stacked.shape[1]))
[docs]
def cross_k_inhom(
points_A: np.ndarray,
points_B: np.ndarray,
lambda_A,
lambda_B,
r_grid: np.ndarray,
*,
domain=None,
edge_correction: str = "epanechnikov",
) -> np.ndarray:
r"""Inhomogeneous cross :math:`K`-function (Baddeley-Moller-Waagepetersen 2000).
For two disjoint label classes :math:`A` and :math:`B` with
intensities :math:`\hat\lambda_A` and :math:`\hat\lambda_B`,
.. math::
\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), :math:`K_{\text{inhom}}^{AB}(r) = \pi r^2` in 2-D —
the *CSR-pair* benchmark.
The edge-correction modes mirror the single-type
:func:`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, points_B
Event coordinates of each label class. No self-exclusion is
applied: the disjoint label hypothesis is assumed.
lambda_A, lambda_B
Reweighting intensities — callables ``X -> lambda(X)`` or
per-point arrays aligned with ``points_A`` / ``points_B``.
r_grid
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
See :func:`k_inhom`.
Returns
-------
np.ndarray
:math:`\hat K_{\text{inhom}}^{AB}(r)` at each lag.
Notes
-----
*Confidence: high* on the independent-pair limit
(:math:`K^{AB} \to \pi r^2` for independent inhomogeneous Poisson
label classes); same plug-in caveat as the single-type estimator.
"""
_validate_edge_correction(edge_correction)
pts_A = _as_points(points_A)
pts_B = _as_points(points_B)
if pts_A.shape[1] != pts_B.shape[1]:
raise ValueError(
"points_A and points_B must share the spatial dimension; got "
f"{pts_A.shape[1]} vs {pts_B.shape[1]}"
)
r_grid = np.asarray(r_grid, dtype=float)
if pts_A.shape[0] == 0 or pts_B.shape[0] == 0:
return np.zeros_like(r_grid)
if domain is None:
domain = _cross_default_domain(pts_A, pts_B)
area = _domain_measure(domain)
lam_A = _intensity_at(lambda_A, pts_A)
lam_B = _intensity_at(lambda_B, pts_B)
d = cdist(pts_A, pts_B)
wgt = 1.0 / (lam_A[:, None] * lam_B[None, :])
if edge_correction == "epanechnikov":
K = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
mask = d <= r
K[k] = float(np.sum(wgt[mask])) / area
return K
_ensure_rectangular_domain(domain, edge_correction)
if edge_correction == "isotropic":
# Symmetric Ripley weight per pair at the pair distance, matching
# the single-type isotropic convention. The weight depends on the
# pair only, not on r, so compute it once for every pair some r in
# the grid reaches (the mask is cumulative, as in k_inhom) instead
# of once per bin. Each bin then sums the same products in the
# same row-major order as before, so K is bit-identical.
reach = np.zeros(d.shape, dtype=bool)
for r in r_grid:
reach |= d <= r
ii, jj = np.nonzero(reach)
w_iso = np.empty(ii.size, dtype=float)
for m in range(ii.size):
rij = d[ii[m], jj[m]]
fi = frac_disc_in_rect(pts_A[ii[m]], rij, domain)
fj = frac_disc_in_rect(pts_B[jj[m]], rij, domain)
w_iso[m] = 0.5 * (1.0 / fi + 1.0 / fj)
d_pair = d[ii, jj]
contrib = wgt[ii, jj] * w_iso
K = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
sel = d_pair <= r
if not sel.any():
K[k] = 0.0
continue
K[k] = float(np.sum(contrib[sel])) / area
return K
if edge_correction == "translation":
offs = pts_A[:, None, :] - pts_B[None, :, :]
# Vectorised translation weight: rectangle / overlap area.
(xlo, xhi), (ylo, yhi) = domain
Lx, Ly = float(xhi - xlo), float(yhi - ylo)
ax = np.abs(offs[..., 0])
ay = np.abs(offs[..., 1])
overlap = np.maximum(Lx - ax, 0.0) * np.maximum(Ly - ay, 0.0)
with np.errstate(divide="ignore", invalid="ignore"):
w_tr = np.where(overlap > 0, (Lx * Ly) / overlap, np.inf)
K = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
mask = d <= r
K[k] = float(np.sum(wgt[mask] * w_tr[mask])) / area
return K
# border: focal A-events with boundary distance >= r.
K = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
usable = border_usable_mask(pts_A, domain, r)
if not usable.any():
K[k] = np.nan
continue
(xlo, xhi), (ylo, yhi) = domain
eff_area = max(0.0, (xhi - xlo) - 2 * r) * max(0.0, (yhi - ylo) - 2 * r)
if eff_area <= 0:
K[k] = np.nan
continue
sub = wgt[usable, :]
sub_d = d[usable, :]
K[k] = float(np.sum(sub[sub_d <= r])) / eff_area
return K
[docs]
def cross_pair_correlation(
points_A: np.ndarray,
points_B: np.ndarray,
lambda_A,
lambda_B,
r_grid: np.ndarray,
*,
bw: float | None = None,
domain=None,
edge_correction: str = "epanechnikov",
) -> np.ndarray:
r"""Cross pair correlation :math:`g_{AB}(r)` for two label classes.
Kernel estimator with the Epanechnikov edge kernel and the
cross-SOIRS reweighting :math:`1 / (\hat\lambda_A(x_i)\,
\hat\lambda_B(x_j))`:
.. math::
\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, :math:`g_{AB}(r) =
1` at every lag; :math:`> 1` indicates *cross-attraction*, ``< 1``
*cross-repulsion*.
Parameters
----------
points_A, points_B, lambda_A, lambda_B, domain
See :func:`cross_k_inhom`.
r_grid
Lags at which to evaluate :math:`g_{AB}`.
bw
Kernel bandwidth; defaults to a Stoyan-style rule of thumb,
identical to the single-type :func:`pair_correlation` default.
edge_correction
See :func:`pair_correlation`.
Returns
-------
np.ndarray
:math:`\hat g_{AB}(r)` at each lag in ``r_grid``.
Notes
-----
*Confidence: high* on the bivariate-Poisson limit (:math:`g_{AB}
\to 1` for independent label classes). Same plug-in caveat as the
single-type estimator — prefer held-out :math:`\hat\lambda_A,
\hat\lambda_B`.
"""
_validate_edge_correction(edge_correction)
pts_A = _as_points(points_A)
pts_B = _as_points(points_B)
if pts_A.shape[1] != 2 or pts_B.shape[1] != 2:
raise ValueError("cross_pair_correlation is implemented for d=2 (planar).")
r_grid = np.asarray(r_grid, dtype=float)
if pts_A.shape[0] == 0 or pts_B.shape[0] == 0:
return np.zeros_like(r_grid)
if domain is None:
domain = _cross_default_domain(pts_A, pts_B)
area = _domain_measure(domain)
if bw is None:
n_total = pts_A.shape[0] + pts_B.shape[0]
bw = float(0.1 / np.sqrt(n_total / area)) if area > 0 else 0.05
bw = max(bw, float(np.mean(np.diff(r_grid))) if len(r_grid) > 1 else 0.05)
lam_A = _intensity_at(lambda_A, pts_A)
lam_B = _intensity_at(lambda_B, pts_B)
d = cdist(pts_A, pts_B)
wgt = 1.0 / (lam_A[:, None] * lam_B[None, :])
if edge_correction == "epanechnikov":
g = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
ring = 2.0 * np.pi * r * area
if ring <= 0:
g[k] = 0.0
continue
ker = epanechnikov((d - r) / bw) / bw
g[k] = float(np.sum(ker * wgt)) / ring
return g
_ensure_rectangular_domain(domain, edge_correction)
if edge_correction == "isotropic":
g = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
ring = 2.0 * np.pi * r * area
if ring <= 0:
g[k] = 0.0
continue
ker = epanechnikov((d - r) / bw) / bw
w_iso = np.zeros_like(d)
ii, jj = np.where(ker > 0)
for m in range(ii.size):
rij = d[ii[m], jj[m]]
fi = frac_disc_in_rect(pts_A[ii[m]], rij, domain)
fj = frac_disc_in_rect(pts_B[jj[m]], rij, domain)
w_iso[ii[m], jj[m]] = 0.5 * (1.0 / fi + 1.0 / fj)
g[k] = float(np.sum(ker * wgt * w_iso)) / ring
return g
if edge_correction == "translation":
offs = pts_A[:, None, :] - pts_B[None, :, :]
(xlo, xhi), (ylo, yhi) = domain
Lx, Ly = float(xhi - xlo), float(yhi - ylo)
ax = np.abs(offs[..., 0])
ay = np.abs(offs[..., 1])
overlap = np.maximum(Lx - ax, 0.0) * np.maximum(Ly - ay, 0.0)
with np.errstate(divide="ignore", invalid="ignore"):
w_tr = np.where(overlap > 0, (Lx * Ly) / overlap, np.inf)
g = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
ring = 2.0 * np.pi * r * area
if ring <= 0:
g[k] = 0.0
continue
ker = epanechnikov((d - r) / bw) / bw
g[k] = float(np.sum(ker * wgt * w_tr)) / ring
return g
# border
g = np.empty_like(r_grid)
for k, r in enumerate(r_grid):
ring = 2.0 * np.pi * r * area
if ring <= 0:
g[k] = 0.0
continue
usable = border_usable_mask(pts_A, domain, r)
if not usable.any():
g[k] = np.nan
continue
(xlo, xhi), (ylo, yhi) = domain
eff_area = max(0.0, (xhi - xlo) - 2 * r) * max(0.0, (yhi - ylo) - 2 * r)
if eff_area <= 0:
g[k] = np.nan
continue
ring_eff = 2.0 * np.pi * r * eff_area
ker = epanechnikov((d[usable, :] - r) / bw) / bw
g[k] = float(np.sum(ker * wgt[usable, :])) / ring_eff
return g
__all__ = [
"pair_correlation",
"k_inhom",
"l_function",
"nearest_neighbour_FGJ",
"global_envelope",
"cross_k_inhom",
"cross_pair_correlation",
"EnvelopeResult",
]