Per-spike diagnostics¶

The paper's method: HPD overlap, KL divergence and the rank-based predictive p-value for every spike, thresholds from a baseline period, and flagging.

events ¶

Per-event diagnostics for marked point-process observations.

Spike trains (and other event streams) are usually decoded with a marked point-process observation model: each event carries a discrete mark (for spike-sorted data, the identity of the unit that fired) and each mark has a state-dependent intensity. The functions here turn such a model into the per-event quantities that the distribution-level diagnostics in :mod:statespacecheck.state_consistency compare:

:func:event_likelihood normalizes one event's mark intensity over the state space, giving the single-event likelihood.
:func:predictive_mark_probabilities gives the predictive probability of each mark for the next event.
:func:event_weighted_predictive gives the state distribution of the next event, weighting the predictive distribution by the total event intensity.
:func:mark_predictive_pvalue evaluates the predictive check exactly by summing over the finite set of marks.
:func:event_diagnostics computes HPD overlap, KL divergence, and the exact predictive p-value for every event in a recording.
:func:baseline_threshold estimates a flagging threshold from a baseline (well-specified) sample of per-event values.
:func:flag_events applies the paper's flagging rule to every event.

Array conventions follow the rest of the package: state distributions are (n_events, ...) or (n_time, ...) where ... is one or more spatial axes, and mark intensity tables are (..., n_marks) over the same spatial axes.

Classes¶

EventDiagnostics ¶

Bases: NamedTuple

Per-event diagnostic values: HPD overlap, KL divergence and the p-value.

Returned by :func:event_diagnostics and :func:~statespacecheck.clusterless_event_diagnostics.

Attributes:

Name	Type	Description
`hpd_overlap`	`(ndarray, shape(n_events))`	HPD overlap between the predictive distribution and each event's likelihood. Low values indicate poor local fit.
`kl_divergence`	`(ndarray, shape(n_events))`	KL divergence from the predictive distribution to each event's likelihood. High values indicate poor local fit.
`predictive_pvalue`	`(ndarray, shape(n_events))`	Rank-based predictive p-value of each event's observed mark: exact from :func:`event_diagnostics` (see :func:`mark_predictive_pvalue`), a Monte Carlo estimate from :func:`~statespacecheck.clusterless_event_diagnostics`. Low values indicate poor local fit.
`likelihood`	`np.ndarray, shape (n_events, ...), or None`	Normalized single-event likelihood of each event over the state space. `None` unless `return_likelihood=True`.

EventFlags ¶

Bases: NamedTuple

Per-event flags returned by :func:flag_events.

Each field is a boolean array of shape (n_events,) marking events whose diagnostic is on the misfit side of its threshold, or None if no threshold was given for that diagnostic.

Attributes:

Name	Type	Description
`hpd_overlap`	`np.ndarray of bool, shape (n_events,), or None`	HPD overlap at or below its threshold.
`kl_divergence`	`np.ndarray of bool, shape (n_events,), or None`	KL divergence at or above its threshold.
`predictive_pvalue`	`np.ndarray of bool, shape (n_events,), or None`	Predictive p-value at or below its cutoff.

Functions:¶

event_likelihood ¶

event_likelihood(event_intensities: ArrayLike) -> DistributionArray

Normalize event intensities over the state space.

For a marked point-process observation model, the likelihood contribution of a single event with mark c is proportional to that mark's intensity lambda_c(x) (or expected count lambda_c(x) * dt). This function normalizes each row to sum to 1 over the spatial axes, giving the single-event likelihood that :func:~statespacecheck.hpd_overlap and :func:~statespacecheck.kl_divergence compare with the predictive distribution.

The Poisson exposure term exp(-sum_c lambda_c(x) * dt) is shared by every event in a time bin, so it is deliberately left out: attaching it to each event would count it once per event when a bin contains several. A common bin width dt cancels on normalization, so rates and expected counts give the same result.

Normalization is done in log space, so rows with tiny but nonzero intensities keep their shape ([1e-20, 2e-20, 4e-20] becomes [1/7, 2/7, 4/7]).

Parameters:

Name	Type	Description	Default
`event_intensities`	`(ndarray, shape(n_events, ...))`	Nonnegative state-dependent intensity (or expected count) of each event's mark, where `...` represents one or more spatial axes.	required

Returns:

Name	Type	Description
`likelihood`	`(ndarray, shape(n_events, ...))`	Single-event likelihood; each row sums to 1 over the spatial axes.

Raises:

Type	Description
`ValueError`	If the input has no spatial axis, contains negative or non-finite values, or has a row that is zero everywhere (no defined likelihood).

Examples:

>>> import numpy as np
>>> from statespacecheck import event_likelihood
>>> event_likelihood(np.array([[1.0, 2.0, 1.0]]))
array([[0.25, 0.5 , 0.25]])

Source code in src/statespacecheck/events.py

def event_likelihood(event_intensities: ArrayLike) -> DistributionArray:
    """Normalize event intensities over the state space.

    For a marked point-process observation model, the likelihood contribution
    of a single event with mark ``c`` is proportional to that mark's intensity
    ``lambda_c(x)`` (or expected count ``lambda_c(x) * dt``). This function
    normalizes each row to sum to 1 over the spatial axes, giving the
    single-event likelihood that :func:`~statespacecheck.hpd_overlap` and
    :func:`~statespacecheck.kl_divergence` compare with the predictive
    distribution.

    The Poisson exposure term ``exp(-sum_c lambda_c(x) * dt)`` is shared by
    every event in a time bin, so it is deliberately left out: attaching it to
    each event would count it once per event when a bin contains several. A
    common bin width ``dt`` cancels on normalization, so rates and expected
    counts give the same result.

    Normalization is done in log space, so rows with tiny but nonzero
    intensities keep their shape (``[1e-20, 2e-20, 4e-20]`` becomes
    ``[1/7, 2/7, 4/7]``).

    Parameters
    ----------
    event_intensities : np.ndarray, shape (n_events, ...)
        Nonnegative state-dependent intensity (or expected count) of each
        event's mark, where ``...`` represents one or more spatial axes.

    Returns
    -------
    likelihood : np.ndarray, shape (n_events, ...)
        Single-event likelihood; each row sums to 1 over the spatial axes.

    Raises
    ------
    ValueError
        If the input has no spatial axis, contains negative or non-finite
        values, or has a row that is zero everywhere (no defined likelihood).

    Examples
    --------
    >>> import numpy as np
    >>> from statespacecheck import event_likelihood
    >>> event_likelihood(np.array([[1.0, 2.0, 1.0]]))
    array([[0.25, 0.5 , 0.25]])
    """
    event_intensities = as_array(event_intensities, "event_intensities", dtype=np.float64)
    if event_intensities.ndim < 2 or np.prod(event_intensities.shape[1:]) == 0:
        msg = (
            "event_intensities must have shape (n_events, ...) with a non-empty spatial "
            f"axis; got shape {event_intensities.shape}"
        )
        raise ValueError(msg)
    if not np.all(np.isfinite(event_intensities)) or np.any(event_intensities < 0.0):
        msg = "event_intensities must contain only finite nonnegative values"
        raise ValueError(msg)
    flat = flatten_time_spatial(event_intensities)
    zero_rows = np.flatnonzero(~(flat > 0.0).any(axis=1))
    if zero_rows.size:
        msg = (
            "Cannot compute an event likelihood for rows that are zero everywhere; "
            f"row indices: {_first(zero_rows)}"
        )
        raise ValueError(msg)
    with np.errstate(divide="ignore"):
        log_intensity = np.log(flat)
    return _normalize_log(log_intensity).reshape(event_intensities.shape)

predictive_mark_probabilities ¶

predictive_mark_probabilities(state_dist: ArrayLike, mark_intensities: ArrayLike) -> DistributionArray

Compute the predictive probability of each mark for the next event.

Mark intensities are averaged over the state distribution and then normalized across marks:

q[c] = sum_x p[x] * lambda_c(x) / sum_d sum_x p[x] * lambda_d(x).

This is the mark distribution of a randomly selected event under the predictive distribution. Normalizing across marks at each state before averaging would drop the weighting by the state-dependent total event rate, and is only equivalent when that total is constant across states.

Parameters:

Name	Type	Description	Default
`state_dist`	`(ndarray, shape(n_events, ...))`	Predictive state distribution for each event, where `...` represents one or more spatial axes.	required
`mark_intensities`	`(ndarray, shape(..., n_marks))`	Nonnegative intensity (or expected count) of every mark at every state.	required

Returns:

Name	Type	Description
`mark_probabilities`	`(ndarray, shape(n_events, n_marks))`	Predictive mark probabilities; each row sums to 1.

Raises:

Type	Description
`ValueError`	If shapes are inconsistent, inputs are negative or non-finite, or a row has zero (or non-finite) total predictive event intensity, for which the mark distribution is undefined.

Examples:

>>> import numpy as np
>>> from statespacecheck import predictive_mark_probabilities
>>> state = np.array([[0.5, 0.5]])
>>> intensities = np.array([[1.0, 0.0], [1.0, 2.0]])  # (n_bins, n_marks)
>>> predictive_mark_probabilities(state, intensities)
array([[0.5, 0.5]])

Source code in src/statespacecheck/events.py

def predictive_mark_probabilities(
    state_dist: ArrayLike, mark_intensities: ArrayLike
) -> DistributionArray:
    """Compute the predictive probability of each mark for the next event.

    Mark intensities are averaged over the state distribution and then
    normalized across marks:

    ``q[c] = sum_x p[x] * lambda_c(x) / sum_d sum_x p[x] * lambda_d(x)``.

    This is the mark distribution of a randomly selected event under the
    predictive distribution. Normalizing across marks at each state before
    averaging would drop the weighting by the state-dependent total event
    rate, and is only equivalent when that total is constant across states.

    Parameters
    ----------
    state_dist : np.ndarray, shape (n_events, ...)
        Predictive state distribution for each event, where ``...`` represents
        one or more spatial axes.
    mark_intensities : np.ndarray, shape (..., n_marks)
        Nonnegative intensity (or expected count) of every mark at every state.

    Returns
    -------
    mark_probabilities : np.ndarray, shape (n_events, n_marks)
        Predictive mark probabilities; each row sums to 1.

    Raises
    ------
    ValueError
        If shapes are inconsistent, inputs are negative or non-finite, or a
        row has zero (or non-finite) total predictive event intensity, for
        which the mark distribution is undefined.

    Examples
    --------
    >>> import numpy as np
    >>> from statespacecheck import predictive_mark_probabilities
    >>> state = np.array([[0.5, 0.5]])
    >>> intensities = np.array([[1.0, 0.0], [1.0, 2.0]])  # (n_bins, n_marks)
    >>> predictive_mark_probabilities(state, intensities)
    array([[0.5, 0.5]])
    """
    return _mark_probabilities(state_dist, mark_intensities)[0]

event_weighted_predictive ¶

event_weighted_predictive(state_dist: ArrayLike, ground_intensity: ArrayLike) -> DistributionArray

Compute the state distribution of the next event.

P_event(x) = Lambda(x) P(x) / sum_u Lambda(u) P(u), where P is the predictive state distribution and Lambda the ground intensity, the total event intensity at each state. A randomly chosen event is more likely to come from states with a higher total event intensity, so the state of an event is distributed as P weighted by Lambda. With a constant, positive ground intensity it is the normalized P.

Parameters:

Name	Type	Description	Default
`state_dist`	`(ndarray, shape(n_events, ...))`	Predictive state distribution for each event, where `...` represents one or more spatial axes. Rows need not be normalized.	required
`ground_intensity`	`(ndarray, shape(...))`	Nonnegative total event intensity at every state, over the same spatial axes. For sorted marks with intensities `mark_intensities` of shape `(..., n_marks)`, this is `mark_intensities.sum(axis=-1)`.	required

Returns:

Name	Type	Description
`event_weighted`	`(ndarray, shape(n_events, ...))`	Event-weighted state distribution; each row sums to 1.

Raises:

Type	Description
`ValueError`	If shapes are inconsistent, inputs are negative or non-finite, or a row has zero total event intensity under the state distribution, for which the event's state is undefined.

mark_predictive_pvalue ¶

mark_predictive_pvalue(state_dist: ArrayLike, mark_intensities: ArrayLike, observed_marks: ArrayLike) -> DistributionArray

Exact predictive p-value of each event's observed mark.

With a finite set of marks, the predictive check can be evaluated exactly rather than by Monte Carlo (compare :func:~statespacecheck.predictive_pvalue). For each event, the p-value is the probability that a mark drawn from the predictive mark distribution q (:func:predictive_mark_probabilities) is no more probable than the observed mark:

p = sum_c q[c] * 1{q[c] <= q[observed]}.

Small values mean the observed mark was unexpected given the predictive state distribution. A large value means only that this mark passes the check: the least probable mark gets a small p-value on every event even under a correct model, and marks the prediction makes equally probable all get p = 1, however unbalanced their observed frequencies. A relative tolerance on the <= comparison, 16 * eps * n_bins times the observed mark's probability, absorbs floating-point rounding (each probability is a sum of n_bins nonnegative terms, accurate to about n_bins * eps of itself), so marks with equal predictive probability receive equal p-values across platforms and at any scale, and a mark more probable than the observed one by more than this tolerance is not counted. Probabilities stay accurate however small the intensities (an event near the smallest normal float, about 2.2e-308, is recomputed at a larger scale), and are resolved to about 1e-30 in absolute terms.

Parameters:

Name	Type	Description	Default
`state_dist`	`(ndarray, shape(n_events, ...))`	Predictive state distribution for each event, where `...` represents one or more spatial axes.	required
`mark_intensities`	`(ndarray, shape(..., n_marks))`	Nonnegative intensity (or expected count) of every mark at every state.	required
`observed_marks`	`(ndarray, shape(n_events))`	Integer index of the observed mark of each event.	required

Returns:

Name	Type	Description
`pvalue`	`(ndarray, shape(n_events))`	Exact predictive p-values in `[0, 1]`.

Raises:

Type	Description
`ValueError`	If shapes are inconsistent, marks are out of range, inputs are masked, or the predictive mark distribution is undefined (see :func:`predictive_mark_probabilities`).

Examples:

>>> import numpy as np
>>> from statespacecheck import mark_predictive_pvalue
>>> state = np.array([[1.0, 0.0], [1.0, 0.0]])
>>> intensities = np.array([[8.0, 1.0, 1.0], [1.0, 1.0, 8.0]])  # (n_bins, n_marks)
>>> mark_predictive_pvalue(state, intensities, np.array([0, 2]))
array([1. , 0.2])

Source code in src/statespacecheck/events.py

def mark_predictive_pvalue(
    state_dist: ArrayLike,
    mark_intensities: ArrayLike,
    observed_marks: ArrayLike,
) -> DistributionArray:
    """Exact predictive p-value of each event's observed mark.

    With a finite set of marks, the predictive check can be evaluated exactly
    rather than by Monte Carlo (compare :func:`~statespacecheck.predictive_pvalue`).
    For each event, the p-value is the probability that a mark drawn from the
    predictive mark distribution ``q`` (:func:`predictive_mark_probabilities`)
    is no more probable than the observed mark:

    ``p = sum_c q[c] * 1{q[c] <= q[observed]}``.

    Small values mean the observed mark was unexpected given the predictive
    state distribution. A large value means only that this mark passes the
    check: the least probable mark gets a small p-value on every event even
    under a correct model, and marks the prediction makes equally probable all
    get p = 1, however unbalanced their observed frequencies. A relative
    tolerance on the ``<=`` comparison, ``16 * eps * n_bins`` times the
    observed mark's probability, absorbs floating-point rounding (each
    probability is a sum of ``n_bins`` nonnegative terms, accurate to about
    ``n_bins * eps`` of itself), so marks with equal predictive probability
    receive equal p-values across platforms and at any scale, and a mark more
    probable than the observed one by more than this tolerance is not counted.
    Probabilities stay accurate however small the intensities (an event near the
    smallest normal float, about ``2.2e-308``, is recomputed at a larger scale),
    and are resolved to about ``1e-30`` in absolute terms.

    Parameters
    ----------
    state_dist : np.ndarray, shape (n_events, ...)
        Predictive state distribution for each event, where ``...`` represents
        one or more spatial axes.
    mark_intensities : np.ndarray, shape (..., n_marks)
        Nonnegative intensity (or expected count) of every mark at every state.
    observed_marks : np.ndarray, shape (n_events,)
        Integer index of the observed mark of each event.

    Returns
    -------
    pvalue : np.ndarray, shape (n_events,)
        Exact predictive p-values in ``[0, 1]``.

    Raises
    ------
    ValueError
        If shapes are inconsistent, marks are out of range, inputs are masked,
        or the predictive mark distribution is undefined (see
        :func:`predictive_mark_probabilities`).

    Examples
    --------
    >>> import numpy as np
    >>> from statespacecheck import mark_predictive_pvalue
    >>> state = np.array([[1.0, 0.0], [1.0, 0.0]])
    >>> intensities = np.array([[8.0, 1.0, 1.0], [1.0, 1.0, 8.0]])  # (n_bins, n_marks)
    >>> mark_predictive_pvalue(state, intensities, np.array([0, 2]))
    array([1. , 0.2])
    """
    mark_probabilities, total_intensity = _mark_probabilities(state_dist, mark_intensities)
    n_events, n_marks = mark_probabilities.shape
    marks = _validate_marks(observed_marks, n_marks, "observed_marks")
    if marks.shape[0] != n_events:
        msg = (
            f"observed_marks must have one entry per event ({n_events}); got {marks.shape[0]}"
        )
        raise ValueError(msg)
    n_bins = int(np.prod(np.shape(state_dist)[1:]))
    observed = mark_probabilities[np.arange(n_events), marks]
    # Relative to the observed probability: the rounding error of each probability
    # scales with the probability itself, so a tolerance scaled by a larger one
    # would count marks many times more probable than a rare observed mark. The
    # absolute term bounds rounding in the subnormal range, up to the smallest
    # subnormal whatever the value's size: n_bins such roundings per expected
    # intensity, in probability units (divided by the total), plus one for the
    # quotient. Rows not rescaled have totals of at least n_bins * tiny / eps, so
    # it is at most about 16 * eps**2 (~8e-31)
    finfo = np.finfo(mark_probabilities.dtype)
    relative_tolerance = float(finfo.eps * n_bins * 16)
    # smallest / total <= 1 (the total is positive), so this cannot overflow
    smallest = float(finfo.smallest_subnormal)
    subnormal_tolerance = 16 * (n_bins * (smallest / total_intensity[:, 0]) + smallest)
    no_more_probable = (
        mark_probabilities
        <= (observed * (1.0 + relative_tolerance) + subnormal_tolerance)[:, None]
    )
    pvalue: DistributionArray = (mark_probabilities * no_more_probable).sum(axis=1)
    # The sum can exceed one by a few ulps; clip only that representational error.
    np.minimum(pvalue, 1.0, out=pvalue)
    return pvalue

event_diagnostics ¶

event_diagnostics(predictive: ArrayLike, mark_intensities: ArrayLike, event_time_ind: ArrayLike, event_marks: ArrayLike, *, coverage: float = DEFAULT_COVERAGE, return_likelihood: bool = False, batch_size: int = DEFAULT_EVENT_BATCH_SIZE) -> EventDiagnostics

Compute HPD overlap, KL divergence, and predictive p-value for every event.

Each event is compared with the one-step predictive distribution of its time bin. The event's likelihood is its mark intensity normalized over the state space (:func:event_likelihood), and the predictive p-value is the exact finite-mark check (:func:mark_predictive_pvalue). Events are processed in batches to bound memory for long recordings.

Parameters:

Name	Type	Description	Default
`predictive`	`(ndarray, shape(n_time, ...))`	One-step predictive state distribution `p(x_t \| y_{1:t-1})` at each time bin, where `...` represents one or more spatial axes. Any other state distribution (for example a smoother) can be substituted.	required
`mark_intensities`	`(ndarray, shape(..., n_marks))`	Nonnegative intensity (or expected count per bin) of every mark at every state, for example each unit's place field.	required
`event_time_ind`	`(ndarray, shape(n_events))`	Time-bin index of each event. If several events share a time bin, list each one separately; all are compared with that bin's predictive distribution, so events of the same mark in the same bin receive identical diagnostics.	required
`event_marks`	`(ndarray, shape(n_events))`	Mark index of each event (for spike-sorted data, the unit that fired).	required
`coverage`	`float`	Coverage probability of the HPD regions.	`0.95`
`return_likelihood`	`bool`	If True, also return each event's normalized likelihood, shape `(n_events, ...)`.	`False`
`batch_size`	`int`	Number of events processed at once.	`50_000`

Returns:

Type	Description
`EventDiagnostics`	Per-event `hpd_overlap`, `kl_divergence`, and `predictive_pvalue` arrays of shape `(n_events,)`, plus `likelihood` if requested.

Raises:

Type	Description
`ValueError`	If shapes are inconsistent, indices are out of range, or inputs are negative, non-finite or masked; if an event's mark has zero intensity at every position; or if the predictive distribution of an event's time bin puts no mass where any mark has intensity (or the total overflows).

Examples:

>>> import numpy as np
>>> from statespacecheck import event_diagnostics
>>> predictive = np.array([[0.7, 0.2, 0.1], [0.1, 0.2, 0.7]])  # (n_time, n_bins)
>>> place_fields = np.array([[5.0, 0.1], [1.0, 1.0], [0.1, 5.0]])  # (n_bins, n_marks)
>>> result = event_diagnostics(
...     predictive, place_fields, np.array([0, 1]), np.array([0, 0])
... )
>>> result.predictive_pvalue.round(3)
array([1.   , 0.172])

Source code in src/statespacecheck/events.py

def event_diagnostics(
    predictive: ArrayLike,
    mark_intensities: ArrayLike,
    event_time_ind: ArrayLike,
    event_marks: ArrayLike,
    *,
    coverage: float = DEFAULT_COVERAGE,
    return_likelihood: bool = False,
    batch_size: int = DEFAULT_EVENT_BATCH_SIZE,
) -> EventDiagnostics:
    """Compute HPD overlap, KL divergence, and predictive p-value for every event.

    Each event is compared with the one-step predictive distribution of its
    time bin. The event's likelihood is its mark intensity normalized over the
    state space (:func:`event_likelihood`), and the predictive p-value is the
    exact finite-mark check (:func:`mark_predictive_pvalue`). Events are
    processed in batches to bound memory for long recordings.

    Parameters
    ----------
    predictive : np.ndarray, shape (n_time, ...)
        One-step predictive state distribution ``p(x_t | y_{1:t-1})`` at each
        time bin, where ``...`` represents one or more spatial axes. Any other
        state distribution (for example a smoother) can be substituted.
    mark_intensities : np.ndarray, shape (..., n_marks)
        Nonnegative intensity (or expected count per bin) of every mark at every
        state, for example each unit's place field.
    event_time_ind : np.ndarray, shape (n_events,)
        Time-bin index of each event. If several events share a time bin, list
        each one separately; all are compared with that bin's predictive
        distribution, so events of the same mark in the same bin receive
        identical diagnostics.
    event_marks : np.ndarray, shape (n_events,)
        Mark index of each event (for spike-sorted data, the unit that fired).
    coverage : float, default 0.95
        Coverage probability of the HPD regions.
    return_likelihood : bool, default False
        If True, also return each event's normalized likelihood,
        shape ``(n_events, ...)``.
    batch_size : int, default 50_000
        Number of events processed at once.

    Returns
    -------
    EventDiagnostics
        Per-event ``hpd_overlap``, ``kl_divergence``, and ``predictive_pvalue``
        arrays of shape ``(n_events,)``, plus ``likelihood`` if requested.

    Raises
    ------
    ValueError
        If shapes are inconsistent, indices are out of range, or inputs are
        negative, non-finite or masked; if an event's mark has zero intensity at every
        position; or if the predictive distribution of an event's time bin
        puts no mass where any mark has intensity (or the total overflows).

    Examples
    --------
    >>> import numpy as np
    >>> from statespacecheck import event_diagnostics
    >>> predictive = np.array([[0.7, 0.2, 0.1], [0.1, 0.2, 0.7]])  # (n_time, n_bins)
    >>> place_fields = np.array([[5.0, 0.1], [1.0, 1.0], [0.1, 5.0]])  # (n_bins, n_marks)
    >>> result = event_diagnostics(
    ...     predictive, place_fields, np.array([0, 1]), np.array([0, 0])
    ... )
    >>> result.predictive_pvalue.round(3)
    array([1.   , 0.172])
    """
    validate_coverage(coverage)
    if batch_size < 1:
        msg = f"batch_size must be at least 1; got {batch_size}"
        raise ValueError(msg)
    predictive = _validate_predictive(predictive)
    spatial_shape = predictive.shape[1:]
    rates = _flatten_mark_intensities(mark_intensities, spatial_shape)
    time_ind = _validate_time_indices(event_time_ind, predictive.shape[0])
    marks = _validate_marks(event_marks, rates.shape[1], "event_marks")
    if time_ind.shape != marks.shape:
        msg = (
            "event_time_ind and event_marks must have the same length; got "
            f"{time_ind.shape[0]} and {marks.shape[0]}"
        )
        raise ValueError(msg)
    predictive_flat = flatten_time_spatial(predictive)
    _check_event_inputs(predictive_flat, rates, time_ind, marks)

    n_events = time_ind.shape[0]
    event_hpd: DistributionArray = np.empty(n_events)
    event_kl: DistributionArray = np.empty(n_events)
    event_pvalue: DistributionArray = np.empty(n_events)
    likelihood: DistributionArray | None = (
        np.empty((n_events, rates.shape[0])) if return_likelihood else None
    )

    for start in range(0, n_events, batch_size):
        stop = min(start + batch_size, n_events)
        batch_marks = marks[start:stop]
        predictive_batch = predictive_flat[time_ind[start:stop]]
        rates_batch = rates[:, batch_marks].T
        likelihood_batch = event_likelihood(rates_batch)

        event_hpd[start:stop] = hpd_overlap(
            predictive_batch, likelihood_batch, coverage=coverage
        )
        event_kl[start:stop] = kl_divergence(predictive_batch, likelihood_batch)
        # Where the likelihood underflowed, its log gives the divergence
        rows = _underflowed_likelihood_rows(
            predictive_batch, likelihood_batch, rates_batch > 0.0
        )
        if rows.size:
            with np.errstate(divide="ignore"):
                event_kl[start + rows] = _log_space_kl(
                    predictive_batch[rows], np.log(rates_batch[rows])
                )
        event_pvalue[start:stop] = mark_predictive_pvalue(predictive_batch, rates, batch_marks)
        if likelihood is not None:
            likelihood[start:stop] = likelihood_batch

    return EventDiagnostics(
        hpd_overlap=event_hpd,
        kl_divergence=event_kl,
        predictive_pvalue=event_pvalue,
        likelihood=None
        if likelihood is None
        else likelihood.reshape(n_events, *spatial_shape),
    )

baseline_threshold ¶

baseline_threshold(baseline_values: ArrayLike, quantile: float) -> float

Estimate a flagging threshold from baseline per-event diagnostic values.

Returns the quantile of values pooled from a period (or simulation) where the model is believed to be well specified. Use a low quantile for diagnostics where small values indicate misfit (HPD overlap, e.g. 0.01) and a high quantile where large values do (KL divergence, e.g. 0.99), then flag events at or beyond the threshold (:func:flag_events). This is the paper's rule. NaN values are ignored.

Because the comparison is inclusive, every value tied at the threshold is flagged, and the flagged fraction of the baseline can exceed quantile. Ties are common: HPD overlap is exactly 1 when one region is nested in the other and 0 when they are disjoint, and KL divergence is exactly 0 when the prediction equals the likelihood. If every baseline HPD overlap is 1, the threshold is 1 and every event at 1 is flagged. Report the flagged fraction of the baseline with the threshold. It describes the baseline the threshold was estimated from, not the rate of false alarms to expect elsewhere.

KL divergence is +inf when the prediction and the likelihood have disjoint support. Such values are allowed: when the requested quantile falls among them the threshold is +inf, and then only infinite values are at or above it.

Parameters:

Name	Type	Description	Default
`baseline_values`	`np.ndarray, shape (n_values,) or any shape`	Baseline diagnostic values; flattened before the quantile is taken.	required
`quantile`	`float`	Quantile in `[0, 1]`.	required

Returns:

Name	Type	Description
`threshold`	`float`	The requested quantile (linear interpolation) of the baseline values.

Raises:

Type	Description
`TypeError`	If the baseline values are complex.
`ValueError`	If `quantile` is outside `[0, 1]`, the baseline contains `-inf`, has no finite values, or is a masked array (use NaN for values to leave out).

Examples:

>>> import numpy as np
>>> from statespacecheck import baseline_threshold
>>> baseline_threshold(np.arange(101.0), 0.99)
99.0

Source code in src/statespacecheck/events.py

def baseline_threshold(baseline_values: ArrayLike, quantile: float) -> float:
    """Estimate a flagging threshold from baseline per-event diagnostic values.

    Returns the ``quantile`` of values pooled from a period (or simulation)
    where the model is believed to be well specified. Use a low quantile for
    diagnostics where small values indicate misfit (HPD overlap, e.g. 0.01)
    and a high quantile where large values do (KL divergence, e.g. 0.99), then
    flag events at or beyond the threshold (:func:`flag_events`). This is the
    paper's rule. NaN values are ignored.

    Because the comparison is inclusive, every value tied at the threshold is
    flagged, and the flagged fraction of the baseline can exceed ``quantile``.
    Ties are common: HPD overlap is exactly 1 when one region is nested in the
    other and 0 when they are disjoint, and KL divergence is exactly 0 when the
    prediction equals the likelihood. If every baseline HPD overlap is 1, the
    threshold is 1 and every event at 1 is flagged. Report the flagged fraction of
    the baseline with the threshold. It describes the baseline the threshold was
    estimated from, not the rate of false alarms to expect elsewhere.

    KL divergence is ``+inf`` when the prediction and the likelihood have
    disjoint support. Such values are allowed: when the requested quantile
    falls among them the threshold is ``+inf``, and then only infinite values
    are at or above it.

    Parameters
    ----------
    baseline_values : np.ndarray, shape (n_values,) or any shape
        Baseline diagnostic values; flattened before the quantile is taken.
    quantile : float
        Quantile in ``[0, 1]``.

    Returns
    -------
    threshold : float
        The requested quantile (linear interpolation) of the baseline values.

    Raises
    ------
    TypeError
        If the baseline values are complex.
    ValueError
        If ``quantile`` is outside ``[0, 1]``, the baseline contains ``-inf``,
        has no finite values, or is a masked array (use NaN for values to leave
        out).

    Examples
    --------
    >>> import numpy as np
    >>> from statespacecheck import baseline_threshold
    >>> baseline_threshold(np.arange(101.0), 0.99)
    99.0
    """
    if not 0.0 <= quantile <= 1.0:
        msg = f"quantile must lie in [0, 1]; got {quantile}"
        raise ValueError(msg)
    values = as_array(
        baseline_values,
        "baseline_values",
        "Pass an ndarray with NaN for values to leave out",
        dtype=float,
    ).ravel()
    if np.any(np.isneginf(values)):
        msg = "baseline_values contains -inf; a threshold cannot be estimated"
        raise ValueError(msg)
    values = values[~np.isnan(values)]
    if not np.any(np.isfinite(values)):
        msg = "baseline_values contains no finite values; the threshold would be undefined"
        raise ValueError(msg)
    # np.quantile interpolates linearly between the order statistics around the
    # quantile's position; interpolating toward +inf gives nan even with zero
    # weight, so check the two order statistics first.
    lower = np.quantile(values, quantile, method="lower")
    higher = np.quantile(values, quantile, method="higher")
    if lower == higher:
        return float(lower)
    if np.isinf(higher):
        return float(np.inf)
    return float(np.quantile(values, quantile))

flag_events ¶

flag_events(diagnostics: EventDiagnostics, *, hpd_overlap_threshold: float | None = None, kl_divergence_threshold: float | None = None, pvalue_threshold: float | None = 0.05) -> EventFlags

Flag events whose diagnostics indicate poor local fit.

Applies the paper's rule to each event independently: an event is flagged when its HPD overlap is at or below hpd_overlap_threshold, its KL divergence is at or above kl_divergence_threshold, or its predictive p-value is at or below pvalue_threshold. Each diagnostic is flagged separately; NaN values are never flagged. Values tied at a threshold are all flagged (see :func:baseline_threshold).

Thresholds for HPD overlap and KL divergence depend on the model and the data, so they have no default. In its simulation the paper sets them from a period where the model is believed to fit, with :func:baseline_threshold (1st percentile of HPD overlap, 99th percentile of KL divergence); for real data without such a period it flags HPD overlap at a fixed 0.05 and sets no KL divergence cutoff. It flags p-values at a fixed 0.05. It recommends HPD overlap and the predictive p-value as the primary diagnostics and KL divergence as a reference, because KL divergence is also large for consistent events when the prediction is broad.

Parameters:

Name	Type	Description	Default
`diagnostics`	`EventDiagnostics`	Per-event diagnostics from :func:`event_diagnostics`.	required
`hpd_overlap_threshold`	`float`	Flag HPD overlap at or below this value. Default None (not flagged).	`None`
`kl_divergence_threshold`	`float`	Flag KL divergence at or above this value. Default None (not flagged).	`None`
`pvalue_threshold`	`float`	Flag predictive p-values at or below this value. Default 0.05; None skips the p-value.	`0.05`

Returns:

Type	Description
`EventFlags`	Boolean flags for each diagnostic that has a threshold, else None.

Raises:

Type	Description
`ValueError`	If a threshold is NaN.

Examples:

Thresholds from a baseline period, as in the paper's simulation:

>>> import numpy as np
>>> from statespacecheck import baseline_threshold, event_diagnostics, flag_events
>>> rng = np.random.default_rng(0)
>>> predictive = rng.dirichlet(np.ones(20), size=200)  # (n_time, n_bins)
>>> place_fields = rng.gamma(2.0, size=(20, 8))  # (n_bins, n_units)
>>> time_ind, units = rng.integers(0, 200, 500), rng.integers(0, 8, 500)
>>> diagnostics = event_diagnostics(predictive, place_fields, time_ind, units)
>>> baseline = time_ind < 100
>>> flags = flag_events(
...     diagnostics,
...     hpd_overlap_threshold=baseline_threshold(diagnostics.hpd_overlap[baseline], 0.01),
...     kl_divergence_threshold=baseline_threshold(
...         diagnostics.kl_divergence[baseline], 0.99
...     ),
... )
>>> flags.predictive_pvalue.shape
(500,)