Coverage for pybmc/sampling_utils.py: 96%
74 statements
« prev ^ index » next coverage.py v7.10.0, created at 2026-07-20 21:03 +0000
1import numpy as np
3from .error_models import VARIANCE_FLOOR
4from .rng import DEFAULT_SEED, get_rng
6#: Default seed for posterior predictive draws (subsampling of posterior
7#: samples and the noise added on top) when the caller does not supply
8#: one explicitly. This is the same seed that drives the package-wide
9#: generator used by the MCMC samplers (see :mod:`pybmc.rng`), so one
10#: constant governs all pybmc randomness. Callers that need independent
11#: draws across repeated calls (e.g. outer Monte Carlo loops) can pass
12#: their own ``seed``, or ``seed=None`` to draw from the shared
13#: package-wide stream instead.
14DEFAULT_PREDICTIVE_SEED = DEFAULT_SEED
17def coverage(percentiles, rndm_m, models_output, truth_column):
18 """
19 Calculates coverage percentages for credible intervals.
21 Args:
22 percentiles (list[int]): Percentiles to evaluate (e.g., `[5, 10, ..., 95]`).
23 rndm_m (numpy.ndarray): Posterior samples of predictions.
24 models_output (pandas.DataFrame): DataFrame containing true values.
25 truth_column (str): Name of column with true values.
27 Returns:
28 list[float]: Coverage percentages for each percentile.
29 """
30 # How often the model's credible intervals actually contain the true value
31 data_total = len(rndm_m.T) # Number of data points
32 M_evals = len(rndm_m) # Number of samples
33 data_true = models_output[truth_column].tolist()
35 coverage_results = []
37 for p in percentiles:
38 count_covered = 0
39 for i in range(data_total):
40 # Sort model evaluations for the i-th data point
41 sorted_evals = np.sort(rndm_m.T[i])
42 # Find indices for lower and upper bounds of the credible interval
43 lower_idx = int((0.5 - p / 200) * M_evals)
44 upper_idx = int((0.5 + p / 200) * M_evals) - 1
45 # Check if the true value y[i] is within this interval
46 if sorted_evals[lower_idx] <= data_true[i] <= sorted_evals[upper_idx]:
47 count_covered += 1
48 coverage_results.append(count_covered / data_total * 100)
50 return coverage_results
53def rndm_m_random_calculator(
54 filtered_model_predictions, samples, Vt_hat, seed=DEFAULT_PREDICTIVE_SEED
55):
56 """
57 Posterior predictive samples for the homoscedastic model.
59 The homoscedastic model is the constant-only special case of the
60 heteroscedastic predictive distribution, so this delegates to
61 `rndm_m_heteroscedastic_calculator` with a constant-only variance
62 basis.
64 Args:
65 filtered_model_predictions (numpy.ndarray): Model predictions.
66 samples (numpy.ndarray): Posterior samples `[beta, sigma^2]`.
67 Vt_hat (numpy.ndarray): Normalized right singular vectors.
68 seed (int | None, optional): Seed for the posterior predictive
69 draws (subsampling of `samples` and the noise added on top).
70 Defaults to `DEFAULT_PREDICTIVE_SEED`; pass None to draw
71 from the shared package-wide stream instead.
73 Returns:
74 tuple[numpy.ndarray, list[numpy.ndarray]]:
75 - `rndm_m` (numpy.ndarray): Posterior predictive samples.
76 - `[lower, median, upper]` (list[numpy.ndarray]): Credible interval arrays.
77 """
78 constant_basis = np.ones((filtered_model_predictions.shape[0], 1))
79 return rndm_m_heteroscedastic_calculator(
80 filtered_model_predictions, samples, Vt_hat, constant_basis, seed=seed
81 )
84def rndm_m_heteroscedastic_calculator(
85 filtered_model_predictions, samples, Vt_hat, variance_basis,
86 seed=DEFAULT_PREDICTIVE_SEED,
87):
88 """
89 Generates posterior predictive samples for any error model.
91 The noise added to each prediction has a per-point variance
92 ``sigma_i^2 = phi_i . theta`` where ``theta`` are the variance
93 parameters of each posterior draw. The homoscedastic model is the
94 special case of a constant-only basis (a single column of ones),
95 handled by the `rndm_m_random_calculator` convenience wrapper.
97 Args:
98 filtered_model_predictions (numpy.ndarray): Model predictions,
99 shape ``(n_points, n_models)``.
100 samples (numpy.ndarray): Posterior samples
101 ``[beta_1..beta_k, theta_1..theta_p]`` from
102 `gibbs_sampler_heteroscedastic`.
103 Vt_hat (numpy.ndarray): Normalized right singular vectors, shape
104 ``(k, n_models)``.
105 variance_basis (numpy.ndarray): Variance design matrix for the
106 prediction points, shape ``(n_points, p)``.
107 seed (int | None, optional): Seed for the posterior predictive
108 draws (subsampling of `samples` and the noise added on top).
109 Defaults to `DEFAULT_PREDICTIVE_SEED`; pass None to draw
110 from the shared package-wide stream instead.
112 Returns:
113 tuple[numpy.ndarray, list[numpy.ndarray]]:
114 - `rndm_m` (numpy.ndarray): Posterior predictive samples.
115 - `[lower, median, upper]` (list[numpy.ndarray]): 95% credible
116 interval bounds and median.
117 """
118 rng = get_rng(seed)
120 n_draws = min(10000, len(samples))
121 replace = len(samples) < 10000
122 theta_rand_selected = rng.choice(samples, n_draws, replace=replace)
124 n_components = Vt_hat.shape[0]
125 betas = theta_rand_selected[:, :n_components]
126 variance_params = theta_rand_selected[:, n_components:]
128 # Model weights and noiseless central predictions.
129 default_weights = np.full(Vt_hat.shape[1], 1 / Vt_hat.shape[1])
130 model_weights_random = betas @ Vt_hat + default_weights
131 yvals_central = model_weights_random @ filtered_model_predictions.T
133 # Per-draw, per-point variances (n_draws, n_points). The floor is
134 # needed here (unlike during sampling) because the normalized metrics
135 # of extrapolated points can be negative, so phi . theta can dip
136 # below zero even though every theta draw is positive.
137 sigma2 = variance_params @ variance_basis.T
138 sigma2 = np.maximum(sigma2, VARIANCE_FLOOR)
140 noise = rng.standard_normal(yvals_central.shape) * np.sqrt(sigma2)
141 rndm_m = yvals_central + noise
143 lower = np.percentile(rndm_m, 2.5, axis=0)
144 median = np.percentile(rndm_m, 50, axis=0)
145 upper = np.percentile(rndm_m, 97.5, axis=0)
147 return rndm_m, [lower, median, upper]
150def coverage_quality(percentiles, coverage_results):
151 """
152 Scalar calibration score: mean |empirical - nominal| coverage.
154 Args:
155 percentiles (array-like): Nominal credible-interval widths in
156 percent (as passed to `coverage`).
157 coverage_results (array-like): Empirical coverage in percent.
159 Returns:
160 float: Mean absolute deviation in percentage points (lower is
161 better; 0 is perfect calibration).
162 """
163 return float(
164 np.mean(np.abs(np.asarray(coverage_results) - np.asarray(percentiles)))
165 )
168def mace(rndm_m, y_true, quantile_levels=None):
169 """
170 Mean Absolute Calibration Error (MACE): a quantile-based calibration score.
172 Complementary to `coverage`/`coverage_quality`'s two-sided central-interval
173 view: for each one-sided quantile level ``q``, compares the empirical
174 fraction of true values at or below the predictive ``q``-th percentile
175 against ``q`` itself.
177 Args:
178 rndm_m (numpy.ndarray): Posterior predictive draws, shape
179 ``(n_samples, n_points)``.
180 y_true (array-like): True/observed values, shape ``(n_points,)``.
181 quantile_levels (array-like, optional): Quantile levels in percent,
182 0-100 (default: 5, 10, ..., 95).
184 Returns:
185 float: Mean absolute deviation between empirical and nominal
186 quantile coverage, in percentage points (lower is better; 0 is
187 perfect calibration).
188 """
189 y_true = np.asarray(y_true, dtype=float)
190 if quantile_levels is None:
191 quantile_levels = np.arange(5, 100, 5)
192 quantile_levels = np.asarray(quantile_levels, dtype=float)
194 predicted_quantiles = np.percentile(rndm_m, quantile_levels, axis=0)
195 empirical = np.mean(y_true[None, :] <= predicted_quantiles, axis=1) * 100.0
196 return float(np.mean(np.abs(empirical - quantile_levels)))
199def reduced_chi_square(rndm_m, y_true):
200 """
201 Reduced chi-squared statistic of the posterior predictive distribution.
203 ``chi^2_red = mean_i[ (y_true_i - mean_i)^2 / var_i ]`` where
204 ``mean_i``/``var_i`` are the posterior predictive mean/variance at
205 point ``i``. A value near 1 indicates well-calibrated predictive
206 uncertainties; > 1 means the intervals are too narrow (overconfident,
207 "underdispersed" in the language of `diagnose_coverage_shape`); < 1
208 means too wide ("overdispersed").
210 Args:
211 rndm_m (numpy.ndarray): Posterior predictive draws, shape
212 ``(n_samples, n_points)``.
213 y_true (array-like): True/observed values, shape ``(n_points,)``.
215 Returns:
216 float: Reduced chi-squared statistic.
218 Raises:
219 ValueError: If the posterior predictive variance is non-positive
220 at any point (a degenerate/zero-noise predictive distribution).
221 """
222 y_true = np.asarray(y_true, dtype=float)
223 pred_mean = np.mean(rndm_m, axis=0)
224 pred_var = np.var(rndm_m, axis=0)
225 if np.any(pred_var <= 0):
226 raise ValueError(
227 "Posterior predictive variance must be positive at every point."
228 )
229 return float(np.mean((y_true - pred_mean) ** 2 / pred_var))
232def diagnose_coverage_shape(
233 percentiles, coverage_results, bias_tolerance=5.0, balance_threshold=0.7
234):
235 """
236 Classifies credible intervals as under-/over-dispersed or calibrated.
238 Args:
239 percentiles (array-like): Nominal interval widths in percent.
240 coverage_results (array-like): Empirical coverage in percent.
241 bias_tolerance (float, optional): Mean absolute deviation
242 (percentage points) below which the model counts as well
243 calibrated (default: 5.0).
244 balance_threshold (float, optional): Fraction of intervals that
245 must lie consistently on one side of nominal to call the
246 direction (default: 0.7).
248 Returns:
249 dict: Keys ``'diagnosis'`` (``'well_calibrated'`` |
250 ``'underdispersed'`` | ``'overdispersed'`` | ``'mixed'``),
251 ``'mean_bias'``, ``'mean_abs_error'``, ``'frac_below'``,
252 ``'frac_above'`` and ``'residuals'``. Underdispersed means the
253 intervals are too narrow (overconfident); overdispersed too wide.
254 """
255 percentiles = np.asarray(percentiles, dtype=float)
256 coverage_results = np.asarray(coverage_results, dtype=float)
258 residuals = coverage_results - percentiles
259 mean_bias = float(np.mean(residuals))
260 mean_abs_error = float(np.mean(np.abs(residuals)))
261 frac_below = float(np.mean(residuals < 0))
262 frac_above = float(np.mean(residuals > 0))
264 if mean_abs_error <= bias_tolerance:
265 diagnosis = "well_calibrated"
266 elif mean_bias < -bias_tolerance and frac_below >= balance_threshold:
267 diagnosis = "underdispersed"
268 elif mean_bias > bias_tolerance and frac_above >= balance_threshold:
269 diagnosis = "overdispersed"
270 else:
271 diagnosis = "mixed"
273 return {
274 "diagnosis": diagnosis,
275 "mean_bias": mean_bias,
276 "mean_abs_error": mean_abs_error,
277 "frac_below": frac_below,
278 "frac_above": frac_above,
279 "residuals": residuals,
280 }