nstat.extras.spatial.marked_gof

Marked / discrete-time-rescaling goodness-of-fit (pure NumPy/SciPy).

A Python-only companion to the discrete-time / marked goodness-of-fit framework. Its centre of gravity is the discrete-time-rescaling correction of Haslinger, Pipa & Brown (2010), which fixes a real bug in the naive KS test at finite bin width.

The problem

Continuous time-rescaling says the rescaled inter-event times \(z_j = \Lambda(t_j) - \Lambda(t_{j-1})\) are i.i.d. Exp(1), so a KS test against Exp(1) (equivalently Unif(0,1) after the CDF transform) is exact. But on discretized time with bin width \(\Delta\) and per-bin spike probabilities \(p_k\), the naive bin-summed compensator is not Exp(1) even for the true model — the discreteness biases the statistic by \(O(\bar\lambda\Delta)\), so the naive KS false-rejects a correct model.

The fix (the randomized PIT)

With i.i.d. \(r_j \sim \text{Unif}(0,1)\), the corrected variate

\[u_j = \Bigl[\prod_{k=k_{j-1}+1}^{k_j - 1}(1 - p_k)\Bigr] \,(1 - r_j\, p_{k_j})\]

is exactly Unif(0,1) under the true model. A KS test on \(\{u_j\}\) (averaged over several \(r_j\) draws, or compared to a Monte-Carlo reference band) has the correct size.

Mark spaces

Marked rescaling factors as (rescaled time) × (conditional mark). A finite mark space (counting reference measure → mark = channel) is the multivariate-rescaling case (Gerhard-Haslinger-Pipa 2011); a continuous mark space (Lebesgue reference measure → e.g. a waveform amplitude) is rescaled through the conditional mark CDF \(F(m \mid \cdot)\). A mis-specified \(p(m\mid\cdot)\) shows up as a mark-axis KS departure while the time axis still passes.

All functions are pure NumPy/SciPy.

References

Haslinger R, Pipa G, Brown E (2010). Discrete time rescaling theorem: determining goodness of fit for discrete time statistical models of neural spiking. Neural Computation 22(10):2477.
Gerhard F, Haslinger R, Pipa G (2011). Applying the multivariate time-rescaling theorem to neural population models. Neural Computation 23(6):1452.
Tao L et al. (2018). Marked time-rescaling for population goodness-of-fit.

Functions

`corrected_rescaled`(spike_bins, p_k[, rng])	Discrete-time-corrected uniform variates (Haslinger-Pipa-Brown 2010).
`marked_time_rescaling`(spike_bins, marks, p_k)	Discrete-time + marked goodness-of-fit in one call.
`multivariate_gof_with_coupling`(...[, ...])	Discrete-time per-channel GoF and population coupling test, in one call.
`multivariate_time_rescaling`(...[, n_draws, ...])	Per-channel discrete-time rescaling for a finite (channel) mark space.
`rescaled_acf`(u_rescaled, *[, n_lags])	Autocorrelation of the rescaled-time variates with a Bartlett band.
`uncorrected_rescaled`(spike_bins, p_k)	Naive bin-summed rescaled times \(1 - e^{-\sum p_k}\).

Classes

`CoupledMarkedGOFResult`(per_channel, population)	Bundled result of `multivariate_gof_with_coupling()`.
`MarkedGOFResult`(ks_uncorrected, ...)	Result of a discrete-time / marked goodness-of-fit run.
`RescaledACFResult`(lags, acf, band, inside_band)	Autocorrelation of the rescaled-time variates with a Bartlett band.

class nstat.extras.spatial.marked_gof.MarkedGOFResult(ks_uncorrected, ks_corrected, ks_mark, ks_band, inside_uncorrected, inside_corrected, inside_mark, n_events, u_corrected, u_uncorrected)[source]

Bases: object

Result of a discrete-time / marked goodness-of-fit run.

Variables:

ks_uncorrected (float) – KS statistic of the naive (bin-summed) rescaled times vs Unif(0,1) — biased upward at finite bin width; expected to false-reject.
ks_corrected (float) – KS statistic of the discrete-time-corrected times (averaged over the randomization draws), vs Unif(0,1) — should pass for the true model.
ks_mark (float | None) – KS statistic on the mark axis (rescaled marks through the conditional mark CDF), vs Unif(0,1). None if no marks given.
ks_band (float) – The two-sided KS critical value at alpha for n events (Kolmogorov asymptotic c(alpha)/sqrt(n)).
inside_uncorrected (bool) – Whether each statistic falls inside the band (< ks_band). inside_* True ⇒ the test does NOT reject.
inside_corrected (bool) – See inside_uncorrected.
inside_mark (bool | None) – See inside_uncorrected.
n_events (int) – Number of inter-event intervals used.
u_corrected (numpy.ndarray) – The corrected variates from the first randomization draw (for plotting the empirical CDF).
u_uncorrected (numpy.ndarray) – The uncorrected variates (for the side-by-side plot).

Parameters:

ks_uncorrected (float)
ks_corrected (float)
ks_mark (float | None)
ks_band (float)
inside_uncorrected (bool)
inside_corrected (bool)
inside_mark (bool | None)
n_events (int)
u_corrected (ndarray)
u_uncorrected (ndarray)

ks_uncorrected: float

ks_corrected: float

ks_mark: float | None

ks_band: float

inside_uncorrected: bool

inside_corrected: bool

inside_mark: bool | None

n_events: int

u_corrected: ndarray

u_uncorrected: ndarray

class nstat.extras.spatial.marked_gof.CoupledMarkedGOFResult(per_channel, population)[source]

Bases: object

Bundled result of multivariate_gof_with_coupling().

Variables:

per_channel (dict[int, nstat.extras.spatial.marked_gof.MarkedGOFResult]) – Mapping channel_index -> MarkedGOFResult returned by multivariate_time_rescaling() — the discrete-time-corrected KS test for each channel in isolation.
population (object) – PopulationTimeRescaleResult returned by nstat.population_time_rescale() — the ground-process KS and marked-region Pearson \(\chi^2\) over the joint \((\tau, k)\) space, which catches cross-channel coupling misfit invisible to the per-channel test.

Parameters:

per_channel (dict[int, MarkedGOFResult])
population (object)

per_channel: dict[int, MarkedGOFResult]

population: object

class nstat.extras.spatial.marked_gof.RescaledACFResult(lags, acf, band, inside_band)[source]

Bases: object

Autocorrelation of the rescaled-time variates with a Bartlett band.

The discrete-time-rescaling KS test of marked_time_rescaling() checks the marginal distribution of the rescaled variates (Unif(0,1)) but is blind to serial dependence (Brown et al. 2002). This dataclass carries the lag-1..``n_lags`` autocorrelation of the normal-score-transformed uniforms, \(z_j = \Phi^{-1}(u_j)\), together with the asymptotic \(\pm 1.96/\sqrt{n}\) Bartlett band (Andersen 1997; Truccolo et al. 2005) and a per-lag in-band flag.

Variables:

lags (numpy.ndarray) – Integer lags 1..n_lags, shape (n_lags,).
acf (numpy.ndarray) – Sample autocorrelation of \(z_j\) at each lag, shape (n_lags,).
band (float) – Two-sided \(\pm 1.96 / \sqrt{n}\) Bartlett band — a scalar, applied symmetrically about zero.
inside_band (numpy.ndarray) – Boolean (n_lags,) mask: True where |acf| < band.

Parameters:

lags (ndarray)
acf (ndarray)
band (float)
inside_band (ndarray)

lags: ndarray

acf: ndarray

band: float

inside_band: ndarray

nstat.extras.spatial.marked_gof.uncorrected_rescaled(spike_bins, p_k)[source]

Naive bin-summed rescaled times \(1 - e^{-\sum p_k}\).

This is the statistic that false-rejects at finite bin width — the bug the discrete-time correction fixes. Provided so a demo can show the uncorrected curve straying outside the KS band.

Parameters:

spike_bins (ndarray) – Sorted integer bin indices at which events occurred.
p_k (ndarray) – Per-bin spike (hazard) probabilities \(p_k = \lambda_k\Delta\), length = number of time bins.

Return type:

ndarray

nstat.extras.spatial.marked_gof.corrected_rescaled(spike_bins, p_k, rng=None)[source]

Discrete-time-corrected uniform variates (Haslinger-Pipa-Brown 2010).

\[u_j = \Bigl[\prod_{k=k_{j-1}+1}^{k_j-1}(1-p_k)\Bigr] (1 - r_j\, p_{k_j}), \quad r_j \sim \text{Unif}(0,1).\]

Exactly Unif(0,1) under the true model.

Parameters:

spike_bins (ndarray) – Sorted integer event-bin indices.
p_k (ndarray) – Per-bin spike probabilities.
rng (Generator | None) – Random generator for the per-interval randomization draw \(r_j\).

Return type:

ndarray

nstat.extras.spatial.marked_gof.marked_time_rescaling(spike_bins, marks, p_k, mark_cdf=None, *, decoded=None, n_draws=25, alpha=0.05, rng=None)[source]

Discrete-time + marked goodness-of-fit in one call.

Runs the time axis both the naive (uncorrected) way and the discrete-time-corrected way (Haslinger-Pipa-Brown 2010), so the bin-width bias is visible side by side, and — if marks and a conditional mark CDF are supplied — the mark axis via \(F(m \mid \cdot)\).

Typical call:

from nstat.extras.spatial import marked_gof
ks = marked_gof.marked_time_rescaling(sb, spike_m, p_k, x_hat)

where x_hat is the decoded position at each event (passed as decoded); the conditional mark CDF is then mark_cdf(m, x).

Parameters:

spike_bins (ndarray) – Sorted integer event-bin indices.
marks (ndarray | None) – Per-event mark values (continuous-mark case), or None to run only the time axis. For a finite (channel) mark space, see multivariate_time_rescaling().
p_k (ndarray) – Per-bin spike probabilities \(p_k = \lambda_k \Delta\).
mark_cdf – Callable giving the conditional mark CDF. Either mark_cdf(m) (mark-only) or mark_cdf(m, x) when decoded is supplied — the latter is the conditional mark CDF F(m | x_hat).
decoded – Optional per-event covariate (e.g. decoded position) forwarded as the second argument to mark_cdf.
n_draws (int) – Number of randomization draws to average the corrected KS over.
alpha (float) – KS test level (sets the band).
rng (Generator | None) – Random generator.

Return type:

MarkedGOFResult

nstat.extras.spatial.marked_gof.multivariate_time_rescaling(spike_bins_per_channel, p_k_per_channel, *, n_draws=25, alpha=0.05, rng=None)[source]

Per-channel discrete-time rescaling for a finite (channel) mark space.

The finite-mark / multivariate case (Gerhard-Haslinger-Pipa 2011): rescale each channel by its own compensator computed under the joint history. Returns one MarkedGOFResult per channel.

Note

A coupling-blind fit can pass every per-channel KS while the joint model is wrong (Gerhard-Haslinger-Pipa 2011, §4): per-channel passing is necessary, not sufficient. The population-level diagnostic that closes that gap is nstat.population_time_rescale() (Tao, Weber, Arai & Eden 2018); see multivariate_gof_with_coupling() for a single call that runs both tests on the same data.

Parameters:

spike_bins_per_channel – Sequence of per-channel sorted event-bin index arrays.
p_k_per_channel – Sequence of per-channel per-bin spike-probability arrays.
n_draws (int) – As in marked_time_rescaling().
alpha (float) – As in marked_time_rescaling().
rng (Generator | None) – As in marked_time_rescaling().

Returns:

Channel index → result.

Return type:

dict[int, MarkedGOFResult]