Comparing distributions

HPD overlap and KL divergence between a state distribution and a likelihood, for each row; used per spike by event_diagnostics and clusterless_event_diagnostics, or per time bin with a whole-bin likelihood (an extension beyond the paper).

state_consistency

Compare a state distribution with a likelihood: HPD overlap and KL divergence.

Each row (a time bin, or an event) pairs a state distribution, such as the one-step predictive distribution, with a likelihood over the same states. :func:hpd_overlap asks whether the two are consistent (their high-probability regions overlap); :func:kl_divergence measures how different they are.

The paper applies both to each spike, with the spike's single-event likelihood; :func:~statespacecheck.event_diagnostics does this. Applying them to whole time bins with a whole-bin likelihood, which also includes the Poisson exposure term and silent units, is an extension beyond the paper.

Functions:

kl_divergence

kl_divergence(state_dist: ArrayLike, likelihood: ArrayLike) -> DistributionArray

Compute Kullback-Leibler divergence between state distribution and likelihood.

Measures how different the likelihood is from the state distribution at each time point, D(state_dist || likelihood). The divergence is large when the two put their mass in different places, but also when the state distribution is broad relative to a consistent likelihood, so the paper uses it as a reference alongside :func:hpd_overlap and the predictive p-value.

Parameters:

Name	Type	Description	Default
state_dist	(ndarray, shape(n_time, ...))	State probability distributions over position at each time point where ... represents arbitrary spatial dimensions. Can be either one-step predictive distribution or smoother output. Non-negative values (NaN allowed to mark invalid bins). Automatically normalized over valid (non-NaN) bins.	required
likelihood	(ndarray, shape(n_time, ...))	Likelihood distributions at each time point. This is the likelihood p(y_t | x_t) across spatial positions. Non-negative values (NaN allowed to mark invalid bins). Automatically normalized over valid (non-NaN) bins. Must have same shape as state_dist.	required

Returns:

Name	Type	Description
kl_divergence	(ndarray, shape(n_time))	Kullback-Leibler divergence D_KL(state_dist || likelihood) at each time point. Values are non-negative, with 0 indicating identical distributions.

Raises:

Type	Description
ValueError	If state_dist and likelihood have different shapes, if distributions contain negative values, or if an input is a masked array (mark bins to exclude with NaN).
TypeError	If an input is complex.

Examples:

>>> import numpy as np
>>> from statespacecheck import kl_divergence
>>> # Identical distributions have zero divergence; the likelihood in the
>>> # second time bin puts its mass where the state distribution does not
>>> state = np.array([[0.3, 0.4, 0.3], [0.3, 0.4, 0.3]])
>>> like = np.array([[0.3, 0.4, 0.3], [0.1, 0.2, 0.7]])
>>> kl_divergence(state, like).round(3)
array([0.   , 0.353])

See Also

hpd_overlap : Compute spatial overlap between HPD regions highest_density_region : Compute highest density region mask

Notes

The KL divergence is computed using scipy.stats.entropy with the formula: D_KL(P || Q) = sum(P * log(P / Q)) where P is the state distribution and Q is the likelihood.

Distributions are automatically normalized over valid (non-NaN) bins. NaN values mark invalid spatial bins (e.g., inaccessible locations); a bin that is NaN (or infinite) in either input is excluded from both, for normalization and for the divergence.

Time slices where distributions have no valid mass return inf for the divergence, as do slices where the state has mass at a bin whose likelihood is zero (disjoint supports). A likelihood that is positive but too small relative to its total to represent once normalized (below about 2.2e-308) is taken from its log instead, so its divergence stays finite and accurate.

Source code in src/statespacecheck/state_consistency.py

def kl_divergence(state_dist: ArrayLike, likelihood: ArrayLike) -> DistributionArray:
    """Compute Kullback-Leibler divergence between state distribution and likelihood.

    Measures how different the likelihood is from the state distribution at each
    time point, D(state_dist || likelihood). The divergence is large when the two
    put their mass in different places, but also when the state distribution is
    broad relative to a consistent likelihood, so the paper uses it as a reference
    alongside :func:`hpd_overlap` and the predictive p-value.

    Parameters
    ----------
    state_dist : np.ndarray, shape (n_time, ...)
        State probability distributions over position at each time point where
        ... represents arbitrary spatial dimensions.
        Can be either one-step predictive distribution or smoother output.
        Non-negative values (NaN allowed to mark invalid bins).
        Automatically normalized over valid (non-NaN) bins.
    likelihood : np.ndarray, shape (n_time, ...)
        Likelihood distributions at each time point. This is the
        likelihood p(y_t | x_t) across spatial positions.
        Non-negative values (NaN allowed to mark invalid bins).
        Automatically normalized over valid (non-NaN) bins.
        Must have same shape as state_dist.

    Returns
    -------
    kl_divergence : np.ndarray, shape (n_time,)
        Kullback-Leibler divergence D_KL(state_dist || likelihood) at each
        time point. Values are non-negative, with 0 indicating identical
        distributions.

    Raises
    ------
    ValueError
        If state_dist and likelihood have different shapes, if distributions
        contain negative values, or if an input is a masked array (mark bins to
        exclude with NaN).
    TypeError
        If an input is complex.

    Examples
    --------
    >>> import numpy as np
    >>> from statespacecheck import kl_divergence
    >>> # Identical distributions have zero divergence; the likelihood in the
    >>> # second time bin puts its mass where the state distribution does not
    >>> state = np.array([[0.3, 0.4, 0.3], [0.3, 0.4, 0.3]])
    >>> like = np.array([[0.3, 0.4, 0.3], [0.1, 0.2, 0.7]])
    >>> kl_divergence(state, like).round(3)
    array([0.   , 0.353])

    See Also
    --------
    hpd_overlap : Compute spatial overlap between HPD regions
    highest_density_region : Compute highest density region mask

    Notes
    -----
    The KL divergence is computed using scipy.stats.entropy with the formula:
    D_KL(P || Q) = sum(P * log(P / Q))
    where P is the state distribution and Q is the likelihood.

    Distributions are automatically normalized over valid (non-NaN) bins.
    NaN values mark invalid spatial bins (e.g., inaccessible locations); a bin
    that is NaN (or infinite) in either input is excluded from both, for
    normalization and for the divergence.

    Time slices where distributions have no valid mass return inf for the divergence,
    as do slices where the state has mass at a bin whose likelihood is zero
    (disjoint supports). A likelihood that is positive but too small relative to
    its total to represent once normalized (below about 2.2e-308) is taken from its
    log instead, so its divergence stays finite and accurate.

    """
    state, like = as_paired_arrays(state_dist, likelihood)
    divergence: DistributionArray = np.empty(state.shape[0])
    for rows in row_chunks(state.shape):
        divergence[rows] = _kl_divergence_rows(state[rows], like[rows])
    return divergence

hpd_overlap

hpd_overlap(state_dist: ArrayLike, likelihood: ArrayLike, *, coverage: float = DEFAULT_COVERAGE) -> DistributionArray

Compute overlap between HPD regions of state distribution and likelihood.

Measures the overlap between the highest probability-density (HPD) regions of the state distribution and the likelihood, as a fraction of the smaller region (the Szymkiewicz-Simpson overlap coefficient). It is 1 when one region lies inside the other, so a broad prediction and a precise, consistent likelihood score 1, and 0 when the regions are disjoint. A high overlap means the likelihood passes this check, not that the prediction is informative: a prediction spread over most of the state space overlaps almost any likelihood.

Parameters:

Name	Type	Description	Default
state_dist	(ndarray, shape(n_time, ...))	State probability distributions over position at each time point where ... represents arbitrary spatial dimensions. Can be either one-step predictive distribution or smoother output. Non-negative values (NaN allowed to mark invalid bins). Automatically normalized over valid (non-NaN) bins.	required
likelihood	(ndarray, shape(n_time, ...))	Likelihood distributions at each time point. This is the likelihood p(y_t | x_t) across spatial positions. Non-negative values (NaN allowed to mark invalid bins). Automatically normalized over valid (non-NaN) bins. Must have same shape as state_dist.	required
coverage	float	Coverage probability for the HPD regions. Must be between 0 and 1. Default is 0.95 for 95% HPD regions.	DEFAULT_COVERAGE

Returns:

Name	Type	Description
hpd_overlap	(ndarray, shape(n_time))	Proportion of overlap between the HPD regions of state_dist and likelihood at each time point. Values range from 0 (no overlap) to 1 (the smaller region lies entirely inside the larger).

Raises:

Type	Description
ValueError	If state_dist and likelihood have different shapes, if coverage is not in (0, 1), if distributions contain negative values, or if an input is a masked array (mark bins to exclude with NaN).
TypeError	If an input is complex.

Examples:

>>> import numpy as np
>>> from statespacecheck import hpd_overlap
>>> # 80% HPD regions: bins {0, 1} for the state distribution and {1, 2}
>>> # for the likelihood share one bin out of the smaller region's two
>>> state = np.array([[0.4, 0.4, 0.2, 0.0, 0.0]])
>>> like = np.array([[0.0, 0.4, 0.4, 0.2, 0.0]])
>>> hpd_overlap(state, like, coverage=0.8)
array([0.5])

See Also

kl_divergence : Measure information divergence between distributions highest_density_region : Compute highest density region mask

Notes

The overlap is computed as: overlap = intersection(HPD_state, HPD_like) / min(size(HPD_state), size(HPD_like))

where a region's size is its number of bins. This equals the paper's region volume when all bins have the same volume; on a nonuniform grid, resample to a uniform one first.

This normalization ensures that: - overlap = 1.0 when one region completely contains the other - overlap = 0.0 when regions don't overlap at all - Values are comparable even when HPD regions have different sizes

When either HPD region is empty (a row with no probability mass), the denominator is 0 and overlap is defined as 0, so such rows read as disagreement. Check for all-zero rows separately if they can occur.

Distributions are automatically normalized over valid (non-NaN) bins. NaN values mark invalid spatial bins (e.g., inaccessible locations); a bin that is NaN (or infinite) in either input is excluded from both HPD regions.

Source code in src/statespacecheck/state_consistency.py

def hpd_overlap(
    state_dist: ArrayLike,
    likelihood: ArrayLike,
    *,
    coverage: float = DEFAULT_COVERAGE,
) -> DistributionArray:
    """Compute overlap between HPD regions of state distribution and likelihood.

    Measures the overlap between the highest probability-density (HPD) regions of
    the state distribution and the likelihood, as a fraction of the smaller region
    (the Szymkiewicz-Simpson overlap coefficient). It is 1 when one region lies
    inside the other, so a broad prediction and a precise, consistent likelihood
    score 1, and 0 when the regions are disjoint. A high overlap means the
    likelihood passes this check, not that the prediction is informative: a
    prediction spread over most of the state space overlaps almost any
    likelihood.

    Parameters
    ----------
    state_dist : np.ndarray, shape (n_time, ...)
        State probability distributions over position at each time point where
        ... represents arbitrary spatial dimensions.
        Can be either one-step predictive distribution or smoother output.
        Non-negative values (NaN allowed to mark invalid bins).
        Automatically normalized over valid (non-NaN) bins.
    likelihood : np.ndarray, shape (n_time, ...)
        Likelihood distributions at each time point. This is the
        likelihood p(y_t | x_t) across spatial positions.
        Non-negative values (NaN allowed to mark invalid bins).
        Automatically normalized over valid (non-NaN) bins.
        Must have same shape as state_dist.
    coverage : float, optional
        Coverage probability for the HPD regions. Must be between 0 and 1.
        Default is 0.95 for 95% HPD regions.

    Returns
    -------
    hpd_overlap : np.ndarray, shape (n_time,)
        Proportion of overlap between the HPD regions of state_dist and
        likelihood at each time point. Values range from 0 (no overlap)
        to 1 (the smaller region lies entirely inside the larger).

    Raises
    ------
    ValueError
        If state_dist and likelihood have different shapes, if coverage
        is not in (0, 1), if distributions contain negative values, or if an
        input is a masked array (mark bins to exclude with NaN).
    TypeError
        If an input is complex.

    Examples
    --------
    >>> import numpy as np
    >>> from statespacecheck import hpd_overlap
    >>> # 80% HPD regions: bins {0, 1} for the state distribution and {1, 2}
    >>> # for the likelihood share one bin out of the smaller region's two
    >>> state = np.array([[0.4, 0.4, 0.2, 0.0, 0.0]])
    >>> like = np.array([[0.0, 0.4, 0.4, 0.2, 0.0]])
    >>> hpd_overlap(state, like, coverage=0.8)
    array([0.5])

    See Also
    --------
    kl_divergence : Measure information divergence between distributions
    highest_density_region : Compute highest density region mask

    Notes
    -----
    The overlap is computed as:
        overlap = intersection(HPD_state, HPD_like) / min(size(HPD_state), size(HPD_like))

    where a region's size is its number of bins. This equals the paper's
    region volume when all bins have the same volume; on a nonuniform grid,
    resample to a uniform one first.

    This normalization ensures that:
    - overlap = 1.0 when one region completely contains the other
    - overlap = 0.0 when regions don't overlap at all
    - Values are comparable even when HPD regions have different sizes

    When either HPD region is empty (a row with no probability mass), the
    denominator is 0 and overlap is defined as 0, so such rows read as
    disagreement. Check for all-zero rows separately if they can occur.

    Distributions are automatically normalized over valid (non-NaN) bins.
    NaN values mark invalid spatial bins (e.g., inaccessible locations); a bin
    that is NaN (or infinite) in either input is excluded from both HPD regions.

    """
    validate_coverage(coverage)
    state, like = as_paired_arrays(state_dist, likelihood)
    overlap: DistributionArray = np.empty(state.shape[0])
    for rows in row_chunks(state.shape):
        overlap[rows] = _hpd_overlap_rows(state[rows], like[rows], coverage)
    return overlap