Predictive densities and Monte Carlo checks¶

An extension beyond the paper: predictive densities of whole time bins and a Monte Carlo predictive p-value with a user-supplied sampler.

predictive_checks ¶

Predictive densities and Monte Carlo predictive checks of whole time bins.

These are extensions beyond the paper. The paper's predictive check is the rank-based predictive p-value of each spike, computed exactly over the units by :func:~statespacecheck.mark_predictive_pvalue and :func:~statespacecheck.event_diagnostics. The functions here compute the predictive density of all observations in a time bin, and a Monte Carlo p-value from a user-supplied sampler.

Functions:¶

predictive_density ¶

predictive_density(state_dist: ArrayLike, observation_likelihood: ArrayLike) -> DistributionArray

Compute predictive density by integrating state dist with obs likelihood.

CRITICAL: This function normalizes state_dist ONLY, NOT likelihood. The likelihood p(y|x) is a likelihood function, not a distribution over x. Normalizing it over x would change its value and mask real model misfit.

Formula: f_predictive(y) = ∑_x p(x) * p(y|x)

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. Non-negative values (NaN allowed to mark invalid bins). Will be normalized over spatial dimensions (everything except time).	required
`observation_likelihood`	`(ndarray, shape(n_time, ...))`	Observation likelihood p(y\|x) of the observed data at each position. Non-negative values; a NaN bin is excluded from both inputs (the state is renormalized over the other bins). Not normalized over positions (unlike the `likelihood` of :func:`~statespacecheck.kl_divergence`): it is a function of x, not a distribution. Must have same shape as state_dist.	required

Returns:

Name	Type	Description
`predictive_density`	`(ndarray, shape(n_time))`	Predictive density at each time point.

Raises:

Type	Description
`ValueError`	If state_dist and observation_likelihood have different shapes, if they contain negative values, if observation_likelihood contains +inf, 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 predictive_density
>>> # Simple 1D example with unnormalized state
>>> state = np.array([[3.0, 4.0, 3.0]])  # Unnormalized (sums to 10, not 1)
>>> like = np.array([[2.0, 3.0, 1.0]])  # Likelihood values (not normalized)
>>> pred = predictive_density(state, like)
>>> pred.shape
(1,)

log_predictive_density ¶

log_predictive_density(state_dist: ArrayLike, observation_likelihood: ArrayLike | None = None, *, log_observation_likelihood: ArrayLike | None = None) -> DistributionArray

Compute log predictive density directly in log-space using logsumexp.

CRITICAL: This function normalizes state_dist ONLY, NOT likelihood. Computes log predictive density natively in log-space for numerical stability. DO NOT compute as np.log(predictive_density(...)) - this loses precision.

Formula: log f_predictive(y) = log ∑_x p(x) * p(y|x) = logsumexp(log p(x) + log p(y|x))

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. Non-negative values (NaN allowed to mark invalid bins). Will be normalized over spatial dimensions (everything except time).	required
`observation_likelihood`	`(ndarray, shape(n_time, ...))`	Observation likelihood p(y\|x) of the observed data at each position. Non-negative values; a NaN bin is excluded from both inputs (the state is renormalized over the other bins). Not normalized over positions. Must have same shape as state_dist. Exactly one of `observation_likelihood` or `log_observation_likelihood` must be provided.	`None`
`log_observation_likelihood`	`(ndarray, shape(n_time, ...))`	Log observation likelihood log p(y\|x), keyword-only. Passing it avoids an exp/log round-trip. -inf is a zero likelihood; a NaN bin is excluded from both inputs. Must have same shape as state_dist.	`None`

Returns:

Name	Type	Description
`log_predictive_density`	`(ndarray, shape(n_time))`	Log predictive density at each time point.

Raises:

Type	Description
`ValueError`	If neither or both of the likelihood arguments are provided, if shapes don't match, if distributions contain negative values, if the (log) observation likelihood contains +inf, 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 log_predictive_density
>>> state = np.array([[1.0, 1.0, 1.0]])
>>> like = np.array([[2.0, 3.0, 4.0]])
>>> log_pred = log_predictive_density(state, like)
>>> log_pred.shape
(1,)

>>> # From a log likelihood, for numerical stability
>>> log_pred2 = log_predictive_density(state, log_observation_likelihood=np.log(like))
>>> np.allclose(log_pred, log_pred2)
True

See Also

predictive_density : Compute predictive density in linear space kl_divergence : Measure information divergence between distributions hpd_overlap : Compute spatial overlap between HPD regions

Notes

This function computes log predictive density directly in log-space using scipy.special.logsumexp for numerical stability. This prevents underflow when working with very small probabilities or peaked distributions.

The computation is: log ∑_x p(x) * p(y|x) = logsumexp(log p(x) + log p(y|x))

Where: - p(x) is the state distribution (normalized to sum to 1) - p(y|x) is the observation likelihood (NOT normalized)

For users who already have the log likelihood, passing it via log_observation_likelihood avoids the exp/log round-trip and is more efficient and numerically stable.

Source code in src/statespacecheck/predictive_checks.py

def log_predictive_density(
    state_dist: ArrayLike,
    observation_likelihood: ArrayLike | None = None,
    *,
    log_observation_likelihood: ArrayLike | None = None,
) -> DistributionArray:
    """Compute log predictive density directly in log-space using logsumexp.

    CRITICAL: This function normalizes state_dist ONLY, NOT likelihood.
    Computes log predictive density natively in log-space for numerical stability.
    DO NOT compute as np.log(predictive_density(...)) - this loses precision.

    Formula: log f_predictive(y) = log ∑_x p(x) * p(y|x)
             = logsumexp(log p(x) + log p(y|x))

    Parameters
    ----------
    state_dist : np.ndarray, shape (n_time, ...)
        State probability distributions over position at each time point where
        ... represents arbitrary spatial dimensions.
        Non-negative values (NaN allowed to mark invalid bins).
        Will be normalized over spatial dimensions (everything except time).
    observation_likelihood : np.ndarray, shape (n_time, ...), optional
        Observation likelihood p(y|x) of the observed data at each position.
        Non-negative values; a NaN bin is excluded from both inputs (the state
        is renormalized over the other bins). Not normalized over positions.
        Must have same shape as state_dist.
        Exactly one of `observation_likelihood` or `log_observation_likelihood`
        must be provided.
    log_observation_likelihood : np.ndarray, shape (n_time, ...), optional
        Log observation likelihood log p(y|x), keyword-only. Passing it avoids
        an exp/log round-trip. -inf is a zero likelihood; a NaN bin is excluded
        from both inputs. Must have same shape as state_dist.

    Returns
    -------
    log_predictive_density : np.ndarray, shape (n_time,)
        Log predictive density at each time point.

    Raises
    ------
    ValueError
        If neither or both of the likelihood arguments are provided,
        if shapes don't match, if distributions contain negative values, if
        the (log) observation likelihood contains +inf, 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 log_predictive_density
    >>> state = np.array([[1.0, 1.0, 1.0]])
    >>> like = np.array([[2.0, 3.0, 4.0]])
    >>> log_pred = log_predictive_density(state, like)
    >>> log_pred.shape
    (1,)

    >>> # From a log likelihood, for numerical stability
    >>> log_pred2 = log_predictive_density(state, log_observation_likelihood=np.log(like))
    >>> np.allclose(log_pred, log_pred2)
    True

    See Also
    --------
    predictive_density : Compute predictive density in linear space
    kl_divergence : Measure information divergence between distributions
    hpd_overlap : Compute spatial overlap between HPD regions

    Notes
    -----
    This function computes log predictive density directly in log-space using
    scipy.special.logsumexp for numerical stability. This prevents underflow
    when working with very small probabilities or peaked distributions.

    The computation is:
        log ∑_x p(x) * p(y|x) = logsumexp(log p(x) + log p(y|x))

    Where:
    - p(x) is the state distribution (normalized to sum to 1)
    - p(y|x) is the observation likelihood (NOT normalized)

    For users who already have the log likelihood, passing it via
    `log_observation_likelihood` avoids the exp/log round-trip and is more
    efficient and numerically stable.
    """
    if (observation_likelihood is None) == (log_observation_likelihood is None):
        msg = (
            "Exactly one of 'observation_likelihood' or 'log_observation_likelihood' "
            "must be provided"
        )
        raise ValueError(msg)

    if observation_likelihood is not None:
        state, like = as_paired_arrays(
            state_dist, observation_likelihood, "observation_likelihood"
        )
    else:
        state = as_array(state_dist, "state_dist", EXCLUDE_WITH_NAN, dtype=float)
        # Validate the log likelihood manually (it's in log-space, can be negative!)
        like = as_array(
            log_observation_likelihood,
            "log_observation_likelihood",
            EXCLUDE_WITH_NAN,
            dtype=float,
        )
        if like.ndim < 2:
            msg = (
                f"log_observation_likelihood must be at least 2D with shape (n_time, ...), "
                f"got shape {like.shape}"
            )
            raise ValueError(msg)
        if state.ndim < 2:
            validate_distribution(state, name="state_dist", min_ndim=2)
        if like.shape != state.shape:
            msg = (
                f"state_dist and log_observation_likelihood must have same shape, "
                f"got {state.shape} vs {like.shape}"
            )
            raise ValueError(msg)

    return _by_chunks_warning_on_zero_rows(
        partial(_log_predictive_density_rows, is_log=observation_likelihood is None),
        state,
        like,
    )

predictive_pvalue ¶

predictive_pvalue(observed_log_pred: ArrayLike, sample_log_pred: Callable[[int], ArrayLike], *, n_samples: int = 1000) -> DistributionArray

Compute predictive p-value via Monte Carlo sampling.

Computes p-values for predictive checks by comparing observed log predictive densities to a distribution of simulated log predictive densities. The p-value at each time point is the proportion of simulated values that are less than or equal to the observed value.

This computes a Monte Carlo predictive check p-value for a user-supplied replicate-generating procedure. If the model is correct and the statistic is continuous, p-values should be approximately uniformly distributed; ties and discreteness can make them conservative. Systematic deviations indicate model misfit.

Parameters:

Name	Type	Description	Default
`observed_log_pred`	`(ndarray, shape(n_time))`	Observed log predictive densities for actual data. Must be 1-dimensional.	required
`sample_log_pred`	`callable`	Function that generates samples of log predictive densities under the model. Must accept a single integer argument `n_samples` and return an array of shape (n_samples, n_time) containing simulated log predictive densities. For reproducibility, use np.random.Generator with a fixed seed internally. Example: `lambda n: rng.normal(loc=model_mean, scale=model_std, size=(n, n_time))`	required
`n_samples`	`int`	Number of Monte Carlo samples to draw for p-value computation. Higher values give more accurate p-value estimates but take longer. Default is 1000.	`1000`

Returns:

Name	Type	Description
`p_values`	`(ndarray, shape(n_time))`	P-value at each time point, computed as the proportion of simulated log predictive densities <= observed value. Values range from 0 to 1; NaN where `observed_log_pred` is NaN.

Raises:

Type	Description
`ValueError`	If observed_log_pred is not 1-dimensional, if n_samples <= 0, if sample_log_pred returns an array with the wrong shape or containing NaN, or if observed_log_pred or the sampler's output is a masked array (use NaN for missing observations).
`TypeError`	If sample_log_pred is not callable, or an input is complex.

Examples:

>>> import numpy as np
>>> from statespacecheck import predictive_pvalue
>>> # Observed log predictive densities
>>> observed = np.array([-2.0, -1.5, -1.0])
>>> # Sampler with internal random state for reproducibility
>>> def sampler(n_samples):
...     rng = np.random.default_rng(42)  # Fixed seed for reproducibility
...     return rng.normal(loc=-1.5, scale=0.5, size=(n_samples, 3))
>>> # Monte Carlo estimates of the exact values 0.16, 0.5 and 0.84
>>> predictive_pvalue(observed, sampler, n_samples=1000).round(2)
array([0.16, 0.5 , 0.86])

See Also

log_predictive_density : Compute log predictive density for observed data predictive_density : Compute predictive density in linear space aggregate_over_period : Aggregate metrics over time periods

Notes

The p-value at time t is computed as: p_value[t] = (1 / n_samples) * sum(simulated[t] <= observed[t])

Interpretation: the statistic is a log predictive density, so a small p-value means the observed data were less probable than nearly all replicates, i.e. unexpected under the model. A p-value near 1 means the observation was among the most probable outcomes, which is good fit, not misfit. Flag small values (for example p <= 0.05, as :func:~statespacecheck.periods.flag_extreme_pvalues does).

The estimate is the fraction of n_samples replicates, so it is a multiple of 1 / n_samples and can be exactly 0; choose n_samples large enough to resolve the cutoff you use.

The p-value is r / B: the fraction of B replicates at most as probable as the observation. It estimates the predictive tail probability, with standard error sqrt(p (1 - p) / B). A rank test at a finite B would use (r + 1) / (B + 1) instead, computed from the returned p as (p * B + 1) / (B + 1). It is valid in finite samples when the observation and the replicates are exchangeable under the model, and possibly conservative: it is never below 1 / (B + 1), and its level equals alpha only when alpha (B + 1) is an integer and there are no ties.

The sampler function should: 1. Generate new data from the model 2. Compute log predictive density for each generated dataset 3. Return array of shape (n_samples, n_time) 4. Use np.random.Generator internally for reproducibility

For reproducible results, create your sampler with a fixed seed: rng = np.random.default_rng(42) sampler = lambda n: rng.normal(size=(n, n_time))

Source code in src/statespacecheck/predictive_checks.py

def predictive_pvalue(
    observed_log_pred: ArrayLike,
    sample_log_pred: Callable[[int], ArrayLike],
    *,
    n_samples: int = 1000,
) -> DistributionArray:
    """Compute predictive p-value via Monte Carlo sampling.

    Computes p-values for predictive checks by comparing observed log predictive
    densities to a distribution of simulated log predictive densities. The p-value
    at each time point is the proportion of simulated values that are less than
    or equal to the observed value.

    This computes a Monte Carlo predictive check p-value for a user-supplied
    replicate-generating procedure. If the model is correct and the statistic is
    continuous, p-values should be approximately uniformly distributed; ties and
    discreteness can make them conservative. Systematic deviations indicate model
    misfit.

    Parameters
    ----------
    observed_log_pred : np.ndarray, shape (n_time,)
        Observed log predictive densities for actual data.
        Must be 1-dimensional.
    sample_log_pred : callable
        Function that generates samples of log predictive densities under the model.
        Must accept a single integer argument `n_samples` and return an array of
        shape (n_samples, n_time) containing simulated log predictive densities.
        For reproducibility, use np.random.Generator with a fixed seed internally.
        Example: `lambda n: rng.normal(loc=model_mean, scale=model_std, size=(n, n_time))`
    n_samples : int, optional
        Number of Monte Carlo samples to draw for p-value computation.
        Higher values give more accurate p-value estimates but take longer.
        Default is 1000.

    Returns
    -------
    p_values : np.ndarray, shape (n_time,)
        P-value at each time point, computed as the proportion of simulated
        log predictive densities <= observed value.
        Values range from 0 to 1; NaN where ``observed_log_pred`` is NaN.

    Raises
    ------
    ValueError
        If observed_log_pred is not 1-dimensional, if n_samples <= 0,
        if sample_log_pred returns an array with the wrong shape or
        containing NaN, or if observed_log_pred or the sampler's output is a
        masked array (use NaN for missing observations).
    TypeError
        If sample_log_pred is not callable, or an input is complex.

    Examples
    --------
    >>> import numpy as np
    >>> from statespacecheck import predictive_pvalue
    >>> # Observed log predictive densities
    >>> observed = np.array([-2.0, -1.5, -1.0])
    >>> # Sampler with internal random state for reproducibility
    >>> def sampler(n_samples):
    ...     rng = np.random.default_rng(42)  # Fixed seed for reproducibility
    ...     return rng.normal(loc=-1.5, scale=0.5, size=(n_samples, 3))
    >>> # Monte Carlo estimates of the exact values 0.16, 0.5 and 0.84
    >>> predictive_pvalue(observed, sampler, n_samples=1000).round(2)
    array([0.16, 0.5 , 0.86])

    See Also
    --------
    log_predictive_density : Compute log predictive density for observed data
    predictive_density : Compute predictive density in linear space
    aggregate_over_period : Aggregate metrics over time periods

    Notes
    -----
    The p-value at time t is computed as:
        p_value[t] = (1 / n_samples) * sum(simulated[t] <= observed[t])

    Interpretation: the statistic is a log predictive density, so a small
    p-value means the observed data were less probable than nearly all
    replicates, i.e. unexpected under the model. A p-value near 1 means the
    observation was among the most probable outcomes, which is good fit, not
    misfit. Flag small values (for example ``p <= 0.05``, as
    :func:`~statespacecheck.periods.flag_extreme_pvalues` does).

    The estimate is the fraction of ``n_samples`` replicates, so it is a
    multiple of ``1 / n_samples`` and can be exactly 0; choose ``n_samples``
    large enough to resolve the cutoff you use.

    The p-value is ``r / B``: the fraction of ``B`` replicates at most as
    probable as the observation. It estimates the predictive tail probability,
    with standard error ``sqrt(p (1 - p) / B)``. A rank test at a finite ``B``
    would use ``(r + 1) / (B + 1)`` instead, computed from the returned ``p`` as
    ``(p * B + 1) / (B + 1)``. It is valid in finite samples when the observation
    and the replicates are exchangeable under the model, and possibly conservative:
    it is never below ``1 / (B + 1)``, and its level equals ``alpha`` only when
    ``alpha (B + 1)`` is an integer and there are no ties.

    The sampler function should:
    1. Generate new data from the model
    2. Compute log predictive density for each generated dataset
    3. Return array of shape (n_samples, n_time)
    4. Use np.random.Generator internally for reproducibility

    For reproducible results, create your sampler with a fixed seed:
        rng = np.random.default_rng(42)
        sampler = lambda n: rng.normal(size=(n, n_time))
    """
    # Validate observed_log_pred
    observed_arr = as_array(
        observed_log_pred,
        "observed_log_pred",
        "Pass an ndarray with NaN for missing values",
        dtype=float,
    )
    if observed_arr.ndim != 1:
        msg = (
            f"observed_log_pred must be 1-dimensional, "
            f"got {observed_arr.ndim}D array with shape {observed_arr.shape}"
        )
        raise ValueError(msg)

    n_time = observed_arr.shape[0]

    # Validate n_samples
    if n_samples <= 0:
        msg = f"n_samples must be positive, got {n_samples}"
        raise ValueError(msg)

    # Generate samples
    simulated = sample_log_pred(n_samples)

    # Validate shape of simulated samples
    simulated_arr = as_array(
        simulated, "sample_log_pred's output", "Return an ndarray", dtype=float
    )
    if simulated_arr.shape != (n_samples, n_time):
        msg = (
            f"sample_log_pred output must have shape (n_samples, n_time) = "
            f"({n_samples}, {n_time}), got shape {simulated_arr.shape}"
        )
        raise ValueError(msg)
    # A NaN sample compares False, which would silently pull the p-value toward 0
    # and read as misfit; it is a sampler error.
    nan_times = np.isnan(simulated_arr).any(axis=0)
    if nan_times.any():
        bad = np.flatnonzero(nan_times)
        msg = f"sample_log_pred returned NaN at time indices: {bad[:10].tolist()}"
        raise ValueError(msg)

    # Proportion of samples <= observed at each time: (n_samples, n_time) -> (n_time,).
    # A NaN observation gives a NaN p-value; +-inf follow the comparison
    # (-inf, impossible under the model, gives 0).
    mask = ~np.isnan(observed_arr)
    p_values: DistributionArray = np.full(n_time, np.nan)
    if np.any(mask):
        p_values[mask] = np.mean(simulated_arr[:, mask] <= observed_arr[mask], axis=0)
    return p_values