Highest-density regions¶

The regions HPD overlap compares.

highest_density ¶

Highest probability-density (HPD) regions of distributions over a grid.

Functions:¶

highest_density_region ¶

highest_density_region(distribution: ArrayLike, *, coverage: float = DEFAULT_COVERAGE) -> NDArray[bool_]

Compute boolean mask indicating highest density region membership.

Vectorized HPD mask for arrays shaped (n_time, *spatial). For each time t, includes all bins with value >= threshold_t, where threshold_t is chosen so cumulative mass >= coverage * total_t (when rounding keeps the cumulative mass just short of that, with coverage just below 1, the region ends at the last bin with mass).

Parameters:

Name	Type	Description	Default
`distribution`	`(ndarray, shape(n_time, ...))`	Probability distributions over position at each time point where ... represents arbitrary spatial dimensions.	required
`coverage`	`float`	Desired coverage probability for the highest density region. Must be between 0 and 1. Default is 0.95 for 95% coverage.	`DEFAULT_COVERAGE`

Returns:

Name	Type	Description
`isin_hd`	`(ndarray, shape(n_time, ...))`	Boolean mask indicating which positions are in the highest density region at each time point, matching input shape.

Raises:

Type	Description
`ValueError`	If coverage is not in the range (0, 1), if distribution is not at least 2-D with spatial bins, if it contains negative values, or if it is a masked array (mark bins to exclude with NaN).
`TypeError`	If an input is complex.

Examples:

>>> import numpy as np
>>> from statespacecheck import highest_density_region
>>> # Simple 1D example with peaked distribution
>>> distribution = np.array([[0.1, 0.6, 0.3], [0.2, 0.5, 0.3]])
>>> region = highest_density_region(distribution, coverage=0.9)
>>> region.shape
(2, 3)
>>> region.dtype
dtype('bool')

See Also

hpd_overlap : Compute overlap between HPD regions of two distributions kl_divergence : Measure information divergence between distributions

Notes

NaN and infinite values are ignored (treated as 0 mass).
If total mass at time t is 0, returns all-False for that t.
Works in unnormalized space to avoid numerical issues.
Vectorized within chunks of time points; the chunks bound memory.
The input is probability mass per bin, so the region is the highest probability-density region only when all bins have the same volume.
Uses >= threshold: all bins with value equal to cutoff are included.
Due to ties, actual coverage may slightly exceed requested coverage.
This ensures consistent behavior across equivalent distributions.

References

.. [1] https://stats.stackexchange.com/questions/240749/how-to-find-95-credible-interval

Source code in src/statespacecheck/highest_density.py

def highest_density_region(
    distribution: ArrayLike, *, coverage: float = DEFAULT_COVERAGE
) -> NDArray[np.bool_]:
    """Compute boolean mask indicating highest density region membership.

    Vectorized HPD mask for arrays shaped (n_time, *spatial). For each time t,
    includes all bins with value >= threshold_t, where threshold_t is chosen so
    cumulative mass >= coverage * total_t (when rounding keeps the cumulative mass
    just short of that, with coverage just below 1, the region ends at the last
    bin with mass).

    Parameters
    ----------
    distribution : np.ndarray, shape (n_time, ...)
        Probability distributions over position at each time point where
        ... represents arbitrary spatial dimensions.
    coverage : float, optional
        Desired coverage probability for the highest density region. Must be between 0 and 1.
        Default is 0.95 for 95% coverage.

    Returns
    -------
    isin_hd : np.ndarray, shape (n_time, ...)
        Boolean mask indicating which positions are in the highest density region at each
        time point, matching input shape.

    Raises
    ------
    ValueError
        If coverage is not in the range (0, 1), if distribution is not at least
        2-D with spatial bins, if it contains negative values, or if it is a
        masked array (mark bins to exclude with NaN).
    TypeError
        If an input is complex.

    Examples
    --------
    >>> import numpy as np
    >>> from statespacecheck import highest_density_region
    >>> # Simple 1D example with peaked distribution
    >>> distribution = np.array([[0.1, 0.6, 0.3], [0.2, 0.5, 0.3]])
    >>> region = highest_density_region(distribution, coverage=0.9)
    >>> region.shape
    (2, 3)
    >>> region.dtype
    dtype('bool')

    See Also
    --------
    hpd_overlap : Compute overlap between HPD regions of two distributions
    kl_divergence : Measure information divergence between distributions

    Notes
    -----
    - NaN and infinite values are ignored (treated as 0 mass).
    - If total mass at time t is 0, returns all-False for that t.
    - Works in unnormalized space to avoid numerical issues.
    - Vectorized within chunks of time points; the chunks bound memory.
    - The input is probability mass per bin, so the region is the highest
      probability-density region only when all bins have the same volume.
    - Uses `>=` threshold: all bins with value equal to cutoff are included.
    - Due to ties, actual coverage may slightly exceed requested coverage.
    - This ensures consistent behavior across equivalent distributions.

    References
    ----------
    .. [1] https://stats.stackexchange.com/questions/240749/how-to-find-95-credible-interval

    """
    validate_coverage(coverage)
    values = as_array(distribution, "distribution", EXCLUDE_WITH_NAN, dtype=float)
    if values.ndim < 2:
        # Raise the usual error, which explains the expected shape
        validate_distribution(values, name="distribution", min_ndim=2, allow_nan=True)
    isin_hd = np.empty(values.shape, dtype=bool)
    for rows in row_chunks(values.shape):
        isin_hd[rows] = _highest_density_region_rows(values[rows], coverage)
    return isin_hd