API: tools

offset

Logit-offset (central tendency) adjustment with audit trail.

Theory — including the King–Zeng / Elkan / Tasche equivalences and the uniqueness of the mode-B root: docs/concepts/offset.md.

References

King & Zeng (2001); Elkan (2001); Tasche (2013) — full records in the documentation.

AuditReport dataclass

AuditReport(delta: float, pre_mean: float, post_mean: float, timestamp: str, guardrails_before: GuardrailReport, guardrails_after: GuardrailReport)

Pre/post record of a logit-offset application, for validators.

ATTRIBUTE	DESCRIPTION
delta	Applied log-odds shift. TYPE: float
pre_mean, post_mean	Portfolio mean probability before and after the shift. TYPE: float
timestamp	ISO-8601 UTC time at which the offset was fitted. TYPE: str
guardrails_before, guardrails_after	The three-flag calibration health summary on the input and output probabilities. TYPE: GuardrailReport

LogitOffset

LogitOffset(delta: float | None = None, target_mean: float | None = None)

Uniform log-odds shift: p' = sigma(logit(p) + delta).

Mode A takes delta explicitly; mode B takes target_mean and solves mean(p') = target_mean for delta by bisection — the portfolio mean is strictly increasing in delta, so the root is unique (stated in the offset chapter and unit-tested). Exactly one of the two arguments must be given.

The offset is deliberately not folded into any calibrator's parameters: CalibratedModel.offset_to appends it as a separate, inspectable pipeline stage.

ATTRIBUTE	DESCRIPTION
delta_	Fitted (or given) shift in log-odds. TYPE: float
pre_mean_, post_mean_	Portfolio mean before and after, recorded at fit time. TYPE: float
timestamp_	ISO-8601 UTC fit time — part of the audit trail. TYPE: str

Source code in src/probcal/offset.py

def __init__(self, delta: float | None = None, target_mean: float | None = None) -> None:
    self.delta = delta
    self.target_mean = target_mean

affine_logit_coeffs_ property

affine_logit_coeffs_: tuple[float, float]

(1, delta): the offset is affine on the logit scale.

fit

fit(p: object, sample_weight: object = None) -> Self

Fix delta (mode A) or solve it against the target mean (mode B).

PARAMETER	DESCRIPTION
p	Current calibrated probabilities of the portfolio. TYPE: array_like
sample_weight	Weights for the portfolio mean. TYPE: array_like or None DEFAULT: None

Source code in src/probcal/offset.py

def fit(self, p: object, sample_weight: object = None) -> Self:
    """Fix ``delta`` (mode A) or solve it against the target mean (mode B).

    Parameters
    ----------
    p : array_like
        Current calibrated probabilities of the portfolio.
    sample_weight : array_like or None
        Weights for the portfolio mean.
    """
    if (self.delta is None) == (self.target_mean is None):
        raise ValueError("LogitOffset: give exactly one of delta or target_mean")
    p_arr = validate_scores(p, name="p")
    w = validate_weights(sample_weight, len(p_arr))
    z = logit(p_arr)
    self.pre_mean_ = float(np.average(p_arr, weights=w))
    if self.delta is not None:
        self.delta_ = float(self.delta)
    else:
        target = float(self.target_mean)  # type: ignore[arg-type]
        if not 0.0 < target < 1.0:
            raise ValueError("target_mean must lie in (0, 1)")

        def gap(d: float) -> float:
            return float(np.average(expit(z + d), weights=w)) - target

        self.delta_ = bisect(gap, -_DELTA_BRACKET, _DELTA_BRACKET, tol=1e-14)
    self.post_mean_ = float(np.average(expit(z + self.delta_), weights=w))
    self.timestamp_ = datetime.now(UTC).isoformat(timespec="seconds")
    self.fitted_ = True
    return self

transform

transform(p: object) -> ndarray

Apply the fitted shift to probabilities.

Source code in src/probcal/offset.py

def transform(self, p: object) -> np.ndarray:
    """Apply the fitted shift to probabilities."""
    if not self.fitted_:
        raise RuntimeError("LogitOffset is not fitted; call fit() first")
    return expit(logit(validate_scores(p, name="p")) + self.delta_)

interpret

interpret() -> Interpretation

Read delta in log-odds, odds-factor, and central-tendency terms.

Source code in src/probcal/offset.py

def interpret(self) -> Interpretation:
    """Read delta in log-odds, odds-factor, and central-tendency terms."""
    if not self.fitted_:
        raise RuntimeError("LogitOffset is not fitted; call fit() first")
    return Interpretation(
        method=type(self).__name__,
        param_names=("delta",),
        param_values=(self.delta_,),
        messages=(
            f"delta = {self.delta_:+.4f} log-odds: every observation's odds are "
            f"multiplied by exp(delta) = {np.exp(self.delta_):.4f} uniformly",
            f"portfolio mean re-anchored from {self.pre_mean_:.5f} to "
            f"{self.post_mean_:.5f} (credit-risk central tendency adjustment)",
            "equivalent to King-Zeng prior correction and to Elkan's base-rate "
            "adjustment (see the offset chapter for the derivations)",
            "ranking is untouched: the shift is strictly increasing",
        ),
    )

interval_inverse

interval_inverse(lo: float, hi: float, *, space: str = 'probability', buffer_logit: float = 0.0) -> tuple[float, float]

Closed-form preimage: subtract delta on the logit scale.

Same protocol as BaseCalibrator.interval_inverse; the offset's output range is the full unit interval, so only a crossed buffer can make a target unattainable.

Source code in src/probcal/offset.py

def interval_inverse(
    self,
    lo: float,
    hi: float,
    *,
    space: str = "probability",
    buffer_logit: float = 0.0,
) -> tuple[float, float]:
    """Closed-form preimage: subtract delta on the logit scale.

    Same protocol as ``BaseCalibrator.interval_inverse``; the offset's
    output range is the full unit interval, so only a crossed buffer can
    make a target unattainable.
    """
    if not self.fitted_:
        raise RuntimeError("LogitOffset is not fitted; call fit() first")
    if not 0.0 <= lo <= hi <= 1.0:
        raise ValueError(f"need 0 <= lo <= hi <= 1, got lo={lo}, hi={hi}")
    from .base import UnattainableTargetError

    lo_b, hi_b = float(lo), float(hi)
    if buffer_logit > 0.0:
        if lo > 0.0:
            lo_b = float(expit(np.array([logit(np.array([lo]))[0] + buffer_logit]))[0])
        if hi < 1.0:
            hi_b = float(expit(np.array([logit(np.array([hi]))[0] - buffer_logit]))[0])
        if lo_b > hi_b:
            raise UnattainableTargetError(
                f"buffer_logit={buffer_logit} empties the calibrated interval [{lo}, {hi}]"
            )
    lo_z = -np.inf if lo_b <= 0.0 else float(logit(np.array([lo_b]))[0]) - self.delta_
    hi_z = np.inf if hi_b >= 1.0 else float(logit(np.array([hi_b]))[0]) - self.delta_
    if space == "logit":
        return lo_z, hi_z
    raw_lo = 0.0 if np.isneginf(lo_z) else float(expit(np.array([lo_z]))[0])
    raw_hi = 1.0 if np.isposinf(hi_z) else float(expit(np.array([hi_z]))[0])
    return raw_lo, raw_hi

audit_report

audit_report(y: object, p: object, *, sample_weight: object = None) -> AuditReport

Pre/post guardrail comparison for the validator's one-table view.

Source code in src/probcal/offset.py

def audit_report(self, y: object, p: object, *, sample_weight: object = None) -> AuditReport:
    """Pre/post guardrail comparison for the validator's one-table view."""
    if not self.fitted_:
        raise RuntimeError("LogitOffset is not fitted; call fit() first")
    before = calibration_guardrails(y, p, sample_weight=sample_weight)
    after = calibration_guardrails(y, self.transform(p), sample_weight=sample_weight)
    return AuditReport(
        delta=self.delta_,
        pre_mean=self.pre_mean_,
        post_mean=self.post_mean_,
        timestamp=self.timestamp_,
        guardrails_before=before,
        guardrails_after=after,
    )

wrapper

CalibratedModel: model-level wrapper with prefit and cross-validation flows.

Theory of the flows (why prefit is the credit-risk canon, why the pooled cv variant is the recommended default): docs/concepts/data-splitting.md.

CalibratedModel

CalibratedModel(model: Any, calibrator: BaseCalibrator, flow: str = 'prefit', cv: int = 5, ensemble: bool = False, random_state: int = 42)

Wrap any scoring model with a probcal calibrator (and optional offsets).

PARAMETER	DESCRIPTION
model	Duck-typed model with predict_proba(X) or decision_function(X). For flow="cv" it must also have fit(X, y) and be clonable. TYPE: object
calibrator	Unfitted calibrator instance (its parameters are cloned per fold in the cv flow via get_params). TYPE: BaseCalibrator
flow	"prefit": the model is already trained; fit(X_cal, y_cal) scores the calibration set and fits the calibrator — the canonical credit-risk flow. "cv": the model is cloned and retrained per fold; every observation is scored by a model that did not train on it. TYPE: (prefit, cv) DEFAULT: "prefit"
cv	Fold count for the cv flow (stratified, seeded). TYPE: int DEFAULT: 5
ensemble	False (recommended default): one calibrator on pooled out-of-fold scores, final model refit on all data — a single auditable mapping. True: keep the per-fold (model, calibrator) pairs and average their predictions. TYPE: bool DEFAULT: False
random_state	Seed for the fold assignment. TYPE: int DEFAULT: 42

ATTRIBUTE	DESCRIPTION
model_	The deployed model (the input model for prefit; the full-data refit for pooled cv). TYPE: object
calibrator_	The fitted calibrator (pooled/prefit flows). TYPE: BaseCalibrator
ensemble_	The fold pairs (ensemble flow only). TYPE: list[tuple[model, BaseCalibrator]]
offsets_	Appended offset stages, each separately inspectable. TYPE: list[LogitOffset]

Source code in src/probcal/wrapper.py

def __init__(
    self,
    model: Any,
    calibrator: BaseCalibrator,
    flow: str = "prefit",
    cv: int = 5,
    ensemble: bool = False,
    random_state: int = 42,
) -> None:
    self.model = model
    self.calibrator = calibrator
    self.flow = flow
    self.cv = cv
    self.ensemble = ensemble
    self.random_state = random_state

is_monotone_ property

is_monotone_: bool

Monotone iff the calibrator stage is (offsets always are).

affine_logit_coeffs_ property

affine_logit_coeffs_: tuple[float, float] | None

Composed (a, b + sum(deltas)) when the calibrator stage is affine.

fit

fit(X: object, y: object, sample_weight: object = None) -> Self

Fit the calibration stage per the configured flow.

Source code in src/probcal/wrapper.py

def fit(self, X: object, y: object, sample_weight: object = None) -> Self:
    """Fit the calibration stage per the configured flow."""
    if self.flow not in ("prefit", "cv"):
        raise ValueError(f"flow must be 'prefit' or 'cv', got {self.flow!r}")
    X_arr = np.asarray(X, dtype=np.float64)
    y_arr = validate_binary_y(y)
    w_arr = validate_weights(sample_weight, len(y_arr))
    self.offsets_: list[LogitOffset] = []
    self.ensemble_: list[tuple[Any, BaseCalibrator]] = []
    if self.flow == "prefit":
        self.model_ = self.model
        s = _model_scores(self.model_, X_arr)
        self.calibrator_ = self._fresh_calibrator().fit(s, y_arr, sample_weight=w_arr)
        self._cal_scores = s
    else:
        self._fit_cv(X_arr, y_arr, w_arr)
    self.fitted_ = True
    return self

predict_proba

predict_proba(X: object) -> ndarray

Calibrated (and offset) probabilities P(y=1) for new inputs.

Source code in src/probcal/wrapper.py

def predict_proba(self, X: object) -> np.ndarray:
    """Calibrated (and offset) probabilities ``P(y=1)`` for new inputs."""
    self._check_fitted()
    p = self._base_predict(np.asarray(X, dtype=np.float64))
    for off in self.offsets_:
        p = off.transform(p)
    return p

predict_proba_2d

predict_proba_2d(X: object) -> ndarray

Sklearn-style (n, 2) probability matrix.

Source code in src/probcal/wrapper.py

def predict_proba_2d(self, X: object) -> np.ndarray:
    """Sklearn-style ``(n, 2)`` probability matrix."""
    p = self.predict_proba(X)
    return np.column_stack([1.0 - p, p])

offset_to

offset_to(target_mean: float | None = None, delta: float | None = None, X: object = None) -> Self

Append an inspectable :class:LogitOffset stage.

Mode B (target_mean) anchors the portfolio mean of the current pipeline output — computed on X when given, else on the stored calibration scores (DECISIONS 48). The offset is never folded into the calibrator's parameters.

Source code in src/probcal/wrapper.py

def offset_to(
    self,
    target_mean: float | None = None,
    delta: float | None = None,
    X: object = None,
) -> Self:
    """Append an inspectable :class:`LogitOffset` stage.

    Mode B (``target_mean``) anchors the portfolio mean of the current
    pipeline output — computed on ``X`` when given, else on the stored
    calibration scores (DECISIONS 48). The offset is never folded into
    the calibrator's parameters.
    """
    self._check_fitted()
    if X is not None:
        p_now = self.predict_proba(X)
    else:
        p_now = self._base_predict_from_scores(self._cal_scores)
    off = LogitOffset(delta=delta, target_mean=target_mean)
    off.fit(p_now)
    self.offsets_.append(off)
    return self

interval_inverse

interval_inverse(lo: float, hi: float, *, space: str = 'probability', buffer_logit: float = 0.0) -> tuple[float, float]

Preimage of a calibrated interval through the full pipeline.

Composes right-to-left: the buffer shrinks the final interval, each offset subtracts its delta on the logit scale, and the calibrator's own inverse finishes the job. Returns bounds on the model's probability output (space="probability") or their logits.

Source code in src/probcal/wrapper.py

def interval_inverse(
    self,
    lo: float,
    hi: float,
    *,
    space: str = "probability",
    buffer_logit: float = 0.0,
) -> tuple[float, float]:
    """Preimage of a calibrated interval through the full pipeline.

    Composes right-to-left: the buffer shrinks the final interval, each
    offset subtracts its delta on the logit scale, and the calibrator's
    own inverse finishes the job. Returns bounds on the model's
    probability output (``space="probability"``) or their logits.
    """
    self._check_fitted()
    if self.ensemble_:
        raise NotImplementedError(
            "interval_inverse is not defined for the ensemble flow (K distinct maps); "
            "use ensemble=False for threshold translation"
        )
    if not 0.0 <= lo <= hi <= 1.0:
        raise ValueError(f"need 0 <= lo <= hi <= 1, got lo={lo}, hi={hi}")
    lo_b, hi_b = float(lo), float(hi)
    if buffer_logit > 0.0:
        if lo > 0.0:
            lo_b = float(expit(np.array([logit(np.array([lo]))[0] + buffer_logit]))[0])
        if hi < 1.0:
            hi_b = float(expit(np.array([logit(np.array([hi]))[0] - buffer_logit]))[0])
        if lo_b > hi_b:
            raise UnattainableTargetError(
                f"buffer_logit={buffer_logit} empties the calibrated interval [{lo}, {hi}]"
            )
    total_delta = sum(off.delta_ for off in self.offsets_)
    if total_delta != 0.0:
        if lo_b > 0.0:
            lo_b = float(expit(np.array([logit(np.array([lo_b]))[0] - total_delta]))[0])
        if hi_b < 1.0:
            hi_b = float(expit(np.array([logit(np.array([hi_b]))[0] - total_delta]))[0])
    return self.calibrator_.interval_inverse(lo_b, hi_b, space=space, buffer_logit=0.0)

interpret

interpret() -> Interpretation

Concatenated interpretation of the calibrator and every offset stage.

Source code in src/probcal/wrapper.py

def interpret(self) -> Interpretation:
    """Concatenated interpretation of the calibrator and every offset stage."""
    self._check_fitted()
    if self.ensemble_:
        parts = [cal.interpret() for _, cal in self.ensemble_]
    else:
        parts = [self.calibrator_.interpret()]
    parts += [off.interpret() for off in self.offsets_]
    names: tuple[str, ...] = ()
    values: tuple[float, ...] = ()
    messages: tuple[str, ...] = ()
    for part in parts:
        names += part.param_names
        values += part.param_values
        messages += part.messages
    return Interpretation(
        method=f"CalibratedModel[{', '.join(p.method for p in parts)}]",
        param_names=names,
        param_values=values,
        messages=messages,
    )

selection

CalibratorSelector: automatic method selection under nested validation.

The selector's scoring path only ever receives out-of-fold predictions — selection on fitting data is an unrepresentable state, not a documented misuse. Protocol, criteria, and report reading: docs/concepts/auto-selection.md.

CalibratorSelector

CalibratorSelector(candidates: dict[str, BaseCalibrator] | None = None, scoring: str = 'log_loss', cv: int = 5, random_state: int = 42)

Choose a calibrator by inner cross-validation on the calibration data.

PARAMETER	DESCRIPTION
candidates	Candidate instances (cloned per fold via get_params). None uses the spec's default menu: platt, temperature, beta_abm, isotonic, cir, histogram_mass, scaling_binning, ivap. The heavier methods (spline, BBQ, ENIR, CVAP) join by explicit opt-in. TYPE: dict[str, BaseCalibrator] or None DEFAULT: None
scoring	Out-of-fold selection criterion, lower is better. Plain ECE and Hosmer–Lemeshow are refused — see the metrics chapter's table. TYPE: (log_loss, brier, ici, smooth_ece, ece_sweep) DEFAULT: "log_loss"
cv	Inner stratified fold count. TYPE: int DEFAULT: 5
random_state	Seed for the fold assignment. TYPE: int DEFAULT: 42

ATTRIBUTE	DESCRIPTION
best_name_	Winning candidate's name. TYPE: str
best_calibrator_	The winner refitted on the full calibration set. TYPE: BaseCalibrator
report_	Ranked table: mean ± sd of the criterion, guardrail flags, chosen marker. TYPE: SelectionReport

Source code in src/probcal/selection.py

def __init__(
    self,
    candidates: dict[str, BaseCalibrator] | None = None,
    scoring: str = "log_loss",
    cv: int = 5,
    random_state: int = 42,
) -> None:
    self.candidates = candidates
    self.scoring = scoring
    self.cv = cv
    self.random_state = random_state

fit

fit(s: object, y: object, sample_weight: object = None) -> CalibratorSelector

Run the nested selection and refit the winner on all data.

Source code in src/probcal/selection.py

def fit(self, s: object, y: object, sample_weight: object = None) -> "CalibratorSelector":
    """Run the nested selection and refit the winner on all data."""
    if self.scoring not in _SCORERS:
        raise ValueError(
            f"scoring must be one of {sorted(_SCORERS)} (proper scores and accepted "
            f"binning-free alternatives), got {self.scoring!r}; plain ECE and "
            "Hosmer-Lemeshow are not selection criteria"
        )
    scorer = _SCORERS[self.scoring]
    menu = self.candidates if self.candidates is not None else _default_candidates()
    s_arr = validate_scores(s)
    y_arr = validate_binary_y(y)
    w_arr = validate_weights(sample_weight, len(s_arr))

    rng = np.random.default_rng(self.random_state)
    folds = np.empty(len(y_arr), dtype=np.int64)
    for cls in (0.0, 1.0):
        idx = np.flatnonzero(y_arr == cls)
        perm = rng.permutation(idx)
        folds[perm] = np.arange(len(perm)) % self.cv

    names = list(menu)
    means = np.empty(len(names))
    sds = np.empty(len(names))
    guards = np.empty(len(names), dtype=bool)
    for i, name in enumerate(names):
        proto = menu[name]
        fold_scores = np.empty(self.cv)
        oof = np.empty(len(y_arr))
        for k in range(self.cv):
            train, held = folds != k, folds == k
            cal = type(proto)(**proto.get_params())
            cal.fit(s_arr[train], y_arr[train], sample_weight=w_arr[train])
            pred_held = cal.predict_proba(s_arr[held])
            oof[held] = pred_held
            fold_scores[k] = scorer(y_arr[held], pred_held, sample_weight=w_arr[held])
        means[i] = fold_scores.mean()
        sds[i] = fold_scores.std(ddof=1)
        guards[i] = calibration_guardrails(y_arr, oof, sample_weight=w_arr).all_ok

    # Parsimony tie-break within one standard error of the best mean.
    best_idx = int(np.argmin(means))
    se_best = sds[best_idx] / np.sqrt(self.cv)
    tied = [i for i in range(len(names)) if means[i] <= means[best_idx] + se_best]
    winner = min(tied, key=lambda i: (_PARSIMONY.get(names[i], _PARSIMONY_UNKNOWN), means[i]))

    order = np.argsort(means, kind="stable")
    chosen = np.zeros(len(names), dtype=bool)
    chosen[winner] = True
    self.report_ = SelectionReport(
        methods=tuple(names[i] for i in order),
        score_mean=means[order],
        score_sd=sds[order],
        guardrails_ok=guards[order],
        chosen=chosen[order],
        criterion=self.scoring,
    )
    self.best_name_ = names[winner]
    proto = menu[self.best_name_]
    self.best_calibrator_ = type(proto)(**proto.get_params())
    self.best_calibrator_.fit(s_arr, y_arr, sample_weight=w_arr)
    return self

predict_proba

predict_proba(s: object) -> ndarray

Delegate to the refitted winner.

Source code in src/probcal/selection.py

def predict_proba(self, s: object) -> np.ndarray:
    """Delegate to the refitted winner."""
    return self.best_calibrator_.predict_proba(s)

interpret

interpret() -> Interpretation

Delegate to the refitted winner.

Source code in src/probcal/selection.py

def interpret(self) -> Interpretation:
    """Delegate to the refitted winner."""
    return self.best_calibrator_.interpret()

curves

Reliability-curve builders and the GiViTI-style calibration belt.

Numpy-only; every result is a frozen dataclass carrying both probability- and logit-scale coordinates, plotting-backend-agnostic (rendering lives in probcal.plots). Theory: docs/concepts/visualization.md.

References

Austin & Steyerberg (2014); Nattino, Finazzi & Bertolini (2014); Nattino, Lemeshow, Phillips, Finazzi & Bertolini (2017) — full records in the documentation. The belt is reimplemented from the papers; no GPL code is used.

EcceCurve dataclass

EcceCurve(frac: ndarray, cumdev: ndarray, sd_null: ndarray, stat_max: float, argmax_frac: float)

Cumulative-deviation walk over predictions sorted ascending (ECCE).

reliability_binned

reliability_binned(y: object, p: object, *, n_bins: int = 10, strategy: str = 'mass', sample_weight: object = None) -> ReliabilityCurve

Binned reliability curve with Wilson confidence intervals.

PARAMETER	DESCRIPTION
y	Outcomes and predicted probabilities. TYPE: array_like
p	Outcomes and predicted probabilities. TYPE: array_like
n_bins	Requested bin count. TYPE: int DEFAULT: 10
strategy	Equal-count (default) or equal-width bins. TYPE: (mass, width) DEFAULT: "mass"
sample_weight	Weights for the bin means; Wilson CIs use raw counts. TYPE: array_like or None DEFAULT: None

RETURNS	DESCRIPTION
ReliabilityCurve	Per-bin mean prediction, event rate, count, Wilson CI, and the logit-scale coordinates.

Source code in src/probcal/curves.py

def reliability_binned(
    y: object,
    p: object,
    *,
    n_bins: int = 10,
    strategy: str = "mass",
    sample_weight: object = None,
) -> ReliabilityCurve:
    """Binned reliability curve with Wilson confidence intervals.

    Parameters
    ----------
    y, p : array_like
        Outcomes and predicted probabilities.
    n_bins : int
        Requested bin count.
    strategy : {"mass", "width"}
        Equal-count (default) or equal-width bins.
    sample_weight : array_like or None
        Weights for the bin means; Wilson CIs use raw counts.

    Returns
    -------
    ReliabilityCurve
        Per-bin mean prediction, event rate, count, Wilson CI, and the
        logit-scale coordinates.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    idx, m = _bin_index(p_arr, n_bins, strategy)
    w_sum = np.bincount(idx, weights=w, minlength=m)
    keep = w_sum > 0
    wy = np.bincount(idx, weights=w * y_arr, minlength=m)[keep]
    wp = np.bincount(idx, weights=w * p_arr, minlength=m)[keep]
    counts = np.bincount(idx, minlength=m)[keep]
    w_kept = w_sum[keep]
    pred_mean = wp / w_kept
    event_rate = wy / w_kept
    ci_low, ci_high = _wilson(event_rate, counts.astype(np.float64))
    # The Wilson interval contains the point estimate analytically; enforce it
    # against floating-point noise at 0/1-rate bins (negative yerr otherwise).
    ci_low = np.minimum(np.clip(ci_low, 0.0, 1.0), event_rate)
    ci_high = np.maximum(np.clip(ci_high, 0.0, 1.0), event_rate)
    return ReliabilityCurve(
        pred_mean=pred_mean,
        event_rate=event_rate,
        count=counts.astype(np.int64),
        ci_low=ci_low,
        ci_high=ci_high,
        pred_mean_logit=logit(pred_mean),
    )

reliability_loess

reliability_loess(y: object, p: object, *, frac: float = 0.75, grid_size: int = 100, sample_weight: object = None) -> SmoothReliabilityCurve

LOESS-smoothed reliability curve on a grid (Austin & Steyerberg, 2014).

Source code in src/probcal/curves.py

def reliability_loess(
    y: object,
    p: object,
    *,
    frac: float = 0.75,
    grid_size: int = 100,
    sample_weight: object = None,
) -> SmoothReliabilityCurve:
    """LOESS-smoothed reliability curve on a grid (Austin & Steyerberg, 2014)."""
    y_arr, p_arr, _ = _prep(y, p, sample_weight)
    grid = _grid(p_arr, grid_size)
    rate = np.clip(loess(p_arr, y_arr, frac=frac, xeval=grid), 0.0, 1.0)
    return SmoothReliabilityCurve(grid_p=grid, grid_logit=logit(grid), event_rate=rate)

reliability_spline

reliability_spline(y: object, p: object, *, grid_size: int = 100, sample_weight: object = None) -> SmoothReliabilityCurve

Spline-smoothed reliability curve on a grid (penalized natural cubic spline of the outcome on the logit prediction).

Source code in src/probcal/curves.py

def reliability_spline(
    y: object,
    p: object,
    *,
    grid_size: int = 100,
    sample_weight: object = None,
) -> SmoothReliabilityCurve:
    """Spline-smoothed reliability curve on a grid (penalized natural cubic
    spline of the outcome on the logit prediction)."""
    from .spline import SplineCalibrator

    y_arr, p_arr, w = _prep(y, p, sample_weight)
    cal = SplineCalibrator()
    cal.fit(p_arr, y_arr, sample_weight=w)
    grid = _grid(p_arr, grid_size)
    return SmoothReliabilityCurve(
        grid_p=grid, grid_logit=logit(grid), event_rate=cal.predict_proba(grid)
    )

ecce_curve

ecce_curve(y: object, p: object, *, sample_weight: object = None) -> EcceCurve

Cumulative-deviation walk for the ECCE plot (Arrieta-Ibarra et al., 2022).

Sorts by prediction and accumulates weighted residuals, mirroring metrics.ecce exactly so stat_max agrees with the metric. sd_null is the pointwise standard deviation of the walk under calibration — an envelope for reading, not a simultaneous band.

Source code in src/probcal/curves.py

def ecce_curve(y: object, p: object, *, sample_weight: object = None) -> EcceCurve:
    """Cumulative-deviation walk for the ECCE plot (Arrieta-Ibarra et al., 2022).

    Sorts by prediction and accumulates weighted residuals, mirroring
    ``metrics.ecce`` exactly so ``stat_max`` agrees with the metric.
    ``sd_null`` is the pointwise standard deviation of the walk under
    calibration — an envelope for reading, not a simultaneous band.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    order = np.argsort(p_arr, kind="stable")
    n = len(p_arr)
    wsum = w.sum()
    cumdev = np.cumsum(w[order] * (y_arr[order] - p_arr[order])) / wsum
    # Pointwise H0 SD; reduces to sqrt(cumsum(p(1-p)))/n for unit weights.
    sd_null = np.sqrt(np.cumsum(w[order] ** 2 * p_arr[order] * (1.0 - p_arr[order]))) / wsum
    frac = np.arange(1, n + 1) / n
    k = int(np.argmax(np.abs(cumdev)))
    return EcceCurve(
        frac=frac,
        cumdev=cumdev,
        sd_null=sd_null,
        stat_max=float(np.abs(cumdev[k])),
        argmax_frac=float(frac[k]),
    )

calibration_belt

calibration_belt(y: object, p: object, *, confidence: tuple[float, float] = (0.8, 0.95), grid_size: int = 100, sample_weight: object = None) -> BeltResult

GiViTI-style calibration belt (Nattino et al., 2014, 2017).

Fits a polynomial logistic recalibration of the outcome on logit(p), selecting the degree by forward likelihood-ratio testing (p < 0.05 to add a term, capped at degree 4), then draws pointwise confidence bands from the information-matrix ellipsoid — a Wald approximation of the LR-region inversion (DECISIONS entry). The associated p-value tests the fitted polynomial against the identity. Where the band excludes the diagonal, the data reject calibration in that region.

Source code in src/probcal/curves.py

def calibration_belt(
    y: object,
    p: object,
    *,
    confidence: tuple[float, float] = (0.8, 0.95),
    grid_size: int = 100,
    sample_weight: object = None,
) -> BeltResult:
    """GiViTI-style calibration belt (Nattino et al., 2014, 2017).

    Fits a polynomial logistic recalibration of the outcome on
    ``logit(p)``, selecting the degree by forward likelihood-ratio testing
    (p < 0.05 to add a term, capped at degree 4), then draws pointwise
    confidence bands from the information-matrix ellipsoid — a Wald
    approximation of the LR-region inversion (DECISIONS entry). The
    associated p-value tests the fitted polynomial against the identity.
    Where the band excludes the diagonal, the data reject calibration in
    that region.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    z = logit(p_arr)

    def design(deg: int, t: np.ndarray) -> np.ndarray:
        return np.column_stack([t**k for k in range(deg + 1)])

    def loglik(beta: np.ndarray, deg: int) -> float:
        prob = np.clip(expit(design(deg, z) @ beta), 1e-12, 1.0 - 1e-12)
        return float(np.sum(w * (y_arr * np.log(prob) + (1.0 - y_arr) * np.log1p(-prob))))

    # Forward LR selection of the polynomial degree.
    degree = 1
    fit = irls_logistic(design(1, z), y_arr, w=w)
    ll = loglik(fit.beta, 1)
    while degree < 4:
        cand = irls_logistic(design(degree + 1, z), y_arr, w=w)
        ll_cand = loglik(cand.beta, degree + 1)
        lr = max(2.0 * (ll_cand - ll), 0.0)
        p_add = 1.0 - float(gammainc_lower(0.5, lr / 2.0))  # chi-square df=1
        if p_add >= 0.05:
            break
        degree += 1
        fit, ll = cand, ll_cand

    # Associated calibration test: fitted polynomial vs the identity map.
    ll_null = float(
        np.sum(
            w
            * (
                y_arr * np.log(np.clip(p_arr, 1e-12, 1))
                + (1.0 - y_arr) * np.log(np.clip(1.0 - p_arr, 1e-12, 1))
            )
        )
    )
    lr_cal = max(2.0 * (ll - ll_null), 0.0)
    df = degree + 1
    p_value = 1.0 - float(gammainc_lower(df / 2.0, lr_cal / 2.0))

    # Pointwise bands from the information-matrix ellipsoid.
    X = design(degree, z)
    mu = expit(np.clip(X @ fit.beta, -30.0, 30.0))
    info = (X * (w * mu * (1.0 - mu))[:, None]).T @ X
    info_inv = np.linalg.inv(info + 1e-10 * np.eye(df))
    grid_p = _grid(p_arr, grid_size)
    grid_z = logit(grid_p)
    Xg = design(degree, grid_z)
    eta = Xg @ fit.beta
    se_sq = np.einsum("ij,jk,ik->i", Xg, info_inv, Xg)
    bands = {}
    for conf in confidence:
        radius = np.sqrt(chi2_ppf(conf, float(df)) * se_sq)
        bands[conf] = (expit(eta - radius), expit(eta + radius))
    lo_80, hi_80 = bands[confidence[0]]
    lo_95, hi_95 = bands[confidence[1]]
    return BeltResult(
        grid_p=grid_p,
        grid_logit=grid_z,
        lower_80=lo_80,
        upper_80=hi_80,
        lower_95=lo_95,
        upper_95=hi_95,
        degree=degree,
        p_value=p_value,
    )

plots

Matplotlib plotting helpers (requires the [viz] extra; import-guarded).

All computation lives in probcal.curves and probcal.metrics; this module only renders. The logit-scale views are the flagship for low-PD portfolios: axis ticks sit at logit positions but are labeled in probabilities, so the low-probability region stays readable. Styling is applied per call via rc_context — global rcParams are never touched. Theory: docs/concepts/visualization.md.

plot_reliability

plot_reliability(curve: ReliabilityCurve, *, smooth: SmoothReliabilityCurve | None = None, scale: str = 'probability', y: object = None, p: object = None, annotate: bool = True, rug: bool = True, counts: bool = False, ax: Any = None) -> Any

Annotated reliability diagram: binned points with Wilson CIs, optional smooth overlay, stats box, and event/non-event rug.

scale="logit" stretches the low-probability region — the recommended view for PD portfolios. Bins whose event rate is exactly 0 or 1 have no finite logit and are omitted from the logit-scale point layer; they remain visible in the rug (or the counts=True margin).

Passing the raw y/p enables the stats box (annotate=True, computed by :func:probcal.metrics.reliability_summary) and the rug (rug=True, events along the top edge, non-events along the bottom, deterministically thinned to at most 1000 marks per class). Both are silently skipped when y/p are absent. counts=True restores the twin-axis count-bar margin.

Source code in src/probcal/plots.py

def plot_reliability(
    curve: ReliabilityCurve,
    *,
    smooth: SmoothReliabilityCurve | None = None,
    scale: str = "probability",
    y: object = None,
    p: object = None,
    annotate: bool = True,
    rug: bool = True,
    counts: bool = False,
    ax: Any = None,
) -> Any:
    """Annotated reliability diagram: binned points with Wilson CIs, optional
    smooth overlay, stats box, and event/non-event rug.

    ``scale="logit"`` stretches the low-probability region — the recommended
    view for PD portfolios. Bins whose event rate is exactly 0 or 1 have no
    finite logit and are omitted from the logit-scale point layer; they remain
    visible in the rug (or the ``counts=True`` margin).

    Passing the raw ``y``/``p`` enables the stats box (``annotate=True``,
    computed by :func:`probcal.metrics.reliability_summary`) and the rug
    (``rug=True``, events along the top edge, non-events along the bottom,
    deterministically thinned to at most 1000 marks per class). Both are
    silently skipped when ``y``/``p`` are absent. ``counts=True`` restores the
    twin-axis count-bar margin.
    """
    _require_mpl()
    if (y is None) != (p is None):
        raise ValueError("y and p must be given together")
    with _plt.rc_context(_STYLE):
        if ax is None:
            _, ax = _plt.subplots(figsize=(6.5, 6))
        if scale == "logit":
            keep = (curve.event_rate > 0.0) & (curve.event_rate < 1.0)
            x, ylo = curve.pred_mean_logit[keep], logit(curve.ci_low[keep])
            yv, yhi = logit(curve.event_rate[keep]), logit(curve.ci_high[keep])
            diag = np.linspace(x.min() - 0.5, x.max() + 0.5, 50)
            ax.plot(diag, diag, ls="--", c=_GREY, lw=1, label="identity")
            ax.errorbar(
                x,
                yv,
                yerr=[np.maximum(yv - ylo, 0.0), np.maximum(yhi - yv, 0.0)],
                fmt="o",
                ms=4,
                capsize=2,
                color=_BLUE,
                label="binned",
            )
            if smooth is not None:
                ax.plot(
                    smooth.grid_logit, logit(smooth.event_rate), lw=1.5, c=_ORANGE, label="smoothed"
                )
            _logit_axis(ax)
            ax.set_xlabel("predicted probability (logit scale)")
            ax.set_ylabel("event rate (logit scale)")
        else:
            ax.plot([0, 1], [0, 1], ls="--", c=_GREY, lw=1, label="identity")
            ax.errorbar(
                curve.pred_mean,
                curve.event_rate,
                yerr=[curve.event_rate - curve.ci_low, curve.ci_high - curve.event_rate],
                fmt="o",
                ms=4,
                capsize=2,
                color=_BLUE,
                label="binned",
            )
            if smooth is not None:
                ax.plot(smooth.grid_p, smooth.event_rate, lw=1.5, c=_ORANGE, label="smoothed")
            ax.set_xlabel("predicted probability")
            ax.set_ylabel("event rate")

        boxed = False
        if y is not None and p is not None:
            y_arr = np.asarray(y, dtype=np.float64)
            p_arr = np.asarray(p, dtype=np.float64)
            if rug:
                ev = _rug_subsample(p_arr[y_arr == 1.0])
                ne = _rug_subsample(p_arr[y_arr == 0.0])
                if scale == "logit":
                    ev, ne = logit(ev), logit(ne)
                tf = ax.get_xaxis_transform()
                ax.plot(
                    ev, np.full(len(ev), 0.99), transform=tf,
                    ls="none", marker="|", ms=7, c=_RED, alpha=0.25,
                )  # fmt: skip
                ax.plot(
                    ne, np.full(len(ne), 0.01), transform=tf,
                    ls="none", marker="|", ms=7, c="#777777", alpha=0.18,
                )  # fmt: skip
            if annotate:
                s = reliability_summary(y_arr, p_arr)
                txt = (
                    f"n = {s.n:,}\n"
                    f"events = {s.events:,}\n"
                    f"intercept = {s.intercept:+.3f}\n"
                    f"slope = {s.slope:.3f}\n"
                    f"ICI = {s.ici:.4f}\n"
                    f"E90 = {s.e90:.4f}\n"
                    f"Spiegelhalter p = {s.spiegelhalter_p:.3f}"
                )
                ax.text(0.03, 0.97, txt, transform=ax.transAxes, va="top", fontsize=9, bbox=_BOX)
                boxed = True
        if counts:
            # Count margin as a twin bar strip along the x-axis.
            ax2 = ax.twinx()
            xs = curve.pred_mean_logit if scale == "logit" else curve.pred_mean
            ax2.bar(xs, curve.count, width=np.ptp(xs) / (3 * len(xs) + 1), alpha=0.15, color=_GREY)
            ax2.set_yticks([])
        ax.legend(loc="lower right" if boxed else "upper left")
        return ax

plot_belt

plot_belt(belt: BeltResult, *, scale: str = 'probability', ax: Any = None) -> Any

GiViTI-style calibration belt with 80/95% bands and the test p-value.

Source code in src/probcal/plots.py

def plot_belt(belt: BeltResult, *, scale: str = "probability", ax: Any = None) -> Any:
    """GiViTI-style calibration belt with 80/95% bands and the test p-value."""
    _require_mpl()
    with _plt.rc_context(_STYLE):
        if ax is None:
            _, ax = _plt.subplots(figsize=(6.5, 6))
        if scale == "logit":
            x = belt.grid_logit
            ax.plot(x, x, ls="--", c=_GREY, lw=1)
            ax.fill_between(
                x, logit(belt.lower_95), logit(belt.upper_95), color=_BLUE, alpha=0.2, label="95%"
            )
            ax.fill_between(
                x, logit(belt.lower_80), logit(belt.upper_80), color=_BLUE, alpha=0.35, label="80%"
            )
            _logit_axis(ax)
        else:
            x = belt.grid_p
            ax.plot(x, x, ls="--", c=_GREY, lw=1)
            ax.fill_between(x, belt.lower_95, belt.upper_95, color=_BLUE, alpha=0.2, label="95%")
            ax.fill_between(x, belt.lower_80, belt.upper_80, color=_BLUE, alpha=0.35, label="80%")
        ax.set_title(f"calibration belt (degree {belt.degree}, p = {belt.p_value:.3g})")
        ax.set_xlabel("predicted probability")
        ax.set_ylabel("event rate")
        ax.legend(loc="upper left")
        return ax

plot_comparison

plot_comparison(before: ReliabilityCurve, after: ReliabilityCurve, *, scale: str = 'probability', labels: tuple[str, str] = ('before', 'after')) -> Any

Side-by-side reliability diagrams (pre/post calibration or offset).

Source code in src/probcal/plots.py

def plot_comparison(
    before: ReliabilityCurve,
    after: ReliabilityCurve,
    *,
    scale: str = "probability",
    labels: tuple[str, str] = ("before", "after"),
) -> Any:
    """Side-by-side reliability diagrams (pre/post calibration or offset)."""
    _require_mpl()
    with _plt.rc_context(_STYLE):
        fig, axes = _plt.subplots(1, 2, figsize=(12, 5.5), sharey=True)
        panel_colors = (_RED, _GREEN)
        for ax, curve, label, color in zip(
            axes, (before, after), labels, panel_colors, strict=True
        ):
            plot_reliability(curve, scale=scale, ax=ax)
            # Recolor the binned series to the panel's before/after semantics.
            for line in ax.lines:
                if line.get_label() == "binned":
                    line.set_color(color)
            for container in ax.containers:
                for artist in container.get_children():
                    artist.set_color(color)
            ax.set_title(label)
        return fig

plot_interval

plot_interval(intervals: ndarray, s: ndarray, *, ax: Any = None) -> Any

Venn–Abers interval widths against the score: where is calibration uncertain?

Source code in src/probcal/plots.py

def plot_interval(intervals: np.ndarray, s: np.ndarray, *, ax: Any = None) -> Any:
    """Venn–Abers interval widths against the score: where is calibration uncertain?"""
    _require_mpl()
    with _plt.rc_context(_STYLE):
        if ax is None:
            _, ax = _plt.subplots(figsize=(6.5, 4.5))
        p0, p1 = intervals[:, 0], intervals[:, 1]
        ax.fill_between(s, p0, p1, color=_BLUE, alpha=0.3, label="Venn–Abers interval")
        ax.plot(s, p1 / (1.0 - p0 + p1), lw=1.2, c=_ORANGE, label="scalarized")
        ax.set_xlabel("score")
        ax.set_ylabel("calibrated probability")
        ax.legend(loc="upper left")
        return ax

plot_selection

plot_selection(report: SelectionReport, *, ax: Any = None) -> Any

SelectionReport as a ranked dot plot with fold-spread whiskers.

Source code in src/probcal/plots.py

def plot_selection(report: SelectionReport, *, ax: Any = None) -> Any:
    """SelectionReport as a ranked dot plot with fold-spread whiskers."""
    _require_mpl()
    with _plt.rc_context(_STYLE):
        if ax is None:
            _, ax = _plt.subplots(figsize=(6.5, 0.6 * len(report.methods) + 1.5))
        order = np.argsort(report.score_mean)
        ys = np.arange(len(order))
        for rank, i in enumerate(order):
            ok = report.guardrails_ok[i]
            marker = "o" if ok else "x"
            color = _GREEN if report.chosen[i] else (_BLUE if ok else _RED)
            ax.errorbar(
                report.score_mean[i],
                rank,
                xerr=report.score_sd[i],
                fmt=marker,
                color=color,
                capsize=3,
            )
        ax.set_yticks(ys)
        ax.set_yticklabels([report.methods[i] for i in order])
        ax.set_xlabel(report.criterion)
        ax.set_title("calibrator selection (chosen in green; x = guardrail flag)")
        return ax

plot_ecce

plot_ecce(curves: Any, *, labels: Any = None, show_band: bool = True, ax: Any = None) -> Any

ECCE cumulative-drift walk(s) from :func:probcal.curves.ecce_curve.

Accepts a single EcceCurve or a sequence (e.g. raw vs calibrated). The grey envelope (show_band=True, from the first curve) is ±2 pointwise standard deviations under calibration — an aid for reading the walk, NOT a simultaneous confidence band; the formal max-statistic test of Arrieta-Ibarra et al. (2022) is out of scope for this release.

Source code in src/probcal/plots.py

def plot_ecce(
    curves: Any,
    *,
    labels: Any = None,
    show_band: bool = True,
    ax: Any = None,
) -> Any:
    """ECCE cumulative-drift walk(s) from :func:`probcal.curves.ecce_curve`.

    Accepts a single ``EcceCurve`` or a sequence (e.g. raw vs calibrated).
    The grey envelope (``show_band=True``, from the first curve) is ±2
    *pointwise* standard deviations under calibration — an aid for reading
    the walk, NOT a simultaneous confidence band; the formal max-statistic
    test of Arrieta-Ibarra et al. (2022) is out of scope for this release.
    """
    _require_mpl()
    if isinstance(curves, EcceCurve):
        curves = [curves]
    curves = list(curves)
    if labels is None:
        labels = [f"curve {i + 1}" for i in range(len(curves))]
    palette = [_RED, _GREEN, _BLUE, _ORANGE]
    with _plt.rc_context(_STYLE):
        if ax is None:
            _, ax = _plt.subplots(figsize=(7.5, 4.8))
        if show_band:
            c0 = curves[0]
            ax.fill_between(
                c0.frac,
                -2.0 * c0.sd_null,
                2.0 * c0.sd_null,
                color=_GREY,
                alpha=0.3,
                label="±2 SD under calibration (pointwise)",
            )
        ax.axhline(0.0, ls="--", c=_GREY, lw=1)
        for i, (c, label) in enumerate(zip(curves, labels, strict=True)):
            color = palette[i % len(palette)]
            ax.plot(
                c.frac, c.cumdev, lw=1.6, c=color, label=f"{label} (max drift {c.stat_max:.4f})"
            )
            ax.axvline(c.argmax_frac, ls=":", c=color, lw=1, alpha=0.7)
        ax.set_xlabel("cumulative share of portfolio (sorted by prediction)")
        ax.set_ylabel("cumulative deviation")
        ax.legend(loc="best")
        return ax

plot_grade_backtest

plot_grade_backtest(result: Any, *, log_scale: bool = True, ax: Any = None) -> Any

Per-grade traffic-light backtest chart (Jeffreys or exact binomial).

Observed default rates as circles colored by the grade's traffic light, grey 90% display intervals (ci_low/ci_high), and the assigned PDs as wide blue dashes. The intervals are display companions only — the verdict is carried by the lights from the unchanged one-sided tests, so no p-values are printed on the canvas. log_scale=True is the right default for PD grades spanning orders of magnitude.

Source code in src/probcal/plots.py

def plot_grade_backtest(result: Any, *, log_scale: bool = True, ax: Any = None) -> Any:
    """Per-grade traffic-light backtest chart (Jeffreys or exact binomial).

    Observed default rates as circles colored by the grade's traffic light,
    grey 90% display intervals (``ci_low``/``ci_high``), and the assigned PDs
    as wide blue dashes. The intervals are display companions only — the
    verdict is carried by the lights from the unchanged one-sided tests, so
    no p-values are printed on the canvas. ``log_scale=True`` is the right
    default for PD grades spanning orders of magnitude.
    """
    _require_mpl()
    light_color = {"green": _GREEN, "yellow": _AMBER, "amber": _AMBER, "red": _RED}
    name = "Jeffreys" if hasattr(result, "p_value") else "exact binomial"
    x = np.arange(len(result.grades))
    rate = result.k / result.n
    colors = [light_color.get(li, _GREY) for li in result.light]
    with _plt.rc_context(_STYLE):
        if ax is None:
            _, ax = _plt.subplots(figsize=(1.1 * len(x) + 3.0, 4.8))
        ax.scatter(x, result.pd, marker="_", s=500, c=_BLUE, zorder=2, label="assigned PD")
        ax.errorbar(
            x,
            rate,
            yerr=[np.maximum(rate - result.ci_low, 0.0), np.maximum(result.ci_high - rate, 0.0)],
            fmt="none",
            ecolor=_GREY,
            capsize=4,
            zorder=2,
        )
        ax.scatter(x, rate, s=90, c=colors, edgecolors="white", zorder=3, label="observed rate")
        for i in range(len(x)):
            ax.annotate(
                f"n={int(result.n[i]):,}\nk={int(result.k[i])}",
                xy=(float(x[i]), float(result.ci_high[i])),
                xytext=(0, 5),
                textcoords="offset points",
                ha="center",
                fontsize=8.5,
                color="#666666",
                clip_on=True,
            )
        if log_scale:
            ax.set_yscale("log")
        # Headroom so the n/k labels never collide with the title.
        lo, hi = ax.get_ylim()
        if log_scale:
            ax.set_ylim(lo, hi * (hi / lo) ** 0.12)
        else:
            ax.set_ylim(lo, hi + 0.12 * (hi - lo))
        ax.set_xticks(x)
        ax.set_xticklabels(result.grades)
        ax.set_xlabel("grade")
        ax.set_ylabel("default rate")
        ax.set_title(f"per-grade backtest ({name}, 90% display intervals)")
        ax.legend(loc="upper left")
        return ax

plot_offset_audit

plot_offset_audit(offset: Any, *, ax: Any = None) -> Any

Audit chart for a fitted :class:probcal.offset.LogitOffset stage.

Draws the offset map t -> t + delta on the logit scale against the identity, marks the pre- and post-adjustment central tendencies, and prints the audit numbers read directly from the fitted attributes. This chart audits the stage, not the outcomes — for the before/after guardrail comparison use LogitOffset.audit_report(y, p).

Source code in src/probcal/plots.py

def plot_offset_audit(offset: Any, *, ax: Any = None) -> Any:
    """Audit chart for a fitted :class:`probcal.offset.LogitOffset` stage.

    Draws the offset map ``t -> t + delta`` on the logit scale against the
    identity, marks the pre- and post-adjustment central tendencies, and
    prints the audit numbers read directly from the fitted attributes. This
    chart audits the *stage*, not the outcomes — for the before/after
    guardrail comparison use ``LogitOffset.audit_report(y, p)``.
    """
    _require_mpl()
    if not getattr(offset, "fitted_", False):
        raise RuntimeError("LogitOffset is not fitted; call fit() first")
    lo_g = float(logit(np.array([0.001]))[0])
    hi_g = float(logit(np.array([0.5]))[0])
    t = np.linspace(lo_g, hi_g, 200)
    lp = float(logit(np.array([offset.pre_mean_]))[0])
    lq = float(logit(np.array([offset.post_mean_]))[0])
    with _plt.rc_context(_STYLE):
        if ax is None:
            _, ax = _plt.subplots(figsize=(6.5, 6))
        ax.plot(t, t, ls="--", c=_GREY, lw=1, label="identity")
        ax.plot(t, t + offset.delta_, c=_BLUE, lw=1.6, label="offset map")
        if offset.target_mean is not None:
            ax.axhline(float(logit(np.array([offset.target_mean]))[0]), c=_GREY, lw=0.8, alpha=0.7)
        pre_xy = (lp, lp)
        post_xy = (lq - offset.delta_, lq)
        ax.scatter(*pre_xy, c=_RED, s=60, zorder=3, label="pre mean")
        ax.scatter(*post_xy, c=_GREEN, s=60, zorder=3, label="post mean")
        ax.annotate(
            "", xy=post_xy, xytext=pre_xy, arrowprops={"arrowstyle": "->", "color": "#555555"}
        )
        ax.annotate(
            f"δ = {offset.delta_:+.3f}",
            xy=((pre_xy[0] + post_xy[0]) / 2.0, (pre_xy[1] + post_xy[1]) / 2.0),
            xytext=(8, 0),
            textcoords="offset points",
            fontsize=9,
            color="#555555",
        )
        txt = (
            f"delta = {offset.delta_:+.4f} log-odds\n"
            f"odds factor = {math.exp(offset.delta_):.3f}\n"
            f"pre mean = {offset.pre_mean_:.4%}\n"
            f"post mean = {offset.post_mean_:.4%}\n"
            f"fitted {offset.timestamp_}"
        )
        ax.text(0.03, 0.97, txt, transform=ax.transAxes, va="top", fontsize=9, bbox=_BOX)
        _logit_axis(ax)
        ax.set_xlabel("input probability (logit scale)")
        ax.set_ylabel("shifted probability (logit scale)")
        ax.set_title("logit offset audit")
        ax.legend(loc="lower right")
        return ax

attribution

SHAP / additive-attribution adjustment to calibrated outputs.

Post-hoc calibration breaks SHAP local accuracy: base + sum(phi) reconstructs the raw score, not the calibrated probability. This module restores additivity on the calibrated scale — exactly for calibrators affine on the logit scale, and by the Aumann–Shapley secant rescaling in general. Theory, identifiability limits, and invariance properties: docs/concepts/shap-calibration.md.

No shap import: plain arrays are accepted, and shap.Explanation objects are duck-typed via their .values / .base_values attributes.

References

Lundberg & Lee (2017); Lundberg et al. (2020); Sundararajan, Taly & Yan (2017); Aumann & Shapley (1974) — full records in the documentation.

AdjustedAttribution dataclass

AdjustedAttribution(phi_adj: ndarray, base_adj: ndarray, target: ndarray, method_used: str, max_reconstruction_error: float)

Attributions rescaled to the calibrated output scale.

ATTRIBUTE	DESCRIPTION
phi_adj	Adjusted per-feature attributions. TYPE: numpy.ndarray of shape (n, d)
base_adj	Adjusted base values. TYPE: numpy.ndarray of shape (n,)
target	The calibrated output each row reconstructs (base_adj + phi_adj.sum(axis=1)). TYPE: numpy.ndarray of shape (n,)
method_used	"affine-exact" (exact Shapley values by linearity) or "aumann-shapley" (exact additivity; nonlinearity distributed proportionally to phi). TYPE: str
max_reconstruction_error	max |base_adj + sum(phi_adj) - target| over rows. TYPE: float

adjust_attributions

adjust_attributions(phi: object, base_value: object, calibrator: Any, *, scale: str = 'logit', method: str = 'auto') -> AdjustedAttribution

Rescale additive attributions so they sum to the calibrated output.

PARAMETER	DESCRIPTION
phi	Raw attributions on the model's score scale (log-odds margins for scale="logit", probabilities for scale="probability"). Objects exposing .values and .base_values are duck-typed; base_value is then ignored. TYPE: array_like of shape (n, d) or shap.Explanation-like
base_value	SHAP base value(s) on the same scale as phi. TYPE: float or array_like of shape (n,)
calibrator	Any object with predict_proba; affine_logit_coeffs_ (when not None) enables the exact affine path. TYPE: fitted calibrator
scale	Working scale of the attributions. Affine-exactness exists only on the logit scale. TYPE: (logit, probability) DEFAULT: "logit"
method	"auto" uses affine-exact when available, else Aumann–Shapley. TYPE: (auto, affine, aumann - shapley) DEFAULT: "auto"

RETURNS	DESCRIPTION
AdjustedAttribution

RAISES	DESCRIPTION
ValueError	If method="affine" is forced for a calibrator that is not affine on the logit scale (or on the probability scale, where no calibrator is affine).

Source code in src/probcal/attribution.py

def adjust_attributions(
    phi: object,
    base_value: object,
    calibrator: Any,
    *,
    scale: str = "logit",
    method: str = "auto",
) -> AdjustedAttribution:
    """Rescale additive attributions so they sum to the calibrated output.

    Parameters
    ----------
    phi : array_like of shape (n, d) or shap.Explanation-like
        Raw attributions on the model's score scale (log-odds margins for
        ``scale="logit"``, probabilities for ``scale="probability"``).
        Objects exposing ``.values`` and ``.base_values`` are duck-typed;
        ``base_value`` is then ignored.
    base_value : float or array_like of shape (n,)
        SHAP base value(s) on the same scale as ``phi``.
    calibrator : fitted calibrator
        Any object with ``predict_proba``; ``affine_logit_coeffs_`` (when not
        None) enables the exact affine path.
    scale : {"logit", "probability"}
        Working scale of the attributions. Affine-exactness exists only on
        the logit scale.
    method : {"auto", "affine", "aumann-shapley"}
        ``"auto"`` uses affine-exact when available, else Aumann–Shapley.

    Returns
    -------
    AdjustedAttribution

    Raises
    ------
    ValueError
        If ``method="affine"`` is forced for a calibrator that is not affine
        on the logit scale (or on the probability scale, where no calibrator
        is affine).
    """
    if scale not in ("logit", "probability"):
        raise ValueError(f"scale must be 'logit' or 'probability', got {scale!r}")
    if method not in ("auto", "affine", "aumann-shapley"):
        raise ValueError(f"unknown method {method!r}")
    phi_arr, base_arr = _extract(phi, base_value)
    s = base_arr + phi_arr.sum(axis=1)

    coeffs = getattr(calibrator, "affine_logit_coeffs_", None)
    affine_available = scale == "logit" and coeffs is not None
    if method == "affine" and not affine_available:
        raise ValueError(
            "method='affine' requires a calibrator affine on the logit scale "
            "(affine_logit_coeffs_ is None, or scale='probability' was requested)"
        )
    use_affine = affine_available and method in ("auto", "affine")

    if scale == "logit":

        def g_work(t: np.ndarray) -> np.ndarray:
            return logit(calibrator.predict_proba(expit(t)))

    else:

        def g_work(t: np.ndarray) -> np.ndarray:
            return calibrator.predict_proba(np.clip(t, 1e-12, 1.0 - 1e-12))

    target = g_work(s)

    if use_affine:
        assert coeffs is not None
        a, b = coeffs
        phi_adj = a * phi_arr
        base_adj = a * base_arr + b
        method_used = "affine-exact"
    else:
        g_s0 = g_work(base_arr)
        diff = s - base_arr
        multiplier = np.empty(len(s))
        regular = np.abs(diff) >= _DEGENERATE_EPS
        multiplier[regular] = (target[regular] - g_s0[regular]) / diff[regular]
        if np.any(~regular):
            b0 = base_arr[~regular]
            h = _CENTRAL_DIFF_H
            multiplier[~regular] = (g_work(b0 + h) - g_work(b0 - h)) / (2.0 * h)
        phi_adj = phi_arr * multiplier[:, None]
        base_adj = target - phi_adj.sum(axis=1)
        # base_adj equals g(s0) exactly on regular rows (telescoping); the
        # assignment above additionally zeroes reconstruction error on
        # degenerate rows where the local-slope multiplier is approximate.
        method_used = "aumann-shapley"

    recon_err = float(np.max(np.abs(base_adj + phi_adj.sum(axis=1) - target), initial=0.0))
    return AdjustedAttribution(
        phi_adj=phi_adj,
        base_adj=base_adj,
        target=target,
        method_used=method_used,
        max_reconstruction_error=recon_err,
    )

thresholds

Calibrated-to-raw interval and masterscale-band mapping.

Thin functional wrappers over the calibrators' interval_inverse protocol (spec §10): numpy-only, arrays and floats, no knowledge of any consumer. The canonical rating-grade workflow — masterscale bands defined on calibrated PD, translated once per recalibration into raw-score intervals — is calibrated_bands_to_raw; its output plugs directly into band-style raw targets of a counterfactual engine.

calibrated_interval_to_raw

calibrated_interval_to_raw(calibrator: object, lo: float, hi: float, *, space: str = 'probability', buffer_logit: float = 0.0) -> tuple[float, float]

Translate one calibrated-probability interval into raw-score bounds.

PARAMETER	DESCRIPTION
calibrator	Any object implementing the duck-typed protocol interval_inverse(lo, hi, *, space, buffer_logit) with is_monotone_. TYPE: fitted calibrator
lo	Calibrated bounds; lo=0 / hi=1 map to the full raw range. TYPE: float
hi	Calibrated bounds; lo=0 / hi=1 map to the full raw range. TYPE: float
space	Scale of the returned bounds. TYPE: (probability, logit) DEFAULT: "probability"
buffer_logit	Robustness margin applied in logit space before inversion. TYPE: float DEFAULT: 0.0

Source code in src/probcal/thresholds.py

def calibrated_interval_to_raw(
    calibrator: object,
    lo: float,
    hi: float,
    *,
    space: str = "probability",
    buffer_logit: float = 0.0,
) -> tuple[float, float]:
    """Translate one calibrated-probability interval into raw-score bounds.

    Parameters
    ----------
    calibrator : fitted calibrator
        Any object implementing the duck-typed protocol
        ``interval_inverse(lo, hi, *, space, buffer_logit)`` with
        ``is_monotone_``.
    lo, hi : float
        Calibrated bounds; ``lo=0`` / ``hi=1`` map to the full raw range.
    space : {"probability", "logit"}
        Scale of the returned bounds.
    buffer_logit : float
        Robustness margin applied in logit space before inversion.
    """
    return calibrator.interval_inverse(lo, hi, space=space, buffer_logit=buffer_logit)  # type: ignore[attr-defined]

calibrated_bands_to_raw

calibrated_bands_to_raw(calibrator: object, bands: dict, *, space: str = 'probability', buffer_logit: float = 0.0) -> dict

Translate a masterscale {grade: (lo, hi)} on calibrated PD to raw intervals.

Grade edges are policy artifacts that outlive model versions; this translation is what changes when the calibrator is refitted.

Source code in src/probcal/thresholds.py

def calibrated_bands_to_raw(
    calibrator: object,
    bands: dict,
    *,
    space: str = "probability",
    buffer_logit: float = 0.0,
) -> dict:
    """Translate a masterscale ``{grade: (lo, hi)}`` on calibrated PD to raw intervals.

    Grade edges are policy artifacts that outlive model versions; this
    translation is what changes when the calibrator is refitted.
    """
    return {
        grade: calibrated_interval_to_raw(
            calibrator, lo, hi, space=space, buffer_logit=buffer_logit
        )
        for grade, (lo, hi) in bands.items()
    }

datasets

Synthetic dataset generators (make_pd_portfolio).

PdPortfolio dataclass

PdPortfolio(scores: ndarray, y: ndarray, p_true: ndarray)

Synthetic PD portfolio: model scores, outcomes, and the true probabilities.

ATTRIBUTE	DESCRIPTION
scores	The model's reported PDs — miscalibrated unless generated with slope=1, asymmetry=0, intercept=0. TYPE: ndarray
y	Bernoulli outcomes drawn from p_true. TYPE: ndarray
p_true	True conditional probabilities (mean anchored at event_rate). TYPE: ndarray

make_pd_portfolio

make_pd_portfolio(n: int = 5000, *, event_rate: float = 0.03, slope: float = 0.7, intercept: float = 0.0, asymmetry: float = 0.4, score_location: float = -3.2, score_scale: float = 1.1, random_state: int = 42) -> PdPortfolio

Generate a synthetic, controllably miscalibrated PD portfolio.

The model's scores are drawn as s = sigma(N(score_location, score_scale)); the true probability follows the beta-calibration family

logit p_true = a_lo * ln(s) - a_hi * ln(1 - s) + c

with a_lo = slope * (1 + asymmetry) (low-PD tail) and a_hi = slope (high tail), so asymmetry != 0 produces exactly the one-sided tail distortion low-event-rate portfolios exhibit, and BetaCalibrator can recover the generative exponents. c absorbs intercept plus a portfolio-level anchor solved so that mean(p_true) == event_rate (unique by monotonicity, via bisection). With slope=1, asymmetry=0, intercept=0 the scores are exactly calibrated.

PARAMETER	DESCRIPTION
n	Portfolio size. TYPE: int DEFAULT: 5000
event_rate	Target mean of p_true (the central tendency), ~3% by default. TYPE: float DEFAULT: 0.03
slope	Base exponent of the distortion; < 1 means the model's scores are too spread out (overconfident). TYPE: float DEFAULT: 0.7
intercept	Additional log-odds shift applied before the mean anchor is solved. TYPE: float DEFAULT: 0.0
asymmetry	Relative extra distortion of the low-PD tail (a_lo/a_hi - 1). TYPE: float DEFAULT: 0.4
score_location	Parameters of the normal generating the score logits. TYPE: float DEFAULT: -3.2
score_scale	Parameters of the normal generating the score logits. TYPE: float DEFAULT: -3.2
random_state	Seed. TYPE: int DEFAULT: 42

RETURNS	DESCRIPTION
PdPortfolio

Source code in src/probcal/datasets.py

def make_pd_portfolio(
    n: int = 5000,
    *,
    event_rate: float = 0.03,
    slope: float = 0.7,
    intercept: float = 0.0,
    asymmetry: float = 0.4,
    score_location: float = -3.2,
    score_scale: float = 1.1,
    random_state: int = 42,
) -> PdPortfolio:
    """Generate a synthetic, controllably miscalibrated PD portfolio.

    The model's scores are drawn as ``s = sigma(N(score_location,
    score_scale))``; the true probability follows the beta-calibration family

    ``logit p_true = a_lo * ln(s) - a_hi * ln(1 - s) + c``

    with ``a_lo = slope * (1 + asymmetry)`` (low-PD tail) and ``a_hi = slope``
    (high tail), so ``asymmetry != 0`` produces exactly the one-sided tail
    distortion low-event-rate portfolios exhibit, and `BetaCalibrator` can
    recover the generative exponents. ``c`` absorbs ``intercept`` plus a
    portfolio-level anchor solved so that ``mean(p_true) == event_rate``
    (unique by monotonicity, via bisection). With ``slope=1, asymmetry=0,
    intercept=0`` the scores are exactly calibrated.

    Parameters
    ----------
    n : int
        Portfolio size.
    event_rate : float
        Target mean of ``p_true`` (the central tendency), ~3% by default.
    slope : float
        Base exponent of the distortion; ``< 1`` means the model's scores are
        too spread out (overconfident).
    intercept : float
        Additional log-odds shift applied before the mean anchor is solved.
    asymmetry : float
        Relative extra distortion of the low-PD tail (``a_lo/a_hi - 1``).
    score_location, score_scale : float
        Parameters of the normal generating the score logits.
    random_state : int
        Seed.

    Returns
    -------
    PdPortfolio
    """
    if not 0.0 < event_rate < 1.0:
        raise ValueError("event_rate must lie in (0, 1)")
    rng = np.random.default_rng(random_state)
    s = expit(rng.normal(score_location, score_scale, n))
    a_lo = slope * (1.0 + asymmetry)
    a_hi = slope
    core = a_lo * np.log(s) - a_hi * np.log1p(-s) + intercept

    identity_case = slope == 1.0 and asymmetry == 0.0 and intercept == 0.0
    if identity_case:
        p_true = s.copy()
    else:

        def gap(c: float) -> float:
            return float(np.mean(expit(core + c))) - event_rate

        c_anchor = bisect(gap, -60.0, 60.0, tol=1e-14)
        p_true = expit(core + c_anchor)
    y = (rng.random(n) < p_true).astype(np.float64)
    return PdPortfolio(scores=s, y=y, p_true=p_true)