Risk: Static Conformal, Streaming ACI, EVT & Epistemic¶
Navigation:
Theory introduction: See the Intro
Related mathematical theory: Risk: conformal, ACI, EVT & epistemic
Implementation only, for the coverage guarantees, the ACI update and the EVT/posterior theory see the math theory.
One source of truth for conformal: the static engine safety.py::SafetyTAM holds the finite-sample quantile and p-value; the streaming loop and the wrappers build on it. The static (i.i.d.) engine is kept deliberately apart from the streaming, non-stationary adaptation in statistics/risk/aci.py.
The static finite-sample quantile (SafetyTAM)¶
def conformal_quantile(self, alpha: Optional[float] = None) -> float:
"""Finite-sample-corrected (1 - alpha) empirical quantile of the calibration scores."""
current_alpha = self.alpha_target if alpha is None else alpha
n = len(self.residuals_calib_)
q_level = np.clip((1 - current_alpha) * (1 + 1 / n), 0, 1)
return float(np.quantile(self.residuals_calib_, q_level))
Adaptive Conformal Inference (streaming, model-agnostic)¶
def update_risk_level(alpha_t: float, error: float, alpha_target: float, gamma: float) -> float:
"""One ACI integrator step: alpha_{t+1} = alpha_t + gamma * (alpha_target - error).
error is the realized miscoverage indicator (1.0 if the point fell outside the interval, else 0.0).
This is the raw, unbounded integrator of Gibbs & Candes: the level is intentionally not clamped here so
cumulative miscoverage under sustained drift accumulates. It is clamped only where used - at the quantile
call, via effective_alpha.
"""
return float(alpha_t + gamma * (alpha_target - error))
Conformalized distributional wrapper (CQR + Mondrian + ACI)¶
class ConformalDistributionalTAM:
"""Wrap a fitted distributional StaticTAM with split-conformal calibration, p-values and ACI."""
def __init__(self, distributional_model, alpha: float = 0.1, strata_col: Optional[str] = None):
"""distributional_model is a fitted dict-formula StaticTAM; alpha the target miscoverage; strata_col
an optional column for Mondrian (per-stratum) calibration."""
if not 0.0 < alpha < 1.0:
raise ValueError(f"alpha must be in (0, 1); got {alpha}")
self.model = distributional_model
self.alpha = float(alpha)
self.strata_col = strata_col
self._studentized_engine: Dict[object, SafetyTAM] = {}
self._cqr_engine: Dict[object, SafetyTAM] = {}
def _strata(self, data: pd.DataFrame) -> np.ndarray:
if self.strata_col is None:
return np.full(len(data), _ALL_STRATA, dtype=object)
return data[self.strata_col].to_numpy()
def _studentized_scores(self, data: pd.DataFrame) -> np.ndarray:
return np.abs(self.model.anomaly_score(data)["z_score"].to_numpy())
def _engine(self, table: Dict[object, SafetyTAM], stratum: object) -> SafetyTAM:
return table.get(stratum, table[_ALL_STRATA])
def calibrate(self, calibration_data: pd.DataFrame) -> "ConformalDistributionalTAM":
"""Fit one SafetyTAM per stratum on the CQR and studentized nonconformity scores."""
lower_quantile = self.model.predict_quantile(calibration_data, self.alpha / 2.0)
upper_quantile = self.model.predict_quantile(calibration_data, 1.0 - self.alpha / 2.0)
observed = calibration_data[self.model.target_col_].to_numpy()
cqr_scores = np.maximum(lower_quantile - observed, observed - upper_quantile)
studentized_scores = self._studentized_scores(calibration_data)
strata = self._strata(calibration_data)
for stratum in list(np.unique(strata)) + [_ALL_STRATA]:
mask = np.ones(len(strata), dtype=bool) if stratum is _ALL_STRATA else (strata == stratum)
self._studentized_engine[stratum] = SafetyTAM(self.alpha).calibrate_scores(studentized_scores[mask])
self._cqr_engine[stratum] = SafetyTAM(self.alpha).calibrate_scores(cqr_scores[mask])
return self
def predict_interval(self, data: pd.DataFrame) -> pd.DataFrame:
"""CQR interval with finite-sample 1 - alpha coverage (per-stratum widths when Mondrian)."""
lower_quantile = self.model.predict_quantile(data, self.alpha / 2.0)
upper_quantile = self.model.predict_quantile(data, 1.0 - self.alpha / 2.0)
widths = np.array([self._engine(self._cqr_engine, s).conformal_quantile(self.alpha)
for s in self._strata(data)])
return pd.DataFrame({"lower": lower_quantile - widths, "upper": upper_quantile + widths}, index=data.index)
def conformal_pvalue(self, data: pd.DataFrame) -> np.ndarray:
"""Distribution-free conformal p-value of each observation (small => anomalous)."""
scores = self._studentized_scores(data)
strata = self._strata(data)
pvalues = np.empty(len(data), dtype=float)
for stratum in np.unique(strata):
mask = strata == stratum
pvalues[mask] = self._engine(self._studentized_engine, stratum).pvalue(scores[mask])
return pvalues
def anomaly(self, data: pd.DataFrame) -> pd.DataFrame:
"""Conformal p-value, the anomaly flag at level alpha, and the calibrated interval."""
pvalue = self.conformal_pvalue(data)
interval = self.predict_interval(data)
return pd.DataFrame({
"conformal_pvalue": pvalue,
"is_anomaly": pvalue < self.alpha,
"lower": interval["lower"].to_numpy(),
"upper": interval["upper"].to_numpy(),
}, index=data.index)
def aci_intervals(self, ordered_data: pd.DataFrame, gamma: float = 0.02) -> pd.DataFrame:
"""Adaptive Conformal Inference: intervals whose level adapts online to hold coverage under drift.
A single global risk level alpha_t adapts row by row (via aci.update_risk_level), while the conformal
radius is drawn from the row's stratum engine (SafetyTAM.conformal_quantile). This is the
Mondrian-stratified, location-scale specialisation of the generic ACI loop.
"""
mu_hat, sigma_hat = self.model._mu_sigma(ordered_data)
observed = self.model._to_model_scale(ordered_data[self.model.target_col_].to_numpy())
strata = self._strata(ordered_data)
alpha_t = self.alpha
lowers, uppers, levels = [], [], []
for position in range(len(ordered_data)):
radius = self._engine(self._studentized_engine, strata[position]).conformal_quantile(
effective_alpha(alpha_t)
)
lower = mu_hat[position] - sigma_hat[position] * radius
upper = mu_hat[position] + sigma_hat[position] * radius
inside = lower <= observed[position] <= upper
lowers.append(np.exp(lower) if self.model._log_target_ else lower)
uppers.append(np.exp(upper) if self.model._log_target_ else upper)
levels.append(alpha_t)
alpha_t = update_risk_level(alpha_t, 0.0 if inside else 1.0, self.alpha, gamma)
return pd.DataFrame({"lower": lowers, "upper": uppers, "alpha_t": levels}, index=ordered_data.index)
Extreme Value Theory tail scoring¶
class GeneralizedParetoTail:
"""A peaks-over-threshold GPD model for the upper tail of a nonnegative magnitude sample."""
def __init__(self, threshold_quantile: float = 0.95, max_shape: float = 1.0):
"""threshold_quantile is the high threshold u (an empirical quantile); max_shape caps the fitted GPD
shape xi (GPD-MLE is sensitive to the top order statistics and xi >= 1 means an infinite-mean tail,
so the shape is clamped and a warning is emitted when the raw estimate reaches the cap)."""
if not 0.5 <= threshold_quantile < 1.0:
raise ValueError(f"threshold_quantile must be in [0.5, 1); got {threshold_quantile}")
self.threshold_quantile = float(threshold_quantile)
self.max_shape = float(max_shape)
self.threshold_: float = 0.0
self.shape_: float = 0.0
self.scale_: float = 1.0
self.exceedance_rate_: float = 0.0
def fit(self, magnitudes: np.ndarray) -> "GeneralizedParetoTail":
"""Fit the GPD (shape, scale) to threshold exceedances by maximum likelihood."""
values = np.asarray(magnitudes, dtype=float)
self.threshold_ = float(np.quantile(values, self.threshold_quantile))
exceedances = values[values > self.threshold_] - self.threshold_
if exceedances.size < 10:
raise ValueError("Too few threshold exceedances to fit a GPD; lower threshold_quantile.")
fitted_shape, _, self.scale_ = genpareto.fit(exceedances, floc=0.0)
if fitted_shape >= self.max_shape:
warnings.warn(
f"GPD shape xi={fitted_shape:.3f} reached the cap {self.max_shape} (an infinite-mean tail). "
"Clamping for stable tail scores - inspect for a dominating outlier or raise threshold_quantile.",
RuntimeWarning,
stacklevel=2,
)
self.shape_ = float(min(fitted_shape, self.max_shape))
self.exceedance_rate_ = float(np.mean(values > self.threshold_))
return self
def tail_probability(self, values: np.ndarray) -> np.ndarray:
"""P(M > value) from the fitted GPD above the threshold; 1 below it (not extreme)."""
excess = np.asarray(values, dtype=float) - self.threshold_
survival = genpareto.sf(np.clip(excess, 0.0, None), self.shape_, loc=0.0, scale=self.scale_)
return np.where(excess > 0.0, self.exceedance_rate_ * survival, 1.0)
def return_level(self, exceedance_probability: float) -> float:
"""The magnitude exceeded with probability exceedance_probability (must be below the rate)."""
ratio = exceedance_probability / self.exceedance_rate_
if abs(self.shape_) < 1e-6:
return self.threshold_ - self.scale_ * np.log(ratio)
return self.threshold_ + self.scale_ / self.shape_ * (ratio ** (-self.shape_) - 1.0)
def anomaly_score(self, values: np.ndarray) -> np.ndarray:
"""EVT surprisal -log P(M > value) - large for extreme observations."""
return -np.log(np.clip(self.tail_probability(values), _TINY, 1.0))
Epistemic (parameter) uncertainty¶
def posterior_prediction(
model,
train_data: pd.DataFrame,
new_data: pd.DataFrame,
level: float = 0.95,
) -> pd.DataFrame:
"""Epistemic (parameter-uncertainty) prediction interval for a fitted single-group StaticTAM.
Returns a frame indexed like new_data with the linear-predictor prediction, its epistemic standard error
epistemic_std, and the `level` credible band (lower, upper).
Note: this models epistemic (parameter) uncertainty only, with the noise as a single global scalar
sigma^2. On heteroscedastic data (what a distributional StaticTAM is built for) this global variance
under-covers in high-variance regions and over-covers in low-variance ones; pair it with a distributional
StaticTAM's sigma(x) for the aleatoric component. The posterior is solved through a Cholesky factor of
(Phi.T Phi + n S) - stabler than an explicit inverse, and it never materialises the dense inverse or the
full posterior covariance (only triangular solves against Phi_new.T).
"""
if model.coefficients_ is None:
raise RuntimeError("Model must be fitted before computing posterior uncertainty.")
x_train, y_train = _prepared_tensors(model, train_data, with_target=True)
if x_train.shape[0] != 1:
raise NotImplementedError("posterior_prediction currently supports single-group models only.")
design_train = model._build_design_matrix(x_train.to(torch.float64))
penalty = model._build_penalty_matrix().to(torch.float64)
coefficients = model.coefficients_.to(torch.float64)
n_samples = x_train.shape[1]
n_coeffs = design_train.shape[-1]
gram = design_train.mT @ design_train
regularized = gram + n_samples * penalty + 1e-6 * n_samples * torch.eye(n_coeffs, dtype=torch.float64)
# Factor (Phi.T Phi + nS) once; reuse via triangular solves so neither the dense inverse nor the full
# posterior covariance is ever materialised. cholesky_ex is non-throwing; on a rare non-PD factorisation
# add a stronger ridge and retry.
cholesky_factor, info = torch.linalg.cholesky_ex(regularized)
if bool((info != 0).any()):
regularized = regularized + 1e-3 * n_samples * torch.eye(n_coeffs, dtype=torch.float64)
cholesky_factor, _ = torch.linalg.cholesky_ex(regularized)
residual = y_train.to(torch.float64) - design_train @ coefficients
residual_sum_of_squares = float((residual ** 2).sum())
# Effective dof = tr((Phi.T Phi + nS)^-1 Phi.T Phi), via a Cholesky solve (no explicit inverse).
effective_dof = float(torch.einsum("gii->g", torch.cholesky_solve(gram, cholesky_factor)).sum())
noise_variance = residual_sum_of_squares / max(n_samples - effective_dof, 1.0)
x_new, _ = _prepared_tensors(model, new_data, with_target=False)
design_new = model._build_design_matrix(x_new.to(torch.float64))
prediction = (design_new @ coefficients).squeeze(-1).squeeze(0).cpu().numpy()
# diag(Phi_new (Phi.T Phi + nS)^-1 Phi_new.T) via a solve against Phi_new.T (no D x D covariance).
covariance_solved = torch.cholesky_solve(design_new.mT, cholesky_factor)
prediction_variance = (noise_variance * torch.einsum(
"gnd,gdn->gn", design_new, covariance_solved
)).clamp_min(0.0)
epistemic_std = prediction_variance.sqrt().squeeze(0).cpu().numpy()
z_value = float(stats.norm.ppf(0.5 + level / 2.0))
return pd.DataFrame({
"prediction": prediction,
"epistemic_std": epistemic_std,
"lower": prediction - z_value * epistemic_std,
"upper": prediction + z_value * epistemic_std,
}, index=new_data.index)