API: metrics and tests¶

metrics ¶

Calibration metrics and statistical tests (flat re-exports).

Implementation lives in the submodules (scores, binned, smooth, regression, kernel, grade) and in private helpers (_evaluate, _summary, _binstats, _grades, _common). The sample-weight convention shared by every metric is documented in :mod:probcal.metrics._common. Selection guidance — what may be optimized and what is report-only — is the table in docs/concepts/metrics.md.

GroupedMetricReport `dataclass` ¶

GroupedMetricReport(pooled: MetricReport, groups: tuple[str, ...], reports: tuple[MetricReport, ...], counts: ndarray)

Bases: _ResultBase

Per-group metric reports plus a pooled report, from metrics.evaluate(by=...).

ATTRIBUTE	DESCRIPTION
`pooled`	Report computed on the full, ungrouped data (the `by=None` report, using `seed` unchanged). TYPE: `MetricReport`
`groups`	Sorted, stringified group labels. TYPE: `tuple of str`
`reports`	Per-group reports, aligned with `groups`. Group `i` (in this sorted order) is computed with `seed + 1000 * i`, so results are reproducible independent of the label values themselves. TYPE: `tuple of MetricReport`
`counts`	Observation count per group, aligned with `groups`. TYPE: `ndarray`

to_frame ¶

to_frame() -> object

Rows as a list of dicts, or a pandas DataFrame when pandas is importable.

Each row is {"group", "metric", "value", "ci_low", "ci_high"}; the pooled report is included under the group label "pooled", which is therefore reserved — a group of your own named "pooled" is indistinguishable from it in this frame.

Source code in src/probcal/_results.py

def to_frame(self) -> object:
    """Rows as a list of dicts, or a pandas DataFrame when pandas is importable.

    Each row is ``{"group", "metric", "value", "ci_low", "ci_high"}``;
    the pooled report is included under the group label ``"pooled"``,
    which is therefore reserved — a group of your own named "pooled"
    is indistinguishable from it in this frame.
    """
    return _to_frame(
        [
            {"group": group, "metric": n, "value": v, "ci_low": lo, "ci_high": hi}
            for group, n, v, lo, hi in self._rows()
        ]
    )

ReliabilitySummary `dataclass` ¶

ReliabilitySummary(n: int, events: int, intercept: float, slope: float, ici: float, e90: float, spiegelhalter_p: float)

Stats-box aggregate for the annotated reliability diagram.

ATTRIBUTE	DESCRIPTION
`n`	Observation count. TYPE: `int`
`events`	Event count (`sum(y)`). TYPE: `int`
`intercept`	Calibration-in-the-large intercept (log-odds). TYPE: `float`
`slope`	Cox calibration slope. TYPE: `float`
`ici`	Integrated calibration index. TYPE: `float`
`e90`	90th percentile of the LOESS distances. TYPE: `float`
`spiegelhalter_p`	Spiegelhalter test p-value. TYPE: `float`

HosmerLemeshowResult `dataclass` ¶

HosmerLemeshowResult(statistic: float, df: int, p_value: float)

Hosmer–Lemeshow chi-square test (report-only; never a selection criterion).

ATTRIBUTE	DESCRIPTION
`statistic`	Chi-square test statistic. TYPE: `float`
`df`	Degrees of freedom: non-empty groups minus `ddof` (`0` when no degree of freedom is left). TYPE: `int`
`p_value`	Upper-tail p-value of the statistic; `nan` when `df == 0`. TYPE: `float`

BinomialGradeResult `dataclass` ¶

BinomialGradeResult(grades: tuple, n: ndarray, k: ndarray, pd: ndarray, p_exact: ndarray, p_normal: ndarray, light: tuple, ci_low: ndarray, ci_high: ndarray)

Bases: _ResultBase

Exact and approximate one-sided binomial backtest per rating grade.

ATTRIBUTE	DESCRIPTION
`grades`	Grade labels: the `order` passed, else best to worst for a `Masterscale`, else sorted. TYPE: `tuple of str`
`n`	Observation count per grade. TYPE: `ndarray`
`k`	Default count per grade. TYPE: `ndarray`
`pd`	Assigned PD per grade (mean of `p` within the grade). TYPE: `ndarray`
`p_exact`	Exact binomial tail p-value per grade (also available as :attr:`p_value`, the name the other tests use). TYPE: `ndarray`
`p_normal`	Normal-approximation p-value per grade. TYPE: `ndarray`
`light`	Traffic light per grade (`"green"`, `"amber"`, or `"red"`), derived from `p_exact`. TYPE: `tuple of str`
`ci_low, ci_high`	90% Clopper-Pearson display interval for the observed rate. TYPE: `ndarray`

p_value `property` ¶

p_value: ndarray

Alias of :attr:p_exact (the exact test decides the traffic light).

HlEResult `dataclass` ¶

HlEResult(e_value: float, p_value: float, grades: tuple[str, ...], e_grade: ndarray, construction: str)

Bases: _ResultBase

Mixture-LR grade e-test result (:func:hl_e_test).

ATTRIBUTE	DESCRIPTION
`e_value`	The test e-value, `exp(log E)` with `log E = sum(log E_g)` over grades -- a product of per-grade e-values, itself a valid e-value for the joint null. TYPE: `float`
`p_value`	`min(1, 1 / e_value)` -- a valid (if generally conservative) p-value derived from the e-value via Markov's inequality. TYPE: `float`
`grades`	Grade labels: best to worst for a `Masterscale`, else sorted. TYPE: `tuple of str`
`e_grade`	Per-grade e-value, aligned with `grades`; `e_value` is their product. TYPE: `ndarray`
`construction`	Always `"mixture-lr"` -- recorded so a serialized or logged result names its own construction unambiguously. TYPE: `str`

Examples:

>>> import numpy as np
>>> from probcal.metrics import hl_e_test
>>> rng = np.random.default_rng(0)
>>> p = np.full(200, 0.1)
>>> y = (rng.random(200) < 0.1).astype(float)
>>> grades = np.array(["A"] * 100 + ["B"] * 100)
>>> res = hl_e_test(y, p, grades)
>>> res.construction
'mixture-lr'
>>> bool(np.isclose(res.e_value, np.prod(res.e_grade), rtol=1e-9))
True
>>> res.p_value == min(1.0, 1.0 / res.e_value)
True

interpret ¶

interpret() -> Interpretation

Read one e-value sentence per grade and the test-level conclusion.

RETURNS	DESCRIPTION
`Interpretation`	`method="HlETest"`, parameters `e_value`/`p_value`, one per-grade message plus a closing sentence on how to read the e-value/p-value pair.

Examples:

>>> import numpy as np
>>> from probcal.metrics import hl_e_test
>>> rng = np.random.default_rng(0)
>>> p = np.full(200, 0.1)
>>> y = (rng.random(200) < 0.1).astype(float)
>>> grades = np.array(["A"] * 100 + ["B"] * 100)
>>> interp = hl_e_test(y, p, grades).interpret()
>>> interp.method
'HlETest'
>>> "grade A: e =" in interp.messages[0]
True

Source code in src/probcal/metrics/_safe.py

def interpret(self) -> Interpretation:
    """Read one e-value sentence per grade and the test-level conclusion.

    Returns
    -------
    Interpretation
        ``method="HlETest"``, parameters ``e_value``/``p_value``, one
        per-grade message plus a closing sentence on how to read the
        e-value/p-value pair.

    Examples
    --------
    >>> import numpy as np
    >>> from probcal.metrics import hl_e_test
    >>> rng = np.random.default_rng(0)
    >>> p = np.full(200, 0.1)
    >>> y = (rng.random(200) < 0.1).astype(float)
    >>> grades = np.array(["A"] * 100 + ["B"] * 100)
    >>> interp = hl_e_test(y, p, grades).interpret()
    >>> interp.method
    'HlETest'
    >>> "grade A: e =" in interp.messages[0]
    True
    """
    grade_messages = tuple(
        f"grade {g}: e = {e:.4g}" for g, e in zip(self.grades, self.e_grade, strict=True)
    )
    closing = (
        f"e = {self.e_value:.4g}; e >= 1/alpha rejects H0 (miscalibration) at level alpha "
        f"(Ville/Markov, single fixed-sample look); p = min(1, 1/e) = {self.p_value:.4g} "
        "is a valid p-value derived from the same evidence."
    )
    return Interpretation(
        method="HlETest",
        param_names=("e_value", "p_value"),
        param_values=(self.e_value, self.p_value),
        messages=grade_messages + (closing,),
    )

JeffreysGradeResult `dataclass` ¶

JeffreysGradeResult(grades: tuple, n: ndarray, k: ndarray, pd: ndarray, p_value: ndarray, light: tuple, ci_low: ndarray, ci_high: ndarray)

Bases: _ResultBase

Jeffreys-posterior backtest per rating grade (ECB IRB practice).

ATTRIBUTE	DESCRIPTION
`grades`	Grade labels: the `order` passed, else best to worst for a `Masterscale`, else sorted. TYPE: `tuple of str`
`n`	Observation count per grade. TYPE: `ndarray`
`k`	Default count per grade. TYPE: `ndarray`
`pd`	Assigned PD per grade (mean of `p` within the grade). TYPE: `ndarray`
`p_value`	Posterior `P(theta <= PD \| k, n)` per grade. TYPE: `ndarray`
`light`	Traffic light per grade (`"green"`, `"amber"`, or `"red"`), derived from `p_value`. TYPE: `tuple of str`
`ci_low, ci_high`	Central 90% Jeffreys posterior display interval. TYPE: `ndarray`

PlutoTascheResult `dataclass` ¶

PlutoTascheResult(grades: tuple[str, ...], n: ndarray, d: ndarray, n_pooled: ndarray, d_pooled: ndarray, pd_upper: ndarray, confidence: float, monotonized: bool)

Bases: _ResultBase

Pluto-Tasche one-period most-prudent PD per rating grade.

ATTRIBUTE	DESCRIPTION
`grades`	Grade labels, best to worst, in the order given to :func:`pluto_tasche` / :func:`pluto_tasche_from_arrays`. TYPE: `tuple of str`
`n`	Own obligor count per grade (weighted sum if fitted from arrays with `sample_weight`). TYPE: `ndarray`
`d`	Own default count per grade (weighted sum likewise). TYPE: `ndarray`
`n_pooled`	Obligor count pooled with all worse grades: `n_pooled[i] = sum(n[i:])`. TYPE: `ndarray`
`d_pooled`	Default count pooled the same way. TYPE: `ndarray`
`pd_upper`	Most-prudent PD per grade: the one-sided Clopper-Pearson upper bound of the pooled default rate at `confidence`. TYPE: `ndarray`
`confidence`	Confidence level used for every grade's bound. TYPE: `float`
`monotonized`	`True` if the raw bounds were not already non-decreasing best to worst and `pd_upper` is their running maximum (a prudent hull that only ever raises a bound). Pooled sets are nested, so this happens only when a noisy grade's own rate exceeds the worse-grade pool it joins. TYPE: `bool`

interpret ¶

interpret() -> Interpretation

Read one audit sentence per grade: own counts, pooling, and the bound.

RETURNS	DESCRIPTION
`Interpretation`	`method="PlutoTasche"`, one `pd_upper.<grade>` parameter and one audit sentence per grade.

Examples:

>>> import numpy as np
>>> from probcal.metrics import pluto_tasche
>>> res = pluto_tasche(
...     np.array([100.0, 400.0, 300.0]),
...     np.array([0.0, 0.0, 0.0]),
...     confidence=0.9,
...     grades=("A", "B", "C"),
... )
>>> msg = res.interpret().messages[0]
>>> "grade A: 0 defaults among 100 obligors" in msg
True
>>> "most-prudent PD at 90% confidence = 0.29%" in msg
True

Source code in src/probcal/metrics/_conservative.py

def interpret(self) -> Interpretation:
    """Read one audit sentence per grade: own counts, pooling, and the bound.

    Returns
    -------
    Interpretation
        ``method="PlutoTasche"``, one ``pd_upper.<grade>`` parameter and
        one audit sentence per grade.

    Examples
    --------
    >>> import numpy as np
    >>> from probcal.metrics import pluto_tasche
    >>> res = pluto_tasche(
    ...     np.array([100.0, 400.0, 300.0]),
    ...     np.array([0.0, 0.0, 0.0]),
    ...     confidence=0.9,
    ...     grades=("A", "B", "C"),
    ... )
    >>> msg = res.interpret().messages[0]
    >>> "grade A: 0 defaults among 100 obligors" in msg
    True
    >>> "most-prudent PD at 90% confidence = 0.29%" in msg
    True
    """
    param_names = tuple(f"pd_upper.{g}" for g in self.grades)
    param_values = tuple(float(v) for v in self.pd_upper)
    messages = tuple(
        f"grade {g}: {self.d[i]:g} defaults among {self.n[i]:g} obligors; "
        f"pooled with worse grades (n*={self.n_pooled[i]:g}, d*={self.d_pooled[i]:g}); "
        f"most-prudent PD at {self.confidence:.0%} confidence = {self.pd_upper[i]:.2%}"
        for i, g in enumerate(self.grades)
    )
    return Interpretation(
        method="PlutoTasche",
        param_names=param_names,
        param_values=param_values,
        messages=messages,
    )

SkceTestResult `dataclass` ¶

SkceTestResult(statistic: float, estimator: str, method: str, p_value: float, p_value_bound: float, bandwidth: float, n_boot: int | None)

One-sided SKCE calibration test (H0: calibrated; large positive rejects).

ATTRIBUTE	DESCRIPTION
`statistic`	SKCE point estimate (`"ul"` for the asymptotic method, `"uq"` for the bootstrap method). TYPE: `float`
`estimator`	Estimator used for `statistic` (`"ul"` or `"uq"`). TYPE: `str`
`method`	Test method used (`"asymptotic"` or `"bootstrap"`). TYPE: `str`
`p_value`	Test p-value. TYPE: `float`
`p_value_bound`	Distribution-free worst-case p-value bound (valid without asymptotics). TYPE: `float`
`bandwidth`	Kernel bandwidth used (resolved from `bandwidth=None` if applicable). TYPE: `float`
`n_boot`	Bootstrap replicate count for the bootstrap method; `None` for the asymptotic method. TYPE: `int or None`

CalibrationTestResult `dataclass` ¶

CalibrationTestResult(statistic: float, p_value: float, alpha: float, beta: float, df: float = 2.0)

Likelihood-ratio test of (intercept, slope) = (0, 1) — the Cox-framed weak calibration test.

ATTRIBUTE	DESCRIPTION
`statistic`	Likelihood-ratio test statistic, referred to `chi2(df)` (the Rao-Scott corrected statistic for non-uniform weights). TYPE: `float`
`p_value`	Upper-tail p-value of the statistic. TYPE: `float`
`alpha`	Fitted intercept (also available as :attr:`intercept`). TYPE: `float`
`beta`	Fitted slope (also available as :attr:`slope`); `nan` when `logit(p)` is constant. TYPE: `float`
`df`	Degrees of freedom of the reference chi-square: `2` (`1` when only the intercept is testable), smaller and fractional after the Rao-Scott correction for non-uniform weights. TYPE: `float`

intercept `property` ¶

intercept: float

Fitted intercept; the name used by the rest of the package for alpha.

slope `property` ¶

slope: float

Fitted slope; the name used by the rest of the package for beta.

GuardrailReport `dataclass` ¶

GuardrailReport(slope: float, intercept: float, spiegelhalter_p: float, slope_ok: bool, intercept_ok: bool, spiegelhalter_ok: bool, all_ok: bool)

Three-flag calibration health summary used across the package.

Thresholds are conventions, not theorems: slope within [0.9, 1.1], intercept within +/-0.1 log-odds, Spiegelhalter p above 0.05.

ATTRIBUTE	DESCRIPTION
`slope`	Fitted Cox calibration slope. TYPE: `float`
`intercept`	Fitted calibration-in-the-large intercept (log-odds). TYPE: `float`
`spiegelhalter_p`	Spiegelhalter test p-value. TYPE: `float`
`slope_ok`	Whether `slope` lies in `[0.9, 1.1]`. TYPE: `bool`
`intercept_ok`	Whether `abs(intercept) <= 0.1`. TYPE: `bool`
`spiegelhalter_ok`	Whether `spiegelhalter_p > 0.05`. TYPE: `bool`
`all_ok`	Conjunction of the three flags above. TYPE: `bool`

LogLossDecomposition `dataclass` ¶

LogLossDecomposition(calibration: float, refinement: float)

Calibration/refinement split of the log loss via a plug-in recalibration curve (LOESS).

ATTRIBUTE	DESCRIPTION
`calibration`	Mean KL divergence between the plug-in and predicted Bernoullis. TYPE: `float`
`refinement`	Mean entropy of the plug-in Bernoulli. TYPE: `float`

MurphyCurve `dataclass` ¶

MurphyCurve(thresholds: ndarray, score: ndarray, n: int)

Bases: _ResultBase

Murphy diagram: mean elementary score of the binary mean functional across thresholds.

ATTRIBUTE	DESCRIPTION
`thresholds`	Threshold grid :math:`\theta \in [0, 1]`. TYPE: `ndarray`
`score`	Weighted mean elementary score :math:`S_\theta` at each threshold. TYPE: `ndarray`
`n`	Number of observations. TYPE: `int`

MurphyDecomposition `dataclass` ¶

MurphyDecomposition(reliability: float, resolution: float, uncertainty: float)

Binned Murphy (1973) partition of the Brier score.

reliability - resolution + uncertainty equals the Brier score exactly when predictions are constant within bins; otherwise the identity holds up to the within-bin variance of p (documented binning bias).

ATTRIBUTE	DESCRIPTION
`reliability`	Mean squared gap between within-bin predicted and observed rates. TYPE: `float`
`resolution`	Mean squared gap between within-bin observed rate and the overall base rate. TYPE: `float`
`uncertainty`	Base-rate variance `y_bar * (1 - y_bar)`. TYPE: `float`

EcceResult `dataclass` ¶

EcceResult(stat_max: float, stat_mean: float)

Empirical cumulative calibration error: Kolmogorov-style max and mean of the cumulative deviation over sorted predictions.

ATTRIBUTE	DESCRIPTION
`stat_max`	Maximum absolute cumulative deviation. TYPE: `float`
`stat_mean`	Weighted mean absolute cumulative deviation (each tied-score block weighted by its total sample weight). TYPE: `float`

SpiegelhalterResult `dataclass` ¶

SpiegelhalterResult(z: float, p_value: float)

Spiegelhalter's z test of forecast unbiasedness (two-sided).

ATTRIBUTE	DESCRIPTION
`z`	Standardized test statistic (`nan` when its null variance is 0). TYPE: `float`
`p_value`	Two-sided p-value under the standard normal approximation (`nan` when `z` is). TYPE: `float`

evaluate ¶

evaluate(y: object, p: object, *, sample_weight: object = None, n_boot: int = 1000, random_state: int = 42, metrics: Sequence[str] | None = None, stratify: bool = True, by: None = None) -> MetricReport

evaluate(y: object, p: object, *, sample_weight: object = None, n_boot: int = 1000, random_state: int = 42, metrics: Sequence[str] | None = None, stratify: bool = True, by: object) -> GroupedMetricReport

evaluate(y: object, p: object, *, sample_weight: object = None, n_boot: int = 1000, random_state: int = 42, metrics: Sequence[str] | None = None, stratify: bool = True, by: object = None) -> MetricReport | GroupedMetricReport

Full metric report with seeded bootstrap percentile confidence intervals.

PARAMETER	DESCRIPTION
`y`	Outcomes and predicted probabilities. TYPE: `array_like`
`p`	Outcomes and predicted probabilities. TYPE: `array_like`
`sample_weight`	Positive observation weights (resampled together with the observations); see the weight convention in :mod:`probcal.metrics._common`. TYPE: `array_like or None` DEFAULT: `None`
`n_boot`	Case-resampling bootstrap replicates (percentile CIs at 2.5/97.5). `0` skips the bootstrap: point estimates only, CI bounds `nan`. TYPE: `int` DEFAULT: `1000`
`random_state`	RNG seed; results are bit-reproducible given the seed. TYPE: `int` DEFAULT: `42`
`metrics`	Subset of catalog names to compute; `None` computes the full catalog. The report follows catalog order regardless of the order given here. TYPE: `sequence of str or None` DEFAULT: `None`
`stratify`	If `True` (default), each bootstrap replicate resamples the negative and positive classes separately (case resampling within strata), preserving the observed class counts exactly — the pROC-style default. This conditions the CI on the observed class balance: it removes the additional variance a plain i.i.d. bootstrap picks up from the replicate-to-replicate event count fluctuating, and it makes every replicate well-defined (never a single-class resample) on rare-event data, at the cost of not propagating sampling variance in the event rate itself. `y` must already contain both classes (checked unconditionally, independent of this flag). If `False`, replicates draw i.i.d. from all `n` rows; a degenerate (single-class) draw is redrawn up to 100 times before raising `RuntimeError`. TYPE: `bool` DEFAULT: `True`
`by`	Optional group labels, one per observation (same length as `y`). `None` (default) is the plain report above, unchanged. Otherwise rows are grouped by their raw label values (`1` and `"1"` are different groups) and a separate report is computed per group, in the order of the groups' display names `str(label)` — group `i` in that order is evaluated with `random_state + 1000 * i`, a fixed offset so results are reproducible independent of the label values or how many groups exist — plus a pooled report on the full data using `random_state` unchanged. Returns a :class:`~probcal._results.GroupedMetricReport` instead of a plain report. Group-conditional statistical testing (formal multiplicity-adjusted comparisons across groups) is out of scope here; see `docs/guide/groups.md`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`MetricReport or GroupedMetricReport`	Point estimates and CI bounds for the requested catalog (`by=None`, the default), or a pooled report plus one report per group (`by` given). Note the caveat from the metrics chapter: a bootstrap CI around a biased estimator (plain ECE) quantifies its variance, not its bias.

RAISES	DESCRIPTION
`ValueError`	If `metrics` contains names outside the metric catalog; if `n_boot` is negative; if `by` is given with a length that does not match `y`, or two distinct `by` values share a display name; or if a group has only one outcome class (the underlying `"y must contain both classes"` error, re-raised naming the group).
`RuntimeError`	If `stratify=False` and 100 consecutive bootstrap draws are all single-class.

Notes

Every metric is recomputed n_boot times; for n > 1e6 reduce n_boot or pass a metrics= subset (the ICI family's LOESS fit is the largest single cost). With by given, the cost is paid once per group plus once for the pooled report. Each replicate is sorted by prediction once, and the LOESS fit and the ECCE walk reuse that order; the point estimates always come off the unsorted path. The measured cost table lives in docs/concepts/metrics.md.

Examples:

>>> import numpy as np
>>> from probcal.metrics import evaluate
>>> rng = np.random.default_rng(0)
>>> p = rng.uniform(0.05, 0.5, 300)
>>> y = (rng.random(300) < p).astype(float)
>>> segment = np.where(p < 0.2, "low", "high")
>>> grouped = evaluate(y, p, n_boot=50, metrics=("brier",), by=segment)
>>> grouped.groups
('high', 'low')
>>> len(grouped.reports) == len(grouped.groups)
True

Source code in src/probcal/metrics/_evaluate.py

def evaluate(
    y: object,
    p: object,
    *,
    sample_weight: object = None,
    n_boot: int = 1000,
    random_state: int = 42,
    metrics: Sequence[str] | None = None,
    stratify: bool = True,
    by: object = None,
) -> MetricReport | GroupedMetricReport:
    """Full metric report with seeded bootstrap percentile confidence intervals.

    Parameters
    ----------
    y, p : array_like
        Outcomes and predicted probabilities.
    sample_weight : array_like or None
        Positive observation weights (resampled together with the
        observations); see the weight convention in
        :mod:`probcal.metrics._common`.
    n_boot : int
        Case-resampling bootstrap replicates (percentile CIs at 2.5/97.5).
        ``0`` skips the bootstrap: point estimates only, CI bounds ``nan``.
    random_state : int
        RNG seed; results are bit-reproducible given the seed.
    metrics : sequence of str or None
        Subset of catalog names to compute; ``None`` computes the full
        catalog. The report follows catalog order regardless of the order
        given here.
    stratify : bool
        If ``True`` (default), each bootstrap replicate resamples the
        negative and positive classes separately (case resampling within
        strata), preserving the observed class counts exactly — the
        pROC-style default. This conditions the CI on the observed class
        balance: it removes the additional variance a plain i.i.d. bootstrap
        picks up from the replicate-to-replicate event *count* fluctuating,
        and it makes every replicate well-defined (never a single-class
        resample) on rare-event data, at the cost of not propagating
        sampling variance in the event rate itself. ``y`` must already
        contain both classes (checked unconditionally, independent of this
        flag). If ``False``, replicates draw i.i.d. from all ``n`` rows; a
        degenerate (single-class) draw is redrawn up to 100 times before
        raising ``RuntimeError``.
    by : array_like or None, keyword-only
        Optional group labels, one per observation (same length as ``y``).
        ``None`` (default) is the plain report above, unchanged. Otherwise
        rows are grouped by their *raw* label values (``1`` and ``"1"`` are
        different groups) and a separate report is computed per group, in
        the order of the groups' display names ``str(label)`` — group ``i``
        in that order is evaluated with ``random_state + 1000 * i``, a fixed
        offset so results are reproducible independent of the label values
        or how many groups exist — plus a pooled report on the full data
        using ``random_state`` unchanged. Returns a
        :class:`~probcal._results.GroupedMetricReport` instead of a plain
        report. Group-conditional statistical *testing* (formal
        multiplicity-adjusted comparisons across groups) is out of scope
        here; see ``docs/guide/groups.md``.

    Returns
    -------
    MetricReport or GroupedMetricReport
        Point estimates and CI bounds for the requested catalog
        (``by=None``, the default), or a pooled report plus one report per
        group (``by`` given). Note the caveat from the metrics chapter: a
        bootstrap CI around a *biased* estimator (plain ECE) quantifies its
        variance, not its bias.

    Raises
    ------
    ValueError
        If ``metrics`` contains names outside the metric catalog; if
        ``n_boot`` is negative; if ``by`` is given with a length that does
        not match ``y``, or two distinct ``by`` values share a display name;
        or if a group has only one outcome class (the underlying
        ``"y must contain both classes"`` error, re-raised naming the
        group).
    RuntimeError
        If ``stratify=False`` and 100 consecutive bootstrap draws are all
        single-class.

    Notes
    -----
    Every metric is recomputed ``n_boot`` times; for n > 1e6 reduce
    ``n_boot`` or pass a ``metrics=`` subset (the ICI family's LOESS fit is
    the largest single cost). With ``by`` given, the cost is paid once per
    group plus once for the pooled report. Each replicate is sorted by
    prediction once, and the LOESS fit and the ECCE walk reuse that order;
    the point estimates always come off the unsorted path. The measured cost
    table lives in ``docs/concepts/metrics.md``.

    Examples
    --------
    >>> import numpy as np
    >>> from probcal.metrics import evaluate
    >>> rng = np.random.default_rng(0)
    >>> p = rng.uniform(0.05, 0.5, 300)
    >>> y = (rng.random(300) < p).astype(float)
    >>> segment = np.where(p < 0.2, "low", "high")
    >>> grouped = evaluate(y, p, n_boot=50, metrics=("brier",), by=segment)
    >>> grouped.groups
    ('high', 'low')
    >>> len(grouped.reports) == len(grouped.groups)
    True
    """
    n_boot = validate_nonneg_int(n_boot, "n_boot")
    if by is not None:
        return _evaluate_grouped(
            y,
            p,
            by,
            sample_weight=sample_weight,
            n_boot=n_boot,
            random_state=random_state,
            metrics=metrics,
            stratify=stratify,
        )
    return _evaluate(y, p, sample_weight, n_boot, random_state, metrics, stratify)

reliability_summary ¶

reliability_summary(y: object, p: object, *, sample_weight: object = None, grid_size: int | None = 512) -> ReliabilitySummary

Assemble the annotated-reliability stats box from existing metrics.

No new math: intercept and slope from the recalibration regression, ICI and E90 from the LOESS distance family, and Spiegelhalter's p-value. Lives here because, like evaluate, it aggregates across submodules; probcal.plots only formats the result. grid_size=None fits the LOESS exactly at every prediction.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`grid_size`	LOESS evaluation grid size for the ICI/E90 terms; `None` fits it exactly at every prediction. TYPE: `(int or None, keyword - only)` DEFAULT: `512`

RETURNS	DESCRIPTION
`ReliabilitySummary`	Stats-box fields for the annotated reliability diagram.

Source code in src/probcal/metrics/_summary.py

def reliability_summary(
    y: object,
    p: object,
    *,
    sample_weight: object = None,
    grid_size: int | None = 512,
) -> ReliabilitySummary:
    """Assemble the annotated-reliability stats box from existing metrics.

    No new math: intercept and slope from the recalibration regression, ICI
    and E90 from the LOESS distance family, and Spiegelhalter's p-value.
    Lives here because, like `evaluate`, it aggregates across submodules;
    ``probcal.plots`` only formats the result. ``grid_size=None`` fits the
    LOESS exactly at every prediction.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size for the ICI/E90 terms; ``None`` fits it
        exactly at every prediction.

    Returns
    -------
    ReliabilitySummary
        Stats-box fields for the annotated reliability diagram.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    fam = _ici_family(y_arr, p_arr, w, ("ici", "e90"), grid_size=grid_size)
    return ReliabilitySummary(
        n=len(y_arr),
        events=int(y_arr.sum()),
        intercept=_intercept(y_arr, p_arr, w),
        slope=float(_slope_fit(y_arr, p_arr, w, "reliability_summary")[1][1]),
        ici=fam["ici"],
        e90=fam["e90"],
        spiegelhalter_p=_spiegelhalter(y_arr, p_arr, w).p_value,
    )

ece ¶

ece(y: object, p: object, *, n_bins: int = 15, strategy: str = 'mass', norm: str = 'l1', sample_weight: object = None) -> float

Expected calibration error; norm="max" gives the MCE.

Binning-sensitive and upward-biased in finite samples — report, never select on it (see the metrics chapter's table).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`n_bins`	Requested number of bins (positive integer). TYPE: `(int, keyword - only)` DEFAULT: `15`
`strategy`	`"mass"` (equal-count, default; the "adaptive ECE" of the literature) or `"width"` (equal-width over [0, 1]). TYPE: `(mass, width)` DEFAULT: `"mass"`
`norm`	`"l1"` (default, the usual ECE), `"l2"`, or `"max"` (the MCE). TYPE: `(l1, l2, max)` DEFAULT: `"l1"`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Weighted binned calibration error under the chosen norm.

Source code in src/probcal/metrics/binned.py

def ece(
    y: object,
    p: object,
    *,
    n_bins: int = 15,
    strategy: str = "mass",
    norm: str = "l1",
    sample_weight: object = None,
) -> float:
    """Expected calibration error; ``norm="max"`` gives the MCE.

    Binning-sensitive and upward-biased in finite samples — report, never
    select on it (see the metrics chapter's table).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    n_bins : int, keyword-only
        Requested number of bins (positive integer).
    strategy : {"mass", "width"}, keyword-only
        ``"mass"`` (equal-count, default; the "adaptive ECE" of the
        literature) or ``"width"`` (equal-width over [0, 1]).
    norm : {"l1", "l2", "max"}, keyword-only
        ``"l1"`` (default, the usual ECE), ``"l2"``, or ``"max"`` (the MCE).
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Weighted binned calibration error under the chosen norm.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    bs = bin_stats(y_arr, p_arr, w, n_bins, strategy)
    return _reduce(*_shares_gaps(bs, float(w.sum())), norm)

ece_debiased ¶

ece_debiased(y: object, p: object, *, n_bins: int = 15, strategy: str = 'mass', sample_weight: object = None) -> float

Bias-corrected ECE, floored at zero.

Per-bin squared gaps minus the within-bin variance of the event rate (correction in the spirit of Bröcker 2009 / Ferro & Fricker 2012): sqrt(max(gap_b**2 - rate_b * (1 - rate_b) / (n_b - 1), 0)) for bins with n_b > 1. n_b is the bin's effective count sum_b(w)**2 / sum_b(w**2) -- the size of an unweighted bin whose event rate has the same variance -- so it follows the weights inside the bin, and is the raw count for unit (or within-bin constant) weights.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`n_bins`	Requested number of bins. TYPE: `(int, keyword - only)` DEFAULT: `15`
`strategy`	`"mass"` (equal-count, default) or `"width"` (equal-width over [0, 1]). TYPE: `(mass, width)` DEFAULT: `"mass"`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Bias-corrected calibration error.

Source code in src/probcal/metrics/binned.py

def ece_debiased(
    y: object,
    p: object,
    *,
    n_bins: int = 15,
    strategy: str = "mass",
    sample_weight: object = None,
) -> float:
    """Bias-corrected ECE, floored at zero.

    Per-bin squared gaps minus the within-bin variance of the event rate
    (correction in the spirit of Bröcker 2009 / Ferro & Fricker 2012):
    ``sqrt(max(gap_b**2 - rate_b * (1 - rate_b) / (n_b - 1), 0))`` for bins
    with ``n_b > 1``. ``n_b`` is the bin's effective count
    ``sum_b(w)**2 / sum_b(w**2)`` -- the size of an unweighted bin whose
    event rate has the same variance -- so it follows the weights inside the
    bin, and is the raw count for unit (or within-bin constant) weights.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    n_bins : int, keyword-only
        Requested number of bins.
    strategy : {"mass", "width"}, keyword-only
        ``"mass"`` (equal-count, default) or ``"width"`` (equal-width over [0, 1]).
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Bias-corrected calibration error.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    return _debiased(bin_stats(y_arr, p_arr, w, n_bins, strategy), w)

ece_sweep ¶

ece_sweep(y: object, p: object, *, norm: str = 'l1', sample_weight: object = None, max_bins: int = 100, rule: str = 'first_violation') -> float

Monotonic-sweep calibration error (Roelofs et al., 2022).

Equal-mass bins; the bin count is the largest b such that b and every smaller bin count give non-decreasing bin event rates — the paper's b* = max{b : for all b' <= b, ybar_1 <= ... <= ybar_b'} (its Algorithm 1: sweep b = 2, 3, ... and stop at the first non-monotone binning). The paper sweeps up to b = n; the sweep here stops at max_bins (default 100) as a cost cap.

rule="largest" instead takes the largest monotone b anywhere in 2..max_bins, skipping over non-monotone counts.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`norm`	Norm of the final binned error at the selected bin count. TYPE: `(l1, l2, max)` DEFAULT: `"l1"`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`max_bins`	Upper end of the sweep (capped at `n`). TYPE: `(int, keyword - only)` DEFAULT: `100`
`rule`	`"first_violation"` (default, the paper's rule) or `"largest"` (the largest monotone count, non-monotone counts skipped). TYPE: `(first_violation, largest)` DEFAULT: `"first_violation"`

RETURNS	DESCRIPTION
`float`	Calibration error at the selected bin count.

Source code in src/probcal/metrics/binned.py

def ece_sweep(
    y: object,
    p: object,
    *,
    norm: str = "l1",
    sample_weight: object = None,
    max_bins: int = 100,
    rule: str = "first_violation",
) -> float:
    """Monotonic-sweep calibration error (Roelofs et al., 2022).

    Equal-mass bins; the bin count is the largest ``b`` such that ``b`` *and
    every smaller bin count* give non-decreasing bin event rates — the
    paper's ``b* = max{b : for all b' <= b, ybar_1 <= ... <= ybar_b'}``
    (its Algorithm 1: sweep ``b = 2, 3, ...`` and stop at the first
    non-monotone binning). The paper sweeps up to ``b = n``; the sweep here
    stops at ``max_bins`` (default 100) as a cost cap.

    ``rule="largest"`` instead takes the largest monotone ``b`` anywhere in
    ``2..max_bins``, skipping over non-monotone counts.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    norm : {"l1", "l2", "max"}, keyword-only
        Norm of the final binned error at the selected bin count.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.
    max_bins : int, keyword-only
        Upper end of the sweep (capped at ``n``).
    rule : {"first_violation", "largest"}, keyword-only
        ``"first_violation"`` (default, the paper's rule) or ``"largest"``
        (the largest monotone count, non-monotone counts skipped).

    Returns
    -------
    float
        Calibration error at the selected bin count.
    """
    if rule not in _SWEEP_RULES:
        raise ValueError(f"rule must be one of {_SWEEP_RULES}, got {rule!r}")
    cap = validate_positive_int(max_bins, "max_bins")
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    bs = bin_stats(y_arr, p_arr, w, _sweep_best_b(y_arr, p_arr, w, cap, rule), "mass")
    return _reduce(*_shares_gaps(bs, float(w.sum())), norm)

hosmer_lemeshow ¶

hosmer_lemeshow(y: object, p: object, *, n_bins: int = 10, sample_weight: object = None, ddof: int = 0) -> HosmerLemeshowResult

Hosmer–Lemeshow goodness-of-fit test on n_bins equal-mass risk groups.

The statistic depends on an essentially arbitrary grouping and its power scales with n — see the metrics chapter for why this is report-only.

Per group g, O_g = sum(w y) and E_g = sum(w p) are compared against V_g = S_g * pt_g * (1 - pt_g), where S_g = sum(w**2) and pt_g = sum(w**2 p) / S_g. V_g is the weighted form of the classical n_g pbar_g (1 - pbar_g): it equals the exact null variance sum(w**2 p (1 - p)) of O_g - E_g plus the same within-group spread term sum(w**2 (p - pt_g)**2) the unweighted statistic carries. The statistic is invariant to rescaling the weights, and unit weights give the unweighted test bit for bit.

The reference distribution is chi-square on groups - ddof degrees of freedom. With p fixed in advance (a model or calibrator fitted on other data, the usual validation setting) the per-group standardized residuals are asymptotically independent, so the statistic is chi2(groups): ddof=0, the default (Hosmer, Lemeshow and Sturdivant, 2013, on external validation). The classical groups - 2 is for a logistic model fitted on the same data; pass ddof=2 for it. On out-of-sample p it is anti-conservative (about 12% rejections at a nominal 5% with 10 groups).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`n_bins`	Requested number of equal-mass risk groups. Tied scores can collapse groups; when no more than `ddof` groups are non-empty the test has no degrees of freedom and returns `df=0`, `p_value=nan`. TYPE: `(int, keyword - only)` DEFAULT: `10`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`ddof`	Degrees of freedom subtracted from the number of non-empty groups: `0` (default) for `p` fixed before seeing `y`, `2` for a logistic model fitted on the same `(y, p)` data. TYPE: `(int, keyword - only)` DEFAULT: `0`

RETURNS	DESCRIPTION
`HosmerLemeshowResult`	Chi-square statistic, degrees of freedom, and p-value.

Source code in src/probcal/metrics/binned.py

def hosmer_lemeshow(
    y: object,
    p: object,
    *,
    n_bins: int = 10,
    sample_weight: object = None,
    ddof: int = 0,
) -> HosmerLemeshowResult:
    """Hosmer–Lemeshow goodness-of-fit test on ``n_bins`` equal-mass risk groups.

    The statistic depends on an essentially arbitrary grouping and its power
    scales with n — see the metrics chapter for why this is report-only.

    Per group ``g``, ``O_g = sum(w y)`` and ``E_g = sum(w p)`` are compared
    against ``V_g = S_g * pt_g * (1 - pt_g)``, where ``S_g = sum(w**2)`` and
    ``pt_g = sum(w**2 p) / S_g``. ``V_g`` is the weighted form of the
    classical ``n_g pbar_g (1 - pbar_g)``: it equals the exact null variance
    ``sum(w**2 p (1 - p))`` of ``O_g - E_g`` plus the same within-group
    spread term ``sum(w**2 (p - pt_g)**2)`` the unweighted statistic carries.
    The statistic is invariant to rescaling the weights, and unit weights
    give the unweighted test bit for bit.

    The reference distribution is chi-square on ``groups - ddof`` degrees of
    freedom. With ``p`` fixed in advance (a model or calibrator fitted on
    other data, the usual validation setting) the per-group standardized
    residuals are asymptotically independent, so the statistic is
    ``chi2(groups)``: ``ddof=0``, the default (Hosmer, Lemeshow and
    Sturdivant, 2013, on external validation). The classical ``groups - 2``
    is for a logistic model fitted on the same data; pass ``ddof=2`` for it.
    On out-of-sample ``p`` it is anti-conservative (about 12% rejections at
    a nominal 5% with 10 groups).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    n_bins : int, keyword-only
        Requested number of equal-mass risk groups. Tied scores can collapse
        groups; when no more than ``ddof`` groups are non-empty the test has
        no degrees of freedom and returns ``df=0``, ``p_value=nan``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.
    ddof : int, keyword-only
        Degrees of freedom subtracted from the number of non-empty groups:
        ``0`` (default) for ``p`` fixed before seeing ``y``, ``2`` for a
        logistic model fitted on the same ``(y, p)`` data.

    Returns
    -------
    HosmerLemeshowResult
        Chi-square statistic, degrees of freedom, and p-value.
    """
    ddof = validate_nonneg_int(ddof, "ddof")
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    idx, m = bin_index(p_arr, n_bins, "mass")
    keep = np.bincount(idx, minlength=m) > 0
    w2 = w * w
    obs = np.bincount(idx, weights=w * y_arr, minlength=m)[keep]
    exp = np.bincount(idx, weights=w * p_arr, minlength=m)[keep]
    s2 = np.bincount(idx, weights=w2, minlength=m)[keep]
    s2p = np.bincount(idx, weights=w2 * p_arr, minlength=m)[keep]
    denom = s2p * (1.0 - s2p / s2)
    ok = denom > 0
    stat = float(np.sum((obs[ok] - exp[ok]) ** 2 / denom[ok]))
    df = max(len(s2) - ddof, 0)
    p_value = chi2_sf(stat, df) if df > 0 else math.nan
    return HosmerLemeshowResult(statistic=stat, df=df, p_value=p_value)

binomial_grade_test ¶

binomial_grade_test(y: object, p: object, grades: object, *, sample_weight: object = None, order: object = None) -> BinomialGradeResult

Exact binomial tail test per grade: P(X >= k | n, PD).

Small p-values flag grades with more defaults than the assigned PD supports. The exact tail uses the incomplete-beta identity P(X >= k) = I_PD(k, n - k + 1); the normal approximation is reported alongside. Traffic lights: green > 0.05, amber > 0.01, red <= 0.01. ci_low/ci_high are 90% Clopper-Pearson display intervals for the observed rate; the traffic light itself remains the one-sided exact test, unchanged.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities (assigned PDs) in `[0, 1]`. TYPE: `array_like`
`grades`	Rating grade label per observation, or a :class:`probcal.Masterscale` that assigns them from `p` (results then come out best to worst). TYPE: `array_like or Masterscale`
`sample_weight`	Not used: grade tests use raw integer counts. A `UserWarning` is emitted if the weights are non-uniform. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`order`	Explicit grade order for the result (a permutation of the grade labels present); `None` keeps the masterscale order or sorts the labels. TYPE: `sequence of str or None, keyword-only` DEFAULT: `None`

RETURNS	DESCRIPTION
`BinomialGradeResult`	Per-grade counts, p-values, traffic lights, and display intervals.

Source code in src/probcal/metrics/grade.py

def binomial_grade_test(
    y: object,
    p: object,
    grades: object,
    *,
    sample_weight: object = None,
    order: object = None,
) -> BinomialGradeResult:
    """Exact binomial tail test per grade: P(X >= k | n, PD).

    Small p-values flag grades with more defaults than the assigned PD
    supports. The exact tail uses the incomplete-beta identity
    ``P(X >= k) = I_PD(k, n - k + 1)``; the normal approximation is reported
    alongside. Traffic lights: green > 0.05, amber > 0.01, red <= 0.01.
    ``ci_low``/``ci_high`` are 90% Clopper-Pearson display intervals for the
    observed rate; the traffic light itself remains the one-sided exact test,
    unchanged.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities (assigned PDs) in ``[0, 1]``.
    grades : array_like or Masterscale
        Rating grade label per observation, or a :class:`probcal.Masterscale`
        that assigns them from ``p`` (results then come out best to worst).
    sample_weight : array_like or None, keyword-only
        Not used: grade tests use raw integer counts. A ``UserWarning`` is
        emitted if the weights are non-uniform.
    order : sequence of str or None, keyword-only
        Explicit grade order for the result (a permutation of the grade
        labels present); ``None`` keeps the masterscale order or sorts the
        labels.

    Returns
    -------
    BinomialGradeResult
        Per-grade counts, p-values, traffic lights, and display intervals.
    """
    labels, n, k, pd = _grade_counts(y, p, grades, sample_weight, order)
    p_exact, p_normal, ci_low, ci_high = (np.empty(len(labels)) for _ in range(4))
    for i in range(len(labels)):
        ki, ni = int(k[i]), int(n[i])
        p_exact[i] = 1.0 if ki == 0 else float(betainc(float(ki), float(ni - ki + 1), pd[i]))
        se = np.sqrt(n[i] * pd[i] * (1.0 - pd[i]))
        p_normal[i] = norm_sf((k[i] - n[i] * pd[i]) / se if se > 0 else 0.0)
        ci_low[i] = 0.0 if ki == 0 else beta_ppf(0.05, float(ki), float(ni - ki + 1))
        ci_high[i] = 1.0 if ki == ni else beta_ppf(0.95, float(ki + 1), float(ni - ki))
    light = tuple(_traffic_light(v) for v in p_exact)
    return BinomialGradeResult(
        grades=labels,
        n=n,
        k=k,
        pd=pd,
        p_exact=p_exact,
        p_normal=p_normal,
        light=light,
        ci_low=ci_low,
        ci_high=ci_high,
    )

hl_e_test ¶

hl_e_test(y: object, p: object, grades: object, *, mixture_grid: tuple[float, ...] = (0.1, 0.25, 0.5, 1.0), sample_weight: object = None) -> HlEResult

Fixed-sample mixture-LR grade e-test (safe Hosmer-Lemeshow analogue).

For grade g, with z_i = logit(p_i):

log E_g = logsumexp_{delta in +/-mixture_grid}(sum_{i in g} log LR_i(sigma(z_i + delta) : p_i)) - log(2 * len(mixture_grid))

i.e. the log-mean Bernoulli log-likelihood-ratio (probcal._math .bern_log_lr) of the grade's observations, averaged over the symmetrized offset grid -- the same mixture construction CalibrationMonitor's offset e-process uses, applied once per grade with no predictable (plug-in) component, since a fixed sample has no strictly-earlier data to learn one from. The test statistic is the product across grades, log E = sum_g log E_g, e_value = exp(log E): grades partition the sample into disjoint observations, each grade's mixture average is an e-value for that grade's null (an average of e-values, each with conditional expectation 1 under H0), and the product of e-values over independent (here: disjoint-observation) factors is itself an e-value. p_value = min(1, 1 / e_value) follows from Markov's inequality and is a valid (generally conservative) p-value.

Sample weights, when given, enter as exponents min(1, w / mean(w)) (mean over the grade) on the Bernoulli factors, as in CalibrationMonitor: every factor stays an e-value for any weights, the test is invariant to rescaling the weights, and unit weights give the unweighted test.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`; both classes must be present overall (grade-level all-zero/all-one subsets are fine). TYPE: `array_like`
`p`	Assigned probabilities (the null) in `[0, 1]`. TYPE: `array_like`
`grades`	Rating grade label per observation, or a :class:`probcal.Masterscale` that assigns them from `p` (grades then come out best to worst). TYPE: `array_like or Masterscale`
`mixture_grid`	Positive logit-scale offsets; symmetrized to `+/-` before averaging (matching `CalibrationMonitor(mixture_grid=...)`). TYPE: `tuple of float, keyword-only` DEFAULT: `(0.1, 0.25, 0.5, 1.0)`
`sample_weight`	Optional positive weights, same length as `y`. Each grade's likelihood ratios are raised to `min(1, w / mean(w))` over that grade, which keeps every factor an e-value for any weights. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`HlEResult`	Per-grade and combined e-values, the derived p-value, and the construction tag.

RAISES	DESCRIPTION
`ValueError`	If `y`/`p`/`grades` are not equal-length 1-D arrays, `y` is not binary or single-class, `p` lies outside `[0, 1]`, or `mixture_grid` is empty or not 1-D.

Examples:

>>> import numpy as np
>>> from probcal.metrics import hl_e_test
>>> rng = np.random.default_rng(1)
>>> p = np.full(400, 0.05)
>>> y = (rng.random(400) < 0.05).astype(float)
>>> grades = np.array(["A"] * 200 + ["B"] * 200)
>>> res = hl_e_test(y, p, grades)
>>> res.grades
('A', 'B')
>>> res.e_value > 0.0
True

Source code in src/probcal/metrics/_safe.py

def hl_e_test(
    y: object,
    p: object,
    grades: object,
    *,
    mixture_grid: tuple[float, ...] = (0.1, 0.25, 0.5, 1.0),
    sample_weight: object = None,
) -> HlEResult:
    """Fixed-sample mixture-LR grade e-test (safe Hosmer-Lemeshow analogue).

    For grade ``g``, with ``z_i = logit(p_i)``:

    ``log E_g = logsumexp_{delta in +/-mixture_grid}(sum_{i in g} log
    LR_i(sigma(z_i + delta) : p_i)) - log(2 * len(mixture_grid))``

    i.e. the log-mean Bernoulli log-likelihood-ratio (``probcal._math
    .bern_log_lr``) of the grade's observations, averaged over the
    symmetrized offset grid -- the same mixture construction
    ``CalibrationMonitor``'s offset e-process uses, applied once per grade
    with no predictable (plug-in) component, since a fixed sample has no
    strictly-earlier data to learn one from. The test statistic is the
    product across grades, ``log E = sum_g log E_g``, ``e_value = exp(log
    E)``: grades partition the sample into disjoint observations, each
    grade's mixture average is an e-value for that grade's null (an average
    of e-values, each with conditional expectation 1 under H0), and the
    product of e-values over independent (here: disjoint-observation)
    factors is itself an e-value. ``p_value = min(1, 1 / e_value)`` follows
    from Markov's inequality and is a valid (generally conservative)
    p-value.

    Sample weights, when given, enter as exponents ``min(1, w / mean(w))``
    (mean over the grade) on the Bernoulli factors, as in
    ``CalibrationMonitor``: every factor stays an e-value for any weights,
    the test is invariant to rescaling the weights, and unit weights give
    the unweighted test.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``; both classes must be present overall
        (grade-level all-zero/all-one subsets are fine).
    p : array_like
        Assigned probabilities (the null) in ``[0, 1]``.
    grades : array_like or Masterscale
        Rating grade label per observation, or a :class:`probcal.Masterscale`
        that assigns them from ``p`` (grades then come out best to worst).
    mixture_grid : tuple of float, keyword-only
        Positive logit-scale offsets; symmetrized to ``+/-`` before
        averaging (matching ``CalibrationMonitor(mixture_grid=...)``).
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``. Each grade's
        likelihood ratios are raised to ``min(1, w / mean(w))`` over that
        grade, which keeps every factor an e-value for any weights.

    Returns
    -------
    HlEResult
        Per-grade and combined e-values, the derived p-value, and the
        construction tag.

    Raises
    ------
    ValueError
        If ``y``/``p``/``grades`` are not equal-length 1-D arrays, ``y``
        is not binary or single-class, ``p`` lies outside ``[0, 1]``, or
        ``mixture_grid`` is empty or not 1-D.

    Examples
    --------
    >>> import numpy as np
    >>> from probcal.metrics import hl_e_test
    >>> rng = np.random.default_rng(1)
    >>> p = np.full(400, 0.05)
    >>> y = (rng.random(400) < 0.05).astype(float)
    >>> grades = np.array(["A"] * 200 + ["B"] * 200)
    >>> res = hl_e_test(y, p, grades)
    >>> res.grades
    ('A', 'B')
    >>> res.e_value > 0.0
    True
    """
    y_arr, p_arr, w_arr = _prep(y, p, sample_weight)
    g_arr, default_order = grade_labels(grades, p_arr, len(y_arr))
    grid = np.asarray(mixture_grid, dtype=np.float64)
    if grid.ndim != 1 or grid.size == 0:
        raise ValueError("mixture_grid must be a non-empty 1-D sequence")
    if not np.all(np.isfinite(grid)):
        raise ValueError("mixture_grid must contain only finite values")

    offsets = np.concatenate([grid, -grid])
    log_norm = float(np.log(2.0 * grid.size))
    z_arr = logit(p_arr)

    labels = resolve_order(g_arr, None, default_order)
    log_e_grade = np.empty(len(labels))
    for i, label in enumerate(labels):
        mask = g_arr == label
        y_g, p_g, z_g = y_arr[mask], p_arr[mask], z_arr[mask]
        w_g = exponent_weights(w_arr[mask])
        log_terms = np.array([bern_log_lr(y_g, p_g, expit(z_g + delta), w_g) for delta in offsets])
        log_e_grade[i] = logsumexp(log_terms) - log_norm

    log_e = float(np.sum(log_e_grade))
    e_grade = np.exp(log_e_grade)
    e_value = float(np.exp(log_e))
    p_value = min(1.0, 1.0 / e_value) if e_value > 0.0 else 1.0

    return HlEResult(
        e_value=e_value,
        p_value=p_value,
        grades=labels,
        e_grade=e_grade,
        construction="mixture-lr",
    )

jeffreys_grade_test ¶

jeffreys_grade_test(y: object, p: object, grades: object, *, sample_weight: object = None, order: object = None) -> JeffreysGradeResult

Jeffreys test per grade: posterior P(theta <= PD | k, n) under Beta(k+1/2, n-k+1/2).

One-sided and conservative by design: a small value flags a grade whose PD is likely understated. Do not read it two-sided (a recurring validation error — see the metrics chapter). ci_low/ci_high are the central 90% Jeffreys posterior display intervals; the traffic light itself remains the one-sided posterior test, unchanged.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities (assigned PDs) in `[0, 1]`. TYPE: `array_like`
`grades`	Rating grade label per observation, or a :class:`probcal.Masterscale` that assigns them from `p` (results then come out best to worst). TYPE: `array_like or Masterscale`
`sample_weight`	Not used: grade tests use raw integer counts. A `UserWarning` is emitted if the weights are non-uniform. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`order`	Explicit grade order for the result (a permutation of the grade labels present); `None` keeps the masterscale order or sorts the labels. TYPE: `sequence of str or None, keyword-only` DEFAULT: `None`

RETURNS	DESCRIPTION
`JeffreysGradeResult`	Per-grade counts, p-values, traffic lights, and display intervals.

Source code in src/probcal/metrics/grade.py

def jeffreys_grade_test(
    y: object,
    p: object,
    grades: object,
    *,
    sample_weight: object = None,
    order: object = None,
) -> JeffreysGradeResult:
    """Jeffreys test per grade: posterior P(theta <= PD | k, n) under Beta(k+1/2, n-k+1/2).

    One-sided and conservative by design: a small value flags a grade whose
    PD is likely understated. Do not read it two-sided (a recurring
    validation error — see the metrics chapter). ``ci_low``/``ci_high`` are
    the central 90% Jeffreys posterior display intervals; the traffic light
    itself remains the one-sided posterior test, unchanged.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities (assigned PDs) in ``[0, 1]``.
    grades : array_like or Masterscale
        Rating grade label per observation, or a :class:`probcal.Masterscale`
        that assigns them from ``p`` (results then come out best to worst).
    sample_weight : array_like or None, keyword-only
        Not used: grade tests use raw integer counts. A ``UserWarning`` is
        emitted if the weights are non-uniform.
    order : sequence of str or None, keyword-only
        Explicit grade order for the result (a permutation of the grade
        labels present); ``None`` keeps the masterscale order or sorts the
        labels.

    Returns
    -------
    JeffreysGradeResult
        Per-grade counts, p-values, traffic lights, and display intervals.
    """
    labels, n, k, pd = _grade_counts(y, p, grades, sample_weight, order)
    p_value, ci_low, ci_high = (np.empty(len(labels)) for _ in range(3))
    for i in range(len(labels)):
        a, b = k[i] + 0.5, n[i] - k[i] + 0.5
        p_value[i] = float(betainc(a, b, pd[i]))
        ci_low[i] = beta_ppf(0.05, a, b)
        ci_high[i] = beta_ppf(0.95, a, b)
    light = tuple(_traffic_light(v) for v in p_value)
    return JeffreysGradeResult(
        grades=labels,
        n=n,
        k=k,
        pd=pd,
        p_value=p_value,
        light=light,
        ci_low=ci_low,
        ci_high=ci_high,
    )

jeffreys_upper_bands ¶

jeffreys_upper_bands(y: object, p: object, grades: object, *, confidence: float = 0.9, order: object = None) -> dict[str, tuple[float, float]]

Jeffreys per-grade upper bounds as a contiguous masterscale band table.

For grade i (best to worst, in order): hi_i is the same one-sided Jeffreys posterior upper bound jeffreys_grade_test reports as its own-grade display interval, beta_ppf(confidence, k_i + 0.5, n_i - k_i + 0.5) under a Beta(k_i + 0.5, n_i - k_i + 0.5) posterior on grade i's own default rate; lo_i is the previous grade's hi (0.0 for the best grade), so the bands are contiguous by construction: (lo_0, hi_0), (hi_0, hi_1), (hi_1, hi_2), .... Unlike :func:pluto_tasche, each grade's bound uses only its own counts (no pooling across grades), so a zero-default grade still gets a strictly positive hi from the Jeffreys prior alone.

The raw hi sequence need not come out increasing (a noisy grade can have a posterior upper bound no larger than a better grade's), which would make the bands overlap or invert. The prudent hull of :func:pluto_tasche -- the running maximum best to worst, which never lowers any grade's bound -- fixes the order but would leave such a grade a zero-width band (hi, hi) that no calibrated PD can fall into. Instead, a grade whose own bound does not exceed the running maximum is merged into the band before it: the merged band is keyed by the joined labels ("A+B", best to worst), keeps the running-maximum hi (so no grade's bound is ever lowered), and a UserWarning names the merged grades. The data cannot order those grades at this confidence; merging them (or collecting more data) is the decision the warning asks for. Every returned band therefore has lo < hi, and the bands are contiguous and strictly increasing. A size-weighted PAVA fit is not used because it can lower a grade's upper bound.

The resulting {grade: (lo, hi)} table is exactly the shape :func:probcal.thresholds.calibrated_bands_to_raw and :class:probcal.Masterscale consume to translate a masterscale defined on calibrated PD into raw-score intervals.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. Unlike most probcal metrics (but like :func:`pluto_tasche_from_arrays`), an all-zero `y` is accepted -- the low- and zero-default portfolio is this module's motivating case. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]` (only used, alongside `grades`, to determine the default best-to-worst `order` when `order` is not given). TYPE: `array_like`
`grades`	Rating grade label per observation, or a :class:`probcal.Masterscale` that assigns them from `p` and supplies the default `order`. TYPE: `array_like or Masterscale`
`confidence`	Confidence level in `(0, 1)` for every grade's Jeffreys upper bound. TYPE: `(float, keyword - only)` DEFAULT: `0.9`
`order`	Explicit best-to-worst grade order; must match the unique labels in `grades` exactly. `None` (default) uses the masterscale's order, or orders label-array grades by their mean `p`, ascending (lowest predicted PD first). TYPE: `sequence of str or None, keyword-only` DEFAULT: `None`

RETURNS	DESCRIPTION
`dict[str, tuple[float, float]]`	Mapping of grade label -- or `"+"`-joined labels of merged grades -- to `(lo, hi)` calibrated-probability bounds, contiguous and strictly increasing best to worst, every band with `lo < hi`.

RAISES	DESCRIPTION
`ValueError`	If `y`/`p`/`grades` are not equal-length 1-D arrays, `confidence` is not in `(0, 1)`, or `order` does not match the unique grade labels.

Examples:

>>> import numpy as np
>>> from probcal.metrics import jeffreys_upper_bands
>>> grades = np.array(["A"] * 100 + ["B"] * 100)
>>> y = np.array([0.0] * 100 + [1.0] * 5 + [0.0] * 95)
>>> p = np.array([0.01] * 100 + [0.05] * 100)
>>> bands = jeffreys_upper_bands(y, p, grades, confidence=0.9)
>>> bands["A"][0]
0.0
>>> bands["A"][1] < bands["B"][1]
True

Source code in src/probcal/metrics/_conservative.py

def jeffreys_upper_bands(
    y: object,
    p: object,
    grades: object,
    *,
    confidence: float = 0.9,
    order: object = None,
) -> dict[str, tuple[float, float]]:
    """Jeffreys per-grade upper bounds as a contiguous masterscale band table.

    For grade ``i`` (best to worst, in ``order``): ``hi_i`` is the same
    one-sided Jeffreys posterior upper bound ``jeffreys_grade_test`` reports
    as its own-grade display interval, ``beta_ppf(confidence, k_i + 0.5,
    n_i - k_i + 0.5)`` under a ``Beta(k_i + 0.5, n_i - k_i + 0.5)`` posterior
    on grade ``i``'s own default rate; ``lo_i`` is the previous grade's
    ``hi`` (``0.0`` for the best grade), so the bands are contiguous by
    construction: ``(lo_0, hi_0), (hi_0, hi_1), (hi_1, hi_2), ...``. Unlike
    :func:`pluto_tasche`, each grade's bound uses only its own counts (no
    pooling across grades), so a zero-default grade still gets a strictly
    positive ``hi`` from the Jeffreys prior alone.

    The raw ``hi`` sequence need not come out increasing (a noisy grade can
    have a posterior upper bound no larger than a better grade's), which
    would make the bands overlap or invert. The prudent hull of
    :func:`pluto_tasche` -- the running maximum best to worst, which never
    lowers any grade's bound -- fixes the order but would leave such a grade
    a zero-width band ``(hi, hi)`` that no calibrated PD can fall into.
    Instead, a grade whose own bound does not exceed the running maximum is
    merged into the band before it: the merged band is keyed by the joined
    labels (``"A+B"``, best to worst), keeps the running-maximum ``hi`` (so
    no grade's bound is ever lowered), and a ``UserWarning`` names the
    merged grades. The data cannot order those grades at this confidence;
    merging them (or collecting more data) is the decision the warning asks
    for. Every returned band therefore has ``lo < hi``, and the bands are
    contiguous and strictly increasing. A size-weighted PAVA fit is not
    used because it can *lower* a grade's upper bound.

    The resulting ``{grade: (lo, hi)}`` table is exactly the shape
    :func:`probcal.thresholds.calibrated_bands_to_raw` and
    :class:`probcal.Masterscale` consume to translate a masterscale defined
    on calibrated PD into raw-score intervals.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``. Unlike most probcal metrics (but like
        :func:`pluto_tasche_from_arrays`), an all-zero ``y`` is accepted --
        the low- and zero-default portfolio is this module's motivating
        case.
    p : array_like
        Predicted probabilities in ``[0, 1]`` (only used, alongside
        ``grades``, to determine the default best-to-worst ``order`` when
        ``order`` is not given).
    grades : array_like or Masterscale
        Rating grade label per observation, or a :class:`probcal.Masterscale`
        that assigns them from ``p`` and supplies the default ``order``.
    confidence : float, keyword-only
        Confidence level in ``(0, 1)`` for every grade's Jeffreys upper
        bound.
    order : sequence of str or None, keyword-only
        Explicit best-to-worst grade order; must match the unique labels in
        ``grades`` exactly. ``None`` (default) uses the masterscale's order,
        or orders label-array grades by their mean ``p``, ascending (lowest
        predicted PD first).

    Returns
    -------
    dict[str, tuple[float, float]]
        Mapping of grade label -- or ``"+"``-joined labels of merged grades
        -- to ``(lo, hi)`` calibrated-probability bounds, contiguous and
        strictly increasing best to worst, every band with ``lo < hi``.

    Raises
    ------
    ValueError
        If ``y``/``p``/``grades`` are not equal-length 1-D arrays,
        ``confidence`` is not in ``(0, 1)``, or ``order`` does not match the unique grade
        labels.

    Examples
    --------
    >>> import numpy as np
    >>> from probcal.metrics import jeffreys_upper_bands
    >>> grades = np.array(["A"] * 100 + ["B"] * 100)
    >>> y = np.array([0.0] * 100 + [1.0] * 5 + [0.0] * 95)
    >>> p = np.array([0.01] * 100 + [0.05] * 100)
    >>> bands = jeffreys_upper_bands(y, p, grades, confidence=0.9)
    >>> bands["A"][0]
    0.0
    >>> bands["A"][1] < bands["B"][1]
    True
    """
    y_arr = validate_binary_y(y, require_both_classes=False)
    p_arr = validate_scores(p, name="p")
    if len(p_arr) != len(y_arr):
        raise ValueError("y and p must have equal length")
    g_arr, default_order = grade_labels(grades, p_arr, len(y_arr))
    conf = validate_level(confidence, "confidence")

    order_t = resolve_order(g_arr, order, default_order)
    n, k, mean_p, _ = aggregate(g_arr, order_t, y_arr, p=p_arr)
    if order is None and default_order is None:  # label array: order by mean p
        assert mean_p is not None
        perm = np.argsort(mean_p, kind="stable")
        order_t, n, k = tuple(order_t[i] for i in perm), n[perm], k[perm]

    hi_raw = np.array([beta_ppf(conf, k[i] + 0.5, n[i] - k[i] + 0.5) for i in range(len(order_t))])
    # Running maximum, merging every grade it would leave with a zero-width
    # band into the preceding band.
    blocks: list[tuple[list[str], float]] = []
    for label, h in zip(order_t, hi_raw, strict=True):
        if blocks and h <= blocks[-1][1]:
            blocks[-1][0].append(str(label))
        else:
            blocks.append(([str(label)], float(h)))
    merged = [labels for labels, _ in blocks if len(labels) > 1]
    if merged:
        warnings.warn(
            "jeffreys_upper_bands: the Jeffreys upper bound of grade(s) "
            + ", ".join(repr(lab) for labels in merged for lab in labels[1:])
            + " does not exceed the running maximum of the better grades; merged into "
            + ", ".join(repr("+".join(labels)) for labels in merged)
            + " so that every band has positive width",
            UserWarning,
            stacklevel=2,
        )
    out: dict[str, tuple[float, float]] = {}
    lo = 0.0
    for labels, h in blocks:
        out["+".join(labels)] = (lo, h)
        lo = h
    return out

pluto_tasche ¶

pluto_tasche(grade_n: object, grade_d: object, *, confidence: float = 0.9, grades: object = None) -> PlutoTascheResult

Pluto & Tasche (2005) one-period most-prudent PD, from per-grade counts.

For grade i (best to worst, in the order given), pool its own obligors and defaults with every worse grade's: n*_i = sum(n[i:]), d*_i = sum(d[i:]). The most-prudent PD is the one-sided Clopper-Pearson upper bound of the pooled rate, p solving I_p(d*_i + 1, n*_i - d*_i) = confidence (beta_ppf(confidence, d*_i + 1, n*_i - d*_i)), i.e. the largest PD under which observing at most d*_i defaults in n*_i obligors still has probability >= 1 - confidence. Pooling with worse grades is the rating-monotonicity assumption doing its work: a grade's own data alone is often uninformative (frequently zero defaults), but the assumption that its true PD cannot exceed a worse grade's lets that grade's defaults bound this one.

PARAMETER	DESCRIPTION
`grade_n`	Obligor count per grade, best to worst. Non-integer (weighted) counts are accepted and pass directly into the Beta shape parameters below. TYPE: `array_like`
`grade_d`	Default count per grade, same order; `grade_d[i] <= grade_n[i]`. TYPE: `array_like`
`confidence`	Confidence level in `(0, 1)` for every grade's upper bound. TYPE: `(float, keyword - only)` DEFAULT: `0.9`
`grades`	Grade labels, best to worst; `None` uses `"1", "2", ..., "K"`. TYPE: `sequence of str or None, keyword-only` DEFAULT: `None`

RETURNS	DESCRIPTION
`PlutoTascheResult`	Per-grade counts, pooled counts, and most-prudent PDs.

RAISES	DESCRIPTION
`ValueError`	If `grade_n`/`grade_d` are not equal-length 1-D arrays, contain negative values, have a default count exceeding the obligor count, `confidence` is not in `(0, 1)`, `grades` does not match the count arrays' length, or a grade's pooled obligor count (`n*_i`) is zero.

Examples:

>>> import numpy as np
>>> from probcal.metrics import pluto_tasche
>>> res = pluto_tasche(
...     np.array([100.0, 400.0, 300.0]),
...     np.array([0.0, 0.0, 0.0]),
...     confidence=0.9,
...     grades=("A", "B", "C"),
... )
>>> res.grades
('A', 'B', 'C')
>>> np.round(res.pd_upper, 4)
array([0.0029, 0.0033, 0.0076])

Source code in src/probcal/metrics/_conservative.py

def pluto_tasche(
    grade_n: object,
    grade_d: object,
    *,
    confidence: float = 0.9,
    grades: object = None,
) -> PlutoTascheResult:
    """Pluto & Tasche (2005) one-period most-prudent PD, from per-grade counts.

    For grade ``i`` (best to worst, in the order given), pool its own
    obligors and defaults with every worse grade's: ``n*_i = sum(n[i:])``,
    ``d*_i = sum(d[i:])``. The most-prudent PD is the one-sided
    Clopper-Pearson upper bound of the pooled rate,
    ``p`` solving ``I_p(d*_i + 1, n*_i - d*_i) = confidence``
    (``beta_ppf(confidence, d*_i + 1, n*_i - d*_i)``), i.e. the largest PD
    under which observing at most ``d*_i`` defaults in ``n*_i`` obligors
    still has probability ``>= 1 - confidence``. Pooling with worse grades is the
    rating-monotonicity assumption doing its work: a grade's own data alone
    is often uninformative (frequently zero defaults), but the assumption
    that its true PD cannot exceed a worse grade's lets that grade's
    defaults bound this one.

    Parameters
    ----------
    grade_n : array_like
        Obligor count per grade, best to worst. Non-integer (weighted)
        counts are accepted and pass directly into the Beta shape
        parameters below.
    grade_d : array_like
        Default count per grade, same order; ``grade_d[i] <= grade_n[i]``.
    confidence : float, keyword-only
        Confidence level in ``(0, 1)`` for every grade's upper bound.
    grades : sequence of str or None, keyword-only
        Grade labels, best to worst; ``None`` uses ``"1", "2", ..., "K"``.

    Returns
    -------
    PlutoTascheResult
        Per-grade counts, pooled counts, and most-prudent PDs.

    Raises
    ------
    ValueError
        If ``grade_n``/``grade_d`` are not equal-length 1-D arrays, contain
        negative values, have a default count exceeding the obligor count,
        ``confidence`` is not in ``(0, 1)``, ``grades`` does not match the
        count arrays' length, or a grade's pooled obligor count
        (``n*_i``) is zero.

    Examples
    --------
    >>> import numpy as np
    >>> from probcal.metrics import pluto_tasche
    >>> res = pluto_tasche(
    ...     np.array([100.0, 400.0, 300.0]),
    ...     np.array([0.0, 0.0, 0.0]),
    ...     confidence=0.9,
    ...     grades=("A", "B", "C"),
    ... )
    >>> res.grades
    ('A', 'B', 'C')
    >>> np.round(res.pd_upper, 4)
    array([0.0029, 0.0033, 0.0076])
    """
    n = np.asarray(grade_n, dtype=np.float64)
    d = np.asarray(grade_d, dtype=np.float64)
    if n.ndim != 1 or d.ndim != 1 or n.shape != d.shape:
        raise ValueError("grade_n and grade_d must be 1-D arrays of equal length")
    if len(n) == 0:
        raise ValueError("pluto_tasche requires at least one grade")
    confidence = validate_level(confidence, "confidence")
    if np.any(n < 0.0) or np.any(d < 0.0):
        raise ValueError("grade_n and grade_d must be non-negative")
    if np.any(d > n):
        raise ValueError("grade_d cannot exceed grade_n in any grade")

    k = len(n)
    if grades is None:
        grade_labels = tuple(str(i + 1) for i in range(k))
    else:
        grade_labels = tuple(str(g) for g in np.asarray(grades).reshape(-1))
        if len(grade_labels) != k:
            raise ValueError("grades must have the same length as grade_n/grade_d")

    n_pooled = np.cumsum(n[::-1])[::-1]
    d_pooled = np.cumsum(d[::-1])[::-1]
    if np.any(n_pooled == 0.0):
        raise ValueError("pluto_tasche: grade pooled with worse grades has zero obligors (n* == 0)")

    # d <= n per grade implies d* <= n*, so the second shape is > 0 unless d* == n*.
    pd_upper = np.array(
        [
            1.0 if ds == ns else beta_ppf(confidence, ds + 1.0, ns - ds)
            for ns, ds in zip(n_pooled, d_pooled, strict=True)
        ]
    )
    # Most-prudent touch-up of a local dip: raise the worse grade, never lower the better.
    hull = np.maximum.accumulate(pd_upper)
    monotonized = bool(np.any(hull != pd_upper))

    return PlutoTascheResult(
        grades=grade_labels,
        n=n,
        d=d,
        n_pooled=n_pooled,
        d_pooled=d_pooled,
        pd_upper=hull,
        confidence=confidence,
        monotonized=monotonized,
    )

pluto_tasche_from_arrays ¶

pluto_tasche_from_arrays(grades: object, y: object, *, order: object = None, p: object = None, confidence: float = 0.9, sample_weight: object = None) -> PlutoTascheResult

Pluto-Tasche most-prudent PD from observation-level grades and outcomes.

Convenience wrapper around :func:pluto_tasche: aggregates y by grades into per-grade obligor/default counts (weighted sums when sample_weight is given) in the explicit order, then applies the same pooling and bound.

PARAMETER	DESCRIPTION
`grades`	Rating grade label per observation, or a :class:`probcal.Masterscale` that assigns them from `p`. TYPE: `array_like or Masterscale`
`y`	Binary outcomes in `{0, 1}`. Unlike most probcal metrics, an all-zero `y` is accepted -- Pluto-Tasche is built for exactly that case. TYPE: `array_like`
`order`	Explicit best-to-worst grade order; must match the unique labels in `grades` exactly (same set, same count, any order raises if it does not correspond to a permutation of the unique labels). Required for a label array; defaults to the masterscale's own order otherwise. TYPE: `sequence of str or None, keyword-only` DEFAULT: `None`
`p`	Predicted probabilities, required when `grades` is a `Masterscale` (labels are assigned from it); ignored otherwise. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`confidence`	Confidence level in `(0, 1)` for every grade's upper bound. TYPE: `(float, keyword - only)` DEFAULT: `0.9`
`sample_weight`	Optional positive weights, same length as `y`; per-grade counts become weighted sums, i.e. weights are read literally as frequency counts (integer weights equal row duplication) -- an exception to the relative-weight convention of :mod:`probcal.metrics._common`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`PlutoTascheResult`	Per-grade counts, pooled counts, and most-prudent PDs.

RAISES	DESCRIPTION
`ValueError`	If `grades` and `y` are not equal-length 1-D arrays, `y` contains values outside `{0, 1}`, or `order` does not match the unique labels in `grades`.

Examples:

>>> import numpy as np
>>> from probcal.metrics import pluto_tasche_from_arrays
>>> grades = np.array(["A"] * 100 + ["B"] * 400 + ["C"] * 300)
>>> y = np.zeros(800)
>>> res = pluto_tasche_from_arrays(grades, y, order=("A", "B", "C"))
>>> res.n
array([100., 400., 300.])

Source code in src/probcal/metrics/_conservative.py

def pluto_tasche_from_arrays(
    grades: object,
    y: object,
    *,
    order: object = None,
    p: object = None,
    confidence: float = 0.9,
    sample_weight: object = None,
) -> PlutoTascheResult:
    """Pluto-Tasche most-prudent PD from observation-level grades and outcomes.

    Convenience wrapper around :func:`pluto_tasche`: aggregates ``y`` by
    ``grades`` into per-grade obligor/default counts (weighted sums when
    ``sample_weight`` is given) in the explicit ``order``, then applies the
    same pooling and bound.

    Parameters
    ----------
    grades : array_like or Masterscale
        Rating grade label per observation, or a :class:`probcal.Masterscale`
        that assigns them from ``p``.
    y : array_like
        Binary outcomes in ``{0, 1}``. Unlike most probcal metrics, an
        all-zero ``y`` is accepted -- Pluto-Tasche is built for exactly that
        case.
    order : sequence of str or None, keyword-only
        Explicit best-to-worst grade order; must match the unique labels in
        ``grades`` exactly (same set, same count, any order raises if it
        does not correspond to a permutation of the unique labels). Required
        for a label array; defaults to the masterscale's own order otherwise.
    p : array_like or None, keyword-only
        Predicted probabilities, required when ``grades`` is a
        ``Masterscale`` (labels are assigned from it); ignored otherwise.
    confidence : float, keyword-only
        Confidence level in ``(0, 1)`` for every grade's upper bound.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``; per-grade counts
        become weighted sums, i.e. weights are read literally as frequency
        counts (integer weights equal row duplication) -- an exception to
        the relative-weight convention of :mod:`probcal.metrics._common`.

    Returns
    -------
    PlutoTascheResult
        Per-grade counts, pooled counts, and most-prudent PDs.

    Raises
    ------
    ValueError
        If ``grades`` and ``y`` are not equal-length 1-D arrays, ``y``
        contains values outside ``{0, 1}``, or ``order`` does not match the
        unique labels in ``grades``.

    Examples
    --------
    >>> import numpy as np
    >>> from probcal.metrics import pluto_tasche_from_arrays
    >>> grades = np.array(["A"] * 100 + ["B"] * 400 + ["C"] * 300)
    >>> y = np.zeros(800)
    >>> res = pluto_tasche_from_arrays(grades, y, order=("A", "B", "C"))
    >>> res.n
    array([100., 400., 300.])
    """
    y_arr = validate_binary_y(y, require_both_classes=False)
    if is_masterscale(grades):
        if p is None:
            raise ValueError(
                "p is required when grades is a Masterscale (labels are assigned from p)"
            )
        p_arr = validate_scores(p, name="p")
        if len(p_arr) != len(y_arr):
            raise ValueError("y and p must have equal length")
    else:
        if order is None:
            raise ValueError("order is required when grades is a label array")
        p_arr = np.empty(0)
    g_str, default_order = grade_labels(grades, p_arr, len(y_arr))
    w_arr = validate_weights(sample_weight, len(y_arr))
    order_t = resolve_order(g_str, order, default_order)
    n, d, _, _ = aggregate(g_str, order_t, y_arr, w=w_arr)

    return pluto_tasche(n, d, confidence=confidence, grades=order_t)

skce ¶

skce(y: object, p: object, *, estimator: str = 'uq', kernel: str = 'laplace', bandwidth: float | None = None, scale: str = 'probability', random_state: int = 42) -> float

Squared kernel calibration error (Widmann et al., 2019, Table 1).

"uq" (default) is the unbiased quadratic estimator (may be negative); "ul" the unbiased linear O(n) estimator over seeded disjoint pairs (random_state controls the pairing); "biased" the nonnegative V-statistic. bandwidth=None uses the deterministic median heuristic; scale="logit" transforms the kernel input only (the low-PD option).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`estimator`	Estimator variant; see above. TYPE: `(uq, ul, biased)` DEFAULT: `"uq"`
`kernel`	Kernel family applied to the (scaled) score distance. TYPE: `(laplace, gaussian)` DEFAULT: `"laplace"`
`bandwidth`	Kernel bandwidth; `None` (default) uses the median-heuristic distance (mean fallback if the median is 0). TYPE: `(float or None, keyword - only)` DEFAULT: `None`
`scale`	Scale on which the kernel input `s` is computed; residuals stay on the probability scale regardless. TYPE: `(probability, logit)` DEFAULT: `"probability"`
`random_state`	Seed for the `"ul"` estimator's disjoint-pairing permutation (unused by `"uq"`/`"biased"`). TYPE: `(int, keyword - only)` DEFAULT: `42`

RETURNS	DESCRIPTION
`float`	SKCE point estimate.

RAISES	DESCRIPTION
`ValueError`	If `estimator` is not one of `"uq"`, `"ul"`, `"biased"`, or if fewer than 2 observations are given.

Source code in src/probcal/metrics/kernel.py

def skce(
    y: object,
    p: object,
    *,
    estimator: str = "uq",
    kernel: str = "laplace",
    bandwidth: float | None = None,
    scale: str = "probability",
    random_state: int = 42,
) -> float:
    """Squared kernel calibration error (Widmann et al., 2019, Table 1).

    ``"uq"`` (default) is the unbiased quadratic estimator (may be negative);
    ``"ul"`` the unbiased linear O(n) estimator over seeded disjoint pairs
    (``random_state`` controls the pairing); ``"biased"`` the nonnegative
    V-statistic. ``bandwidth=None`` uses the deterministic median heuristic;
    ``scale="logit"`` transforms the kernel input only (the low-PD option).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    estimator : {"uq", "ul", "biased"}, keyword-only
        Estimator variant; see above.
    kernel : {"laplace", "gaussian"}, keyword-only
        Kernel family applied to the (scaled) score distance.
    bandwidth : float or None, keyword-only
        Kernel bandwidth; ``None`` (default) uses the median-heuristic
        distance (mean fallback if the median is 0).
    scale : {"probability", "logit"}, keyword-only
        Scale on which the kernel input ``s`` is computed; residuals stay on
        the probability scale regardless.
    random_state : int, keyword-only
        Seed for the ``"ul"`` estimator's disjoint-pairing permutation
        (unused by ``"uq"``/``"biased"``).

    Returns
    -------
    float
        SKCE point estimate.

    Raises
    ------
    ValueError
        If ``estimator`` is not one of ``"uq"``, ``"ul"``, ``"biased"``, or if
        fewer than 2 observations are given.
    """
    if estimator not in ("uq", "ul", "biased"):
        raise ValueError(f"estimator must be 'uq', 'ul', or 'biased', got {estimator!r}")
    y_arr, p_arr, s, bw = _setup(y, p, scale, bandwidth, 2, "skce")
    if estimator == "ul":
        return float(np.mean(_ul_terms(y_arr, p_arr, s, kernel, bw, random_state)))
    h = _h_full(y_arr, p_arr, s, kernel, bw)
    if estimator == "biased":
        return float(h.sum() / len(p_arr) ** 2)
    return _uq(h)

skce_test ¶

skce_test(y: object, p: object, *, method: str = 'bootstrap', n_boot: int = 999, kernel: str = 'laplace', bandwidth: float | None = None, scale: str = 'probability', random_state: int = 42) -> SkceTestResult

Calibration test on the SKCE (Widmann et al., 2019, Sec. 6 / App. G).

"bootstrap" (default): quadratic statistic with Arcones–Giné centered resampling; O(n_boot * n^2) — the more powerful choice. "asymptotic": linear statistic, normal approximation (Corollary G.3); O(n), preferred for n >~ 20 000, but a single random pairing can miss slope-type miscalibration that the bootstrap test rejects (the paper's documented power gap). p_value_bound is the distribution-free worst case.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`method`	Test method; see above. TYPE: `(bootstrap, asymptotic)` DEFAULT: `"bootstrap"`
`n_boot`	Bootstrap replicate count (`"bootstrap"` method only). TYPE: `(int, keyword - only)` DEFAULT: `999`
`kernel`	Kernel family applied to the (scaled) score distance. TYPE: `(laplace, gaussian)` DEFAULT: `"laplace"`
`bandwidth`	Kernel bandwidth; `None` (default) uses the median-heuristic distance (mean fallback if the median is 0). TYPE: `(float or None, keyword - only)` DEFAULT: `None`
`scale`	Scale on which the kernel input `s` is computed; residuals stay on the probability scale regardless. TYPE: `(probability, logit)` DEFAULT: `"probability"`
`random_state`	Seed for the resampling (`"bootstrap"`) or disjoint-pairing (`"asymptotic"`) randomness. TYPE: `(int, keyword - only)` DEFAULT: `42`

RETURNS	DESCRIPTION
`SkceTestResult`	Test statistic, method, p-value, and worst-case bound.

Source code in src/probcal/metrics/kernel.py

def skce_test(
    y: object,
    p: object,
    *,
    method: str = "bootstrap",
    n_boot: int = 999,
    kernel: str = "laplace",
    bandwidth: float | None = None,
    scale: str = "probability",
    random_state: int = 42,
) -> SkceTestResult:
    """Calibration test on the SKCE (Widmann et al., 2019, Sec. 6 / App. G).

    ``"bootstrap"`` (default): quadratic statistic with Arcones–Giné centered
    resampling; O(n_boot * n^2) — the more powerful choice. ``"asymptotic"``:
    linear statistic, normal approximation (Corollary G.3); O(n), preferred
    for n >~ 20 000, but a single random pairing can miss slope-type
    miscalibration that the bootstrap test rejects (the paper's documented
    power gap). ``p_value_bound`` is the distribution-free worst case.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    method : {"bootstrap", "asymptotic"}, keyword-only
        Test method; see above.
    n_boot : int, keyword-only
        Bootstrap replicate count (``"bootstrap"`` method only).
    kernel : {"laplace", "gaussian"}, keyword-only
        Kernel family applied to the (scaled) score distance.
    bandwidth : float or None, keyword-only
        Kernel bandwidth; ``None`` (default) uses the median-heuristic
        distance (mean fallback if the median is 0).
    scale : {"probability", "logit"}, keyword-only
        Scale on which the kernel input ``s`` is computed; residuals stay on
        the probability scale regardless.
    random_state : int, keyword-only
        Seed for the resampling (``"bootstrap"``) or disjoint-pairing
        (``"asymptotic"``) randomness.

    Returns
    -------
    SkceTestResult
        Test statistic, method, p-value, and worst-case bound.
    """
    if method not in ("bootstrap", "asymptotic"):
        raise ValueError(f"method must be 'bootstrap' or 'asymptotic', got {method!r}")
    n_boot = validate_positive_int(n_boot, "n_boot")
    y_arr, p_arr, s, bw = _setup(y, p, scale, bandwidth, 4, "skce_test")
    n = len(p_arr)

    if method == "asymptotic":
        terms = _ul_terms(y_arr, p_arr, s, kernel, bw, random_state)
        stat = float(np.mean(terms))
        sd = float(np.std(terms, ddof=1))
        if sd == 0.0:
            p_value = 1.0 if stat <= 0.0 else 0.0
        else:
            z = math.sqrt(len(terms)) * stat / sd
            p_value = norm_sf(z)
        return SkceTestResult(
            statistic=stat,
            estimator="ul",
            method="asymptotic",
            p_value=p_value,
            p_value_bound=_p_value_bound(stat, n),
            bandwidth=bw,
            n_boot=None,
        )

    h = _h_full(y_arr, p_arr, s, kernel, bw)
    stat = _uq(h)
    t_obs = n * stat
    c = h.mean(axis=1)
    h_tilde = h - c[:, None] - c[None, :] + h.mean()
    rng = np.random.default_rng(random_state)
    counts = rng.multinomial(n, np.full(n, 1.0 / n), size=n_boot).astype(np.float64)
    quad = np.einsum("bi,ij,bj->b", counts, h_tilde, counts)
    # Same scale as t_obs = n * U_n = S / (n - 1), with S the off-diagonal sum:
    # the resampled off-diagonal sum over (n - 1). Dividing by n instead would
    # shrink every replicate by (n - 1) / n and make p-values anti-conservative.
    t_b = (quad - counts @ np.diag(h_tilde)) / (n - 1)
    p_value = float((1 + int(np.sum(t_b >= t_obs))) / (n_boot + 1))
    return SkceTestResult(
        statistic=stat,
        estimator="uq",
        method="bootstrap",
        p_value=p_value,
        p_value_bound=_p_value_bound(stat, n),
        bandwidth=bw,
        n_boot=n_boot,
    )

calibration_guardrails ¶

calibration_guardrails(y: object, p: object, *, sample_weight: object = None) -> GuardrailReport

Evaluate the three guardrail flags.

Printed in selection reports and offset audit reports.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`GuardrailReport`	Slope, intercept, and Spiegelhalter-p values with pass/fail flags.

Source code in src/probcal/metrics/regression.py

def calibration_guardrails(
    y: object, p: object, *, sample_weight: object = None
) -> GuardrailReport:
    """Evaluate the three guardrail flags.

    Printed in selection reports and offset audit reports.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    GuardrailReport
        Slope, intercept, and Spiegelhalter-p values with pass/fail flags.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    slope = float(_slope_fit(y_arr, p_arr, w, "calibration_guardrails")[1][1])
    intercept = _intercept(y_arr, p_arr, w)
    sp = _spiegelhalter(y_arr, p_arr, w)
    slope_ok = 0.9 <= slope <= 1.1
    intercept_ok = abs(intercept) <= 0.1
    sp_ok = sp.p_value > 0.05
    return GuardrailReport(
        slope=slope,
        intercept=intercept,
        spiegelhalter_p=sp.p_value,
        slope_ok=slope_ok,
        intercept_ok=intercept_ok,
        spiegelhalter_ok=sp_ok,
        all_ok=slope_ok and intercept_ok and sp_ok,
    )

calibration_intercept ¶

calibration_intercept(y: object, p: object, *, sample_weight: object = None) -> float

Calibration-in-the-large in log-odds.

Logistic intercept with the slope fixed at 1 (offset regression on logit(p)).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Fitted intercept in log-odds units.

Source code in src/probcal/metrics/regression.py

def calibration_intercept(y: object, p: object, *, sample_weight: object = None) -> float:
    """Calibration-in-the-large in log-odds.

    Logistic intercept with the slope fixed at 1 (offset regression on
    logit(p)).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Fitted intercept in log-odds units.
    """
    return _intercept(*_prep(y, p, sample_weight))

calibration_slope ¶

calibration_slope(y: object, p: object, *, sample_weight: object = None) -> float

Cox calibration slope.

< 1 means overfitting/overconfident spread, > 1 underfitting. When every prediction has the same logit(p) (a constant p, or a single PD level) the slope is not identified and nan is returned with a UserWarning.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Fitted slope on the logit scale (`nan` when `logit(p)` is constant).

Source code in src/probcal/metrics/regression.py

def calibration_slope(y: object, p: object, *, sample_weight: object = None) -> float:
    """Cox calibration slope.

    ``< 1`` means overfitting/overconfident spread, ``> 1`` underfitting.
    When every prediction has the same ``logit(p)`` (a constant ``p``, or a
    single PD level) the slope is not identified and ``nan`` is returned
    with a ``UserWarning``.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Fitted slope on the logit scale (``nan`` when ``logit(p)`` is
        constant).
    """
    return float(_slope_fit(*_prep(y, p, sample_weight))[1][1])

calibration_test ¶

calibration_test(y: object, p: object, *, sample_weight: object = None) -> CalibrationTestResult

Likelihood-ratio test of joint calibration (alpha, beta) = (0, 1).

With weights, the fit and the likelihood ratio use sample_weight as given, and the weighted LR statistic -- distributed under calibration as a weighted sum of chi2_1 variables, not as chi2_2 -- is referred to a second-order Rao-Scott corrected chi-square (:func:probcal.metrics._common.rao_scott; the design effects use the exact null variance sum(w**2 p (1 - p) x x') of the weighted score). The test is invariant to rescaling the weights, and unit (or constant) weights give the plain 2-df likelihood-ratio test.

When every prediction has the same logit(p) the slope is not identified: only calibration-in-the-large is tested (slope fixed at 1, df = 1), beta is nan, and a UserWarning is emitted.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`CalibrationTestResult`	Test statistic, p-value, and fitted intercept/slope.

Source code in src/probcal/metrics/regression.py

def calibration_test(
    y: object, p: object, *, sample_weight: object = None
) -> CalibrationTestResult:
    """Likelihood-ratio test of joint calibration (alpha, beta) = (0, 1).

    With weights, the fit and the likelihood ratio use ``sample_weight`` as
    given, and the weighted LR statistic -- distributed under calibration as
    a weighted sum of ``chi2_1`` variables, not as ``chi2_2`` -- is referred
    to a second-order Rao-Scott corrected chi-square
    (:func:`probcal.metrics._common.rao_scott`; the design effects use the
    exact null variance ``sum(w**2 p (1 - p) x x')`` of the weighted score).
    The test is invariant to rescaling the weights, and unit (or constant)
    weights give the plain 2-df likelihood-ratio test.

    When every prediction has the same ``logit(p)`` the slope is not
    identified: only calibration-in-the-large is tested (slope fixed at 1,
    ``df = 1``), ``beta`` is ``nan``, and a ``UserWarning`` is emitted.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    CalibrationTestResult
        Test statistic, p-value, and fitted intercept/slope.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    X, beta_hat = _slope_fit(y_arr, p_arr, w, "calibration_test")
    alpha, beta = float(beta_hat[0]), float(beta_hat[1])

    if math.isnan(beta):  # constant logit(p): intercept-only (offset) test
        X = X[:, :1]
        fitted = expit(alpha + logit(p_arr))
    else:
        fitted = expit(X @ beta_hat)
    lr = max(2.0 * (bern_loglik(y_arr, fitted, w) - bern_loglik(y_arr, p_arr, w)), 0.0)
    stat, df = rao_scott(lr, X, p_arr, w)
    p_value = chi2_sf(stat, df)
    return CalibrationTestResult(statistic=stat, p_value=p_value, alpha=alpha, beta=beta, df=df)

brier_score ¶

brier_score(y: object, p: object, *, sample_weight: object = None) -> float

Weighted mean squared error of the probability forecast (strictly proper).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Weighted mean squared error.

Source code in src/probcal/metrics/scores.py

def brier_score(y: object, p: object, *, sample_weight: object = None) -> float:
    """Weighted mean squared error of the probability forecast (strictly proper).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Weighted mean squared error.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    return float(np.average((p_arr - y_arr) ** 2, weights=w))

brier_skill_score ¶

brier_skill_score(y: object, p: object, *, sample_weight: object = None) -> float

Brier skill score vs the climatology forecast p = mean(y).

Positive values beat the base rate; 0 equals it.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Skill score relative to the weighted base rate.

Source code in src/probcal/metrics/scores.py

def brier_skill_score(y: object, p: object, *, sample_weight: object = None) -> float:
    """Brier skill score vs the climatology forecast ``p = mean(y)``.

    Positive values beat the base rate; 0 equals it.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Skill score relative to the weighted base rate.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    base = float(np.average(y_arr, weights=w))
    bs_ref = float(np.average((base - y_arr) ** 2, weights=w))
    bs = float(np.average((p_arr - y_arr) ** 2, weights=w))
    return 1.0 - bs / bs_ref

log_loss ¶

log_loss(y: object, p: object, *, sample_weight: object = None) -> float

Weighted mean negative log-likelihood.

Strictly proper; the default selection criterion.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Weighted mean negative log-likelihood.

Source code in src/probcal/metrics/scores.py

def log_loss(y: object, p: object, *, sample_weight: object = None) -> float:
    """Weighted mean negative log-likelihood.

    Strictly proper; the default selection criterion.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Weighted mean negative log-likelihood.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    ll = y_arr * np.log(p_arr) + (1.0 - y_arr) * np.log1p(-p_arr)
    return float(-np.average(ll, weights=w))

logloss_calibration_refinement ¶

logloss_calibration_refinement(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None) -> LogLossDecomposition

Split the log loss into a calibration (KL) and refinement (entropy) part.

The conditional event rate c(p) is estimated by a LOESS smoother of the outcome on the prediction; calibration is the mean KL(Bernoulli(c) || Bernoulli(p)) and refinement the mean entropy of Bernoulli(c). Only as good as the plug-in estimate of c.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`frac`	LOESS smoothing fraction passed through to the recalibration curve. TYPE: `(float, keyword - only)` DEFAULT: `0.75`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`LogLossDecomposition`	Calibration and refinement terms.

Source code in src/probcal/metrics/scores.py

def logloss_calibration_refinement(
    y: object,
    p: object,
    *,
    frac: float = 0.75,
    sample_weight: object = None,
) -> LogLossDecomposition:
    """Split the log loss into a calibration (KL) and refinement (entropy) part.

    The conditional event rate ``c(p)`` is estimated by a LOESS smoother of
    the outcome on the prediction; calibration is the mean
    ``KL(Bernoulli(c) || Bernoulli(p))`` and refinement the mean entropy of
    ``Bernoulli(c)``. Only as good as the plug-in estimate of ``c``.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    frac : float, keyword-only
        LOESS smoothing fraction passed through to the recalibration curve.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    LogLossDecomposition
        Calibration and refinement terms.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    c = np.clip(loess(p_arr, y_arr, frac=frac), EPS, 1.0 - EPS)
    kl = c * (np.log(c) - np.log(p_arr)) + (1.0 - c) * (np.log1p(-c) - np.log1p(-p_arr))
    ent = -(c * np.log(c) + (1.0 - c) * np.log1p(-c))
    return LogLossDecomposition(
        calibration=float(np.average(kl, weights=w)),
        refinement=float(np.average(ent, weights=w)),
    )

murphy_curve ¶

murphy_curve(y: object, p: object, *, thresholds: object = 513, sample_weight: object = None) -> MurphyCurve

Murphy diagram data: mean elementary score of the mean functional across a threshold grid.

Uses the Ehm, Gneiting, Jordan & Krüger (2016) elementary score for the mean functional of a binary outcome,

S_theta(p, y) = theta * 1{p > theta, y = 0} + (1 - theta) * 1{p <= theta, y = 1},

whose weighted mean at each theta is this curve's score. 2 * integral(S_theta, theta in [0, 1]) equals the Brier score. S_theta is linear in theta between consecutive breakpoints u = sorted(unique(p) | {0, 1}) and jumps only at them, so the midpoint rule over u reproduces the Brier score to machine precision, while sampling at the breakpoints themselves converges only as n grows. Isotonic (PAV) recalibration never increases the score at any threshold (Ehm et al., 2016), so comparing the two curves diagnoses the value of recalibration. Cost O(n log n + T log n) for T thresholds.

PARAMETER	DESCRIPTION
`y`	Outcomes and predicted probabilities. TYPE: `array_like`
`p`	Outcomes and predicted probabilities. TYPE: `array_like`
`thresholds`	Either the number of equally spaced points in `[0, 1]` (`numpy.linspace(0, 1, thresholds)`; default 513, the same default resolution used elsewhere in the package for dense grids) or an explicit 1-D array of thresholds in `[0, 1]` (sorted internally). TYPE: `(int or array_like, keyword - only)` DEFAULT: `513`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`MurphyCurve`	Threshold grid, weighted mean elementary score, and observation count.

RAISES	DESCRIPTION
`ValueError`	If `thresholds` is not 1-D, or contains values outside `[0, 1]`.

Examples:

>>> import numpy as np
>>> from probcal.metrics import murphy_curve
>>> rng = np.random.default_rng(0)
>>> p = rng.uniform(0.05, 0.5, 300)
>>> y = (rng.random(300) < p).astype(float)
>>> curve = murphy_curve(y, p, thresholds=101)
>>> curve.score.shape
(101,)

Source code in src/probcal/metrics/scores.py

def murphy_curve(
    y: object,
    p: object,
    *,
    thresholds: object = 513,
    sample_weight: object = None,
) -> MurphyCurve:
    """Murphy diagram data: mean elementary score of the mean functional across a threshold grid.

    Uses the Ehm, Gneiting, Jordan & Krüger (2016) elementary score for the
    mean functional of a binary outcome,

    ``S_theta(p, y) = theta * 1{p > theta, y = 0} + (1 - theta) * 1{p <= theta, y = 1}``,

    whose weighted mean at each ``theta`` is this curve's ``score``.
    ``2 * integral(S_theta, theta in [0, 1])`` equals the Brier score.
    ``S_theta`` is linear in ``theta`` between consecutive breakpoints
    ``u = sorted(unique(p) | {0, 1})`` and jumps only at them, so the
    midpoint rule over ``u`` reproduces the Brier score to machine
    precision, while sampling at the breakpoints themselves converges only
    as ``n`` grows. Isotonic (PAV) recalibration never increases the score
    at any threshold (Ehm et al., 2016), so comparing the two curves
    diagnoses the value of recalibration. Cost ``O(n log n + T log n)`` for
    ``T`` thresholds.

    Parameters
    ----------
    y, p : array_like
        Outcomes and predicted probabilities.
    thresholds : int or array_like, keyword-only
        Either the number of equally spaced points in ``[0, 1]``
        (``numpy.linspace(0, 1, thresholds)``; default 513, the same
        default resolution used elsewhere in the package for dense grids)
        or an explicit 1-D array of thresholds in ``[0, 1]`` (sorted
        internally).
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    MurphyCurve
        Threshold grid, weighted mean elementary score, and observation
        count.

    Raises
    ------
    ValueError
        If ``thresholds`` is not 1-D, or contains values outside ``[0, 1]``.

    Examples
    --------
    >>> import numpy as np
    >>> from probcal.metrics import murphy_curve
    >>> rng = np.random.default_rng(0)
    >>> p = rng.uniform(0.05, 0.5, 300)
    >>> y = (rng.random(300) < p).astype(float)
    >>> curve = murphy_curve(y, p, thresholds=101)
    >>> curve.score.shape
    (101,)
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    if isinstance(thresholds, (int, np.integer)):
        theta = np.linspace(0.0, 1.0, int(thresholds))
    else:
        theta = np.sort(np.asarray(thresholds, dtype=np.float64))
        if theta.ndim != 1:
            raise ValueError(f"thresholds must be a 1-D array, got shape {theta.shape}")
        if np.any(theta < 0.0) or np.any(theta > 1.0):
            raise ValueError("thresholds must lie in [0, 1]")

    # Sort p once; cumulative weighted sums over the y==0/y==1 subsets let
    # searchsorted read off, for every threshold at once, how much weight
    # lies on each side — avoids an O(n * T) mask per threshold.
    order = np.argsort(p_arr, kind="stable")
    p_sorted = p_arr[order]
    is_event = y_arr[order] == 1.0
    w_sorted = w[order]
    cum_w0 = np.concatenate(([0.0], np.cumsum(np.where(is_event, 0.0, w_sorted))))
    cum_w1 = np.concatenate(([0.0], np.cumsum(np.where(is_event, w_sorted, 0.0))))

    idx = np.searchsorted(p_sorted, theta, side="right")
    w0_above = cum_w0[-1] - cum_w0[idx]  # y=0, p > theta
    w1_at_or_below = cum_w1[idx]  # y=1, p <= theta

    score = (theta * w0_above + (1.0 - theta) * w1_at_or_below) / w.sum()
    return MurphyCurve(thresholds=theta, score=score, n=len(y_arr))

murphy_decomposition ¶

murphy_decomposition(y: object, p: object, *, n_bins: int = 10, strategy: str = 'mass', bias_corrected: bool = False, sample_weight: object = None) -> MurphyDecomposition

Binned reliability/resolution/uncertainty split of the Brier score.

bias_corrected=True subtracts the within-bin variance of the event rate from the squared-gap terms (within-bin variance corrections in the manner of Ferro & Fricker, 2012); the naive plug-in otherwise. The decomposition inherits the binning choice — see the metrics chapter.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`n_bins`	Requested number of bins. TYPE: `(int, keyword - only)` DEFAULT: `10`
`strategy`	`"mass"` (equal-count, default) or `"width"` (equal-width over [0, 1]). TYPE: `(mass, width)` DEFAULT: `"mass"`
`bias_corrected`	If `True` (default `False`), apply the Ferro & Fricker (2012) within-bin variance correction to the reliability and resolution terms. The variance uses each bin's effective count `sum_b(w)2 / sum_b(w2)` (see :mod:`probcal.metrics._common`; the raw count for unit weights). TYPE: `(bool, keyword - only)` DEFAULT: `False`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`MurphyDecomposition`	Reliability, resolution, and uncertainty terms.

Source code in src/probcal/metrics/scores.py

def murphy_decomposition(
    y: object,
    p: object,
    *,
    n_bins: int = 10,
    strategy: str = "mass",
    bias_corrected: bool = False,
    sample_weight: object = None,
) -> MurphyDecomposition:
    """Binned reliability/resolution/uncertainty split of the Brier score.

    ``bias_corrected=True`` subtracts the within-bin variance of the event
    rate from the squared-gap terms (within-bin variance corrections in the
    manner of Ferro & Fricker, 2012); the naive plug-in otherwise. The
    decomposition inherits the binning choice — see the metrics chapter.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    n_bins : int, keyword-only
        Requested number of bins.
    strategy : {"mass", "width"}, keyword-only
        ``"mass"`` (equal-count, default) or ``"width"`` (equal-width over [0, 1]).
    bias_corrected : bool, keyword-only
        If ``True`` (default ``False``), apply the Ferro & Fricker (2012)
        within-bin variance correction to the reliability and resolution terms.
        The variance uses each bin's effective count ``sum_b(w)**2 /
        sum_b(w**2)`` (see :mod:`probcal.metrics._common`; the raw count for
        unit weights).
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    MurphyDecomposition
        Reliability, resolution, and uncertainty terms.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    bs = bin_stats(y_arr, p_arr, w, n_bins, strategy)
    share = bs.w_sum / float(w.sum())
    y_bar = float(np.average(y_arr, weights=w))
    yb = bs.rate
    rel_term = (bs.p_mean - yb) ** 2
    res_term = (yb - y_bar) ** 2
    if bias_corrected:
        rel_term, res_term = bs.debias(rel_term), bs.debias(res_term)
    rel = float(np.sum(share * rel_term))
    res = float(np.sum(share * res_term))
    unc = y_bar * (1.0 - y_bar)
    return MurphyDecomposition(reliability=rel, resolution=res, uncertainty=unc)

e50 ¶

e50(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None, grid_size: int | None = 512) -> float

Median of the |LOESS(y|p) - p| distances.

grid_size=None fits the LOESS exactly at every prediction. The LOESS distances are always unweighted; sample_weight, when given and not uniform, weights only the quantile step (see :func:~probcal._math.weighted_quantile).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`frac`	LOESS smoothing fraction. TYPE: `(float, keyword - only)` DEFAULT: `0.75`
`sample_weight`	Optional positive weights; used only for the quantile step. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`grid_size`	LOESS evaluation grid size; `None` fits exactly at every prediction. TYPE: `(int or None, keyword - only)` DEFAULT: `512`

RETURNS	DESCRIPTION
`float`	Median of the LOESS distances.

Source code in src/probcal/metrics/smooth.py

def e50(
    y: object,
    p: object,
    *,
    frac: float = 0.75,
    sample_weight: object = None,
    grid_size: int | None = 512,
) -> float:
    """Median of the |LOESS(y|p) - p| distances.

    ``grid_size=None`` fits the LOESS exactly at every prediction. The LOESS distances are
    always unweighted; ``sample_weight``, when given and
    not uniform, weights only the quantile step (see
    :func:`~probcal._math.weighted_quantile`).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    frac : float, keyword-only
        LOESS smoothing fraction.
    sample_weight : array_like or None, keyword-only
        Optional positive weights; used only for the quantile step.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` fits exactly at every prediction.

    Returns
    -------
    float
        Median of the LOESS distances.
    """
    return _ici_public("e50", y, p, frac, sample_weight, grid_size)

e90 ¶

e90(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None, grid_size: int | None = 512) -> float

90th percentile of the |LOESS(y|p) - p| distances.

grid_size=None fits the LOESS exactly at every prediction. The LOESS distances are always unweighted; sample_weight, when given and not uniform, weights only the quantile step (see :func:~probcal._math.weighted_quantile).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`frac`	LOESS smoothing fraction. TYPE: `(float, keyword - only)` DEFAULT: `0.75`
`sample_weight`	Optional positive weights; used only for the quantile step. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`grid_size`	LOESS evaluation grid size; `None` fits exactly at every prediction. TYPE: `(int or None, keyword - only)` DEFAULT: `512`

RETURNS	DESCRIPTION
`float`	90th percentile of the LOESS distances.

Source code in src/probcal/metrics/smooth.py

def e90(
    y: object,
    p: object,
    *,
    frac: float = 0.75,
    sample_weight: object = None,
    grid_size: int | None = 512,
) -> float:
    """90th percentile of the |LOESS(y|p) - p| distances.

    ``grid_size=None`` fits the LOESS exactly at every prediction. The LOESS distances are
    always unweighted; ``sample_weight``, when given and
    not uniform, weights only the quantile step (see
    :func:`~probcal._math.weighted_quantile`).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    frac : float, keyword-only
        LOESS smoothing fraction.
    sample_weight : array_like or None, keyword-only
        Optional positive weights; used only for the quantile step.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` fits exactly at every prediction.

    Returns
    -------
    float
        90th percentile of the LOESS distances.
    """
    return _ici_public("e90", y, p, frac, sample_weight, grid_size)

ecce ¶

ecce(y: object, p: object, *, sample_weight: object = None) -> EcceResult

Cumulative-deviation calibration error (Arrieta-Ibarra et al., 2022).

Sort by prediction and walk the cumulative sum of weighted residuals, normalized by the total weight; under calibration the walk hovers near zero, and drift localizes miscalibration without any smoothing parameter. The walk is read only at the end of each block of tied predictions, so the result does not depend on the input order of tied rows; stat_mean weights each block by its total sample weight. Without ties and with unit weights the walk is read after every row.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`EcceResult`	Max and mean absolute cumulative deviation.

Source code in src/probcal/metrics/smooth.py

def ecce(y: object, p: object, *, sample_weight: object = None) -> EcceResult:
    """Cumulative-deviation calibration error (Arrieta-Ibarra et al., 2022).

    Sort by prediction and walk the cumulative sum of weighted residuals,
    normalized by the total weight; under calibration the walk hovers near
    zero, and drift localizes miscalibration without any smoothing
    parameter. The walk is read only at the end of each block of tied
    predictions, so the result does not depend on the input order of tied
    rows; ``stat_mean`` weights each block by its total sample weight.
    Without ties and with unit weights the walk is read after every row.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    EcceResult
        Max and mean absolute cumulative deviation.
    """
    return _ecce_stats(_ecce_walk(*_prep(y, p, sample_weight)))

emax ¶

emax(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None, grid_size: int | None = 512) -> float

Maximum of the |LOESS(y|p) - p| distances.

grid_size=None fits the LOESS exactly at every prediction.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`frac`	LOESS smoothing fraction. TYPE: `(float, keyword - only)` DEFAULT: `0.75`
`sample_weight`	Optional positive weights, validated like everywhere else. The weighted maximum over strictly positive weights is the plain maximum, so the value does not depend on them. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`grid_size`	LOESS evaluation grid size; `None` fits exactly at every prediction. TYPE: `(int or None, keyword - only)` DEFAULT: `512`

RETURNS	DESCRIPTION
`float`	Maximum of the LOESS distances.

Source code in src/probcal/metrics/smooth.py

def emax(
    y: object,
    p: object,
    *,
    frac: float = 0.75,
    sample_weight: object = None,
    grid_size: int | None = 512,
) -> float:
    """Maximum of the |LOESS(y|p) - p| distances.

    ``grid_size=None`` fits the LOESS exactly at every prediction.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    frac : float, keyword-only
        LOESS smoothing fraction.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, validated like everywhere else. The
        weighted maximum over strictly positive weights is the plain maximum,
        so the value does not depend on them.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` fits exactly at every prediction.

    Returns
    -------
    float
        Maximum of the LOESS distances.
    """
    return _ici_public("emax", y, p, frac, sample_weight, grid_size)

ici ¶

ici(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None, grid_size: int | None = 512) -> float

Integrated calibration index: weighted mean |LOESS(y|p) - p| (Austin & Steyerberg, 2019).

The LOESS stage itself is unweighted. grid_size=None fits the LOESS exactly at every prediction.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`frac`	LOESS smoothing fraction. TYPE: `(float, keyword - only)` DEFAULT: `0.75`
`sample_weight`	Optional positive weights, same length as `y`; weights only the final averaging step, not the LOESS fit. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`grid_size`	LOESS evaluation grid size; `None` fits exactly at every prediction. TYPE: `(int or None, keyword - only)` DEFAULT: `512`

RETURNS	DESCRIPTION
`float`	Weighted mean absolute LOESS-to-prediction distance.

Source code in src/probcal/metrics/smooth.py

def ici(
    y: object,
    p: object,
    *,
    frac: float = 0.75,
    sample_weight: object = None,
    grid_size: int | None = 512,
) -> float:
    """Integrated calibration index: weighted mean |LOESS(y|p) - p|
    (Austin & Steyerberg, 2019).

    The LOESS stage itself is unweighted.
    ``grid_size=None`` fits the LOESS exactly at every prediction.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    frac : float, keyword-only
        LOESS smoothing fraction.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``; weights only the
        final averaging step, not the LOESS fit.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` fits exactly at every prediction.

    Returns
    -------
    float
        Weighted mean absolute LOESS-to-prediction distance.
    """
    return _ici_public("ici", y, p, frac, sample_weight, grid_size)

smooth_ece ¶

smooth_ece(y: object, p: object, *, sample_weight: object = None, bins: int | None = 8192) -> float

Kernel-smoothed ECE with a self-consistent bandwidth (Błasiok–Nakkiran).

Residuals are smoothed with a Gaussian kernel on the logit scale (the paper's reflected kernel is a boundary device for [0, 1]; on the unbounded logit scale no reflection is needed), and the reported value is the fixed point smECE(sigma) = sigma found by bisection.

bins aggregates the weighted residual measure onto an equal-width logit lattice and solves the fixed point there by direct convolution (cost independent of n and sigma). If the solved sigma spans fewer than 8 cells, the solve is repeated once on a lattice refined to 8 cells per sigma; if that would exceed 2**20 cells or is still under-resolved -- or with bins=None or a degenerate logit range -- the exact computation (a 257-point grid evaluator, O(n) per bisection step) runs instead. The lattice value differs from the exact one by ~1e-4 on typical portfolios; on wide clipped-logit data the exact path's fixed grid under-resolves small-sigma kernels and the lattice value is the better one.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`bins`	Number of lattice bins for the fast path (default 8192); `None` forces the exact O(n) computation. TYPE: `(int or None, keyword - only)` DEFAULT: `8192`

RETURNS	DESCRIPTION
`float`	smECE: the fixed point `sigma` solving `smECE(sigma) = sigma`.

Source code in src/probcal/metrics/smooth.py

def smooth_ece(
    y: object, p: object, *, sample_weight: object = None, bins: int | None = 8192
) -> float:
    """Kernel-smoothed ECE with a self-consistent bandwidth (Błasiok–Nakkiran).

    Residuals are smoothed with a Gaussian kernel on the logit scale (the
    paper's reflected kernel is a boundary device for [0, 1]; on the
    unbounded logit scale no reflection is needed), and the
    reported value is the fixed point ``smECE(sigma) = sigma`` found by
    bisection.

    ``bins`` aggregates the weighted residual measure onto an equal-width
    logit lattice and solves the fixed point there by direct convolution
    (cost independent of n and sigma). If the solved ``sigma`` spans fewer
    than 8 cells, the solve is repeated once on a lattice refined to 8 cells
    per sigma; if that would exceed ``2**20`` cells or is still
    under-resolved -- or with ``bins=None`` or a degenerate logit range --
    the exact computation (a 257-point grid evaluator, O(n) per bisection
    step) runs instead. The lattice value differs from the exact one by
    ~1e-4 on typical portfolios; on wide clipped-logit data the exact path's
    fixed grid under-resolves small-sigma kernels and the lattice value is
    the better one.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.
    bins : int or None, keyword-only
        Number of lattice bins for the fast path (default 8192); ``None``
        forces the exact O(n) computation.

    Returns
    -------
    float
        smECE: the fixed point ``sigma`` solving ``smECE(sigma) = sigma``.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    t = logit(p_arr)
    mass = (w / w.sum()) * (y_arr - p_arr)
    return _smece_solve(t, mass, bins)[0]

spiegelhalter_z ¶

spiegelhalter_z(y: object, p: object, *, sample_weight: object = None) -> SpiegelhalterResult

Spiegelhalter (1986) z statistic built on the Brier score.

The numerator has expectation zero under calibration; the statistic is asymptotically standard normal. No binning, no smoothing; aggregates the whole range, so compensating regional errors can cancel.

With weights the numerator is sum(w (y - p) (1 - 2p)) and the null variance its exact variance under calibration, sum(w**2 (1 - 2p)^2 p (1 - p)) (package convention, :mod:probcal.metrics._common): z is invariant to rescaling the weights and unit weights give the unweighted statistic. A Kish-rescaled effective sample size is not used because it misstates the variance when the weights vary with p (1 - p).

The null variance is zero when every prediction is exactly 0.5 (or clipped to 0/1); the statistic is then undefined and z and p_value are nan.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`SpiegelhalterResult`	Z statistic and two-sided p-value.

Source code in src/probcal/metrics/smooth.py

def spiegelhalter_z(y: object, p: object, *, sample_weight: object = None) -> SpiegelhalterResult:
    """Spiegelhalter (1986) z statistic built on the Brier score.

    The numerator has expectation zero under calibration; the statistic is
    asymptotically standard normal. No binning, no smoothing; aggregates the
    whole range, so compensating regional errors can cancel.

    With weights the numerator is ``sum(w (y - p) (1 - 2p))`` and the null
    variance its exact variance under calibration, ``sum(w**2 (1 - 2p)^2 p
    (1 - p))`` (package convention, :mod:`probcal.metrics._common`): ``z``
    is invariant to rescaling the weights and unit weights give the
    unweighted statistic. A Kish-rescaled effective sample size is not used
    because it misstates the variance when the weights vary with ``p (1 - p)``.

    The null variance is zero when every prediction is exactly 0.5 (or
    clipped to 0/1); the statistic is then undefined and ``z`` and
    ``p_value`` are ``nan``.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    SpiegelhalterResult
        Z statistic and two-sided p-value.
    """
    return _spiegelhalter(*_prep(y, p, sample_weight))

scores ¶

Proper scoring rules: log loss, Brier score, and their decompositions.

Theory, formulas, and pathologies: docs/concepts/metrics.md.

MurphyDecomposition `dataclass` ¶

MurphyDecomposition(reliability: float, resolution: float, uncertainty: float)

Binned Murphy (1973) partition of the Brier score.

reliability - resolution + uncertainty equals the Brier score exactly when predictions are constant within bins; otherwise the identity holds up to the within-bin variance of p (documented binning bias).

ATTRIBUTE	DESCRIPTION
`reliability`	Mean squared gap between within-bin predicted and observed rates. TYPE: `float`
`resolution`	Mean squared gap between within-bin observed rate and the overall base rate. TYPE: `float`
`uncertainty`	Base-rate variance `y_bar * (1 - y_bar)`. TYPE: `float`

MurphyCurve `dataclass` ¶

MurphyCurve(thresholds: ndarray, score: ndarray, n: int)

Bases: _ResultBase

Murphy diagram: mean elementary score of the binary mean functional across thresholds.

ATTRIBUTE	DESCRIPTION
`thresholds`	Threshold grid :math:`\theta \in [0, 1]`. TYPE: `ndarray`
`score`	Weighted mean elementary score :math:`S_\theta` at each threshold. TYPE: `ndarray`
`n`	Number of observations. TYPE: `int`

LogLossDecomposition `dataclass` ¶

LogLossDecomposition(calibration: float, refinement: float)

Calibration/refinement split of the log loss via a plug-in recalibration curve (LOESS).

ATTRIBUTE	DESCRIPTION
`calibration`	Mean KL divergence between the plug-in and predicted Bernoullis. TYPE: `float`
`refinement`	Mean entropy of the plug-in Bernoulli. TYPE: `float`

log_loss ¶

log_loss(y: object, p: object, *, sample_weight: object = None) -> float

Weighted mean negative log-likelihood.

Strictly proper; the default selection criterion.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Weighted mean negative log-likelihood.

Source code in src/probcal/metrics/scores.py

def log_loss(y: object, p: object, *, sample_weight: object = None) -> float:
    """Weighted mean negative log-likelihood.

    Strictly proper; the default selection criterion.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Weighted mean negative log-likelihood.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    ll = y_arr * np.log(p_arr) + (1.0 - y_arr) * np.log1p(-p_arr)
    return float(-np.average(ll, weights=w))

brier_score ¶

brier_score(y: object, p: object, *, sample_weight: object = None) -> float

Weighted mean squared error of the probability forecast (strictly proper).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Weighted mean squared error.

Source code in src/probcal/metrics/scores.py

def brier_score(y: object, p: object, *, sample_weight: object = None) -> float:
    """Weighted mean squared error of the probability forecast (strictly proper).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Weighted mean squared error.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    return float(np.average((p_arr - y_arr) ** 2, weights=w))

brier_skill_score ¶

brier_skill_score(y: object, p: object, *, sample_weight: object = None) -> float

Brier skill score vs the climatology forecast p = mean(y).

Positive values beat the base rate; 0 equals it.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Skill score relative to the weighted base rate.

Source code in src/probcal/metrics/scores.py

def brier_skill_score(y: object, p: object, *, sample_weight: object = None) -> float:
    """Brier skill score vs the climatology forecast ``p = mean(y)``.

    Positive values beat the base rate; 0 equals it.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Skill score relative to the weighted base rate.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    base = float(np.average(y_arr, weights=w))
    bs_ref = float(np.average((base - y_arr) ** 2, weights=w))
    bs = float(np.average((p_arr - y_arr) ** 2, weights=w))
    return 1.0 - bs / bs_ref

murphy_decomposition ¶

murphy_decomposition(y: object, p: object, *, n_bins: int = 10, strategy: str = 'mass', bias_corrected: bool = False, sample_weight: object = None) -> MurphyDecomposition

Binned reliability/resolution/uncertainty split of the Brier score.

bias_corrected=True subtracts the within-bin variance of the event rate from the squared-gap terms (within-bin variance corrections in the manner of Ferro & Fricker, 2012); the naive plug-in otherwise. The decomposition inherits the binning choice — see the metrics chapter.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`n_bins`	Requested number of bins. TYPE: `(int, keyword - only)` DEFAULT: `10`
`strategy`	`"mass"` (equal-count, default) or `"width"` (equal-width over [0, 1]). TYPE: `(mass, width)` DEFAULT: `"mass"`
`bias_corrected`	If `True` (default `False`), apply the Ferro & Fricker (2012) within-bin variance correction to the reliability and resolution terms. The variance uses each bin's effective count `sum_b(w)2 / sum_b(w2)` (see :mod:`probcal.metrics._common`; the raw count for unit weights). TYPE: `(bool, keyword - only)` DEFAULT: `False`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`MurphyDecomposition`	Reliability, resolution, and uncertainty terms.

Source code in src/probcal/metrics/scores.py

def murphy_decomposition(
    y: object,
    p: object,
    *,
    n_bins: int = 10,
    strategy: str = "mass",
    bias_corrected: bool = False,
    sample_weight: object = None,
) -> MurphyDecomposition:
    """Binned reliability/resolution/uncertainty split of the Brier score.

    ``bias_corrected=True`` subtracts the within-bin variance of the event
    rate from the squared-gap terms (within-bin variance corrections in the
    manner of Ferro & Fricker, 2012); the naive plug-in otherwise. The
    decomposition inherits the binning choice — see the metrics chapter.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    n_bins : int, keyword-only
        Requested number of bins.
    strategy : {"mass", "width"}, keyword-only
        ``"mass"`` (equal-count, default) or ``"width"`` (equal-width over [0, 1]).
    bias_corrected : bool, keyword-only
        If ``True`` (default ``False``), apply the Ferro & Fricker (2012)
        within-bin variance correction to the reliability and resolution terms.
        The variance uses each bin's effective count ``sum_b(w)**2 /
        sum_b(w**2)`` (see :mod:`probcal.metrics._common`; the raw count for
        unit weights).
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    MurphyDecomposition
        Reliability, resolution, and uncertainty terms.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    bs = bin_stats(y_arr, p_arr, w, n_bins, strategy)
    share = bs.w_sum / float(w.sum())
    y_bar = float(np.average(y_arr, weights=w))
    yb = bs.rate
    rel_term = (bs.p_mean - yb) ** 2
    res_term = (yb - y_bar) ** 2
    if bias_corrected:
        rel_term, res_term = bs.debias(rel_term), bs.debias(res_term)
    rel = float(np.sum(share * rel_term))
    res = float(np.sum(share * res_term))
    unc = y_bar * (1.0 - y_bar)
    return MurphyDecomposition(reliability=rel, resolution=res, uncertainty=unc)

murphy_curve ¶

murphy_curve(y: object, p: object, *, thresholds: object = 513, sample_weight: object = None) -> MurphyCurve

Murphy diagram data: mean elementary score of the mean functional across a threshold grid.

Uses the Ehm, Gneiting, Jordan & Krüger (2016) elementary score for the mean functional of a binary outcome,

S_theta(p, y) = theta * 1{p > theta, y = 0} + (1 - theta) * 1{p <= theta, y = 1},

whose weighted mean at each theta is this curve's score. 2 * integral(S_theta, theta in [0, 1]) equals the Brier score. S_theta is linear in theta between consecutive breakpoints u = sorted(unique(p) | {0, 1}) and jumps only at them, so the midpoint rule over u reproduces the Brier score to machine precision, while sampling at the breakpoints themselves converges only as n grows. Isotonic (PAV) recalibration never increases the score at any threshold (Ehm et al., 2016), so comparing the two curves diagnoses the value of recalibration. Cost O(n log n + T log n) for T thresholds.

PARAMETER	DESCRIPTION
`y`	Outcomes and predicted probabilities. TYPE: `array_like`
`p`	Outcomes and predicted probabilities. TYPE: `array_like`
`thresholds`	Either the number of equally spaced points in `[0, 1]` (`numpy.linspace(0, 1, thresholds)`; default 513, the same default resolution used elsewhere in the package for dense grids) or an explicit 1-D array of thresholds in `[0, 1]` (sorted internally). TYPE: `(int or array_like, keyword - only)` DEFAULT: `513`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`MurphyCurve`	Threshold grid, weighted mean elementary score, and observation count.

RAISES	DESCRIPTION
`ValueError`	If `thresholds` is not 1-D, or contains values outside `[0, 1]`.

Examples:

>>> import numpy as np
>>> from probcal.metrics import murphy_curve
>>> rng = np.random.default_rng(0)
>>> p = rng.uniform(0.05, 0.5, 300)
>>> y = (rng.random(300) < p).astype(float)
>>> curve = murphy_curve(y, p, thresholds=101)
>>> curve.score.shape
(101,)

Source code in src/probcal/metrics/scores.py

def murphy_curve(
    y: object,
    p: object,
    *,
    thresholds: object = 513,
    sample_weight: object = None,
) -> MurphyCurve:
    """Murphy diagram data: mean elementary score of the mean functional across a threshold grid.

    Uses the Ehm, Gneiting, Jordan & Krüger (2016) elementary score for the
    mean functional of a binary outcome,

    ``S_theta(p, y) = theta * 1{p > theta, y = 0} + (1 - theta) * 1{p <= theta, y = 1}``,

    whose weighted mean at each ``theta`` is this curve's ``score``.
    ``2 * integral(S_theta, theta in [0, 1])`` equals the Brier score.
    ``S_theta`` is linear in ``theta`` between consecutive breakpoints
    ``u = sorted(unique(p) | {0, 1})`` and jumps only at them, so the
    midpoint rule over ``u`` reproduces the Brier score to machine
    precision, while sampling at the breakpoints themselves converges only
    as ``n`` grows. Isotonic (PAV) recalibration never increases the score
    at any threshold (Ehm et al., 2016), so comparing the two curves
    diagnoses the value of recalibration. Cost ``O(n log n + T log n)`` for
    ``T`` thresholds.

    Parameters
    ----------
    y, p : array_like
        Outcomes and predicted probabilities.
    thresholds : int or array_like, keyword-only
        Either the number of equally spaced points in ``[0, 1]``
        (``numpy.linspace(0, 1, thresholds)``; default 513, the same
        default resolution used elsewhere in the package for dense grids)
        or an explicit 1-D array of thresholds in ``[0, 1]`` (sorted
        internally).
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    MurphyCurve
        Threshold grid, weighted mean elementary score, and observation
        count.

    Raises
    ------
    ValueError
        If ``thresholds`` is not 1-D, or contains values outside ``[0, 1]``.

    Examples
    --------
    >>> import numpy as np
    >>> from probcal.metrics import murphy_curve
    >>> rng = np.random.default_rng(0)
    >>> p = rng.uniform(0.05, 0.5, 300)
    >>> y = (rng.random(300) < p).astype(float)
    >>> curve = murphy_curve(y, p, thresholds=101)
    >>> curve.score.shape
    (101,)
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    if isinstance(thresholds, (int, np.integer)):
        theta = np.linspace(0.0, 1.0, int(thresholds))
    else:
        theta = np.sort(np.asarray(thresholds, dtype=np.float64))
        if theta.ndim != 1:
            raise ValueError(f"thresholds must be a 1-D array, got shape {theta.shape}")
        if np.any(theta < 0.0) or np.any(theta > 1.0):
            raise ValueError("thresholds must lie in [0, 1]")

    # Sort p once; cumulative weighted sums over the y==0/y==1 subsets let
    # searchsorted read off, for every threshold at once, how much weight
    # lies on each side — avoids an O(n * T) mask per threshold.
    order = np.argsort(p_arr, kind="stable")
    p_sorted = p_arr[order]
    is_event = y_arr[order] == 1.0
    w_sorted = w[order]
    cum_w0 = np.concatenate(([0.0], np.cumsum(np.where(is_event, 0.0, w_sorted))))
    cum_w1 = np.concatenate(([0.0], np.cumsum(np.where(is_event, w_sorted, 0.0))))

    idx = np.searchsorted(p_sorted, theta, side="right")
    w0_above = cum_w0[-1] - cum_w0[idx]  # y=0, p > theta
    w1_at_or_below = cum_w1[idx]  # y=1, p <= theta

    score = (theta * w0_above + (1.0 - theta) * w1_at_or_below) / w.sum()
    return MurphyCurve(thresholds=theta, score=score, n=len(y_arr))

logloss_calibration_refinement ¶

logloss_calibration_refinement(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None) -> LogLossDecomposition

Split the log loss into a calibration (KL) and refinement (entropy) part.

The conditional event rate c(p) is estimated by a LOESS smoother of the outcome on the prediction; calibration is the mean KL(Bernoulli(c) || Bernoulli(p)) and refinement the mean entropy of Bernoulli(c). Only as good as the plug-in estimate of c.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`frac`	LOESS smoothing fraction passed through to the recalibration curve. TYPE: `(float, keyword - only)` DEFAULT: `0.75`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`LogLossDecomposition`	Calibration and refinement terms.

Source code in src/probcal/metrics/scores.py

def logloss_calibration_refinement(
    y: object,
    p: object,
    *,
    frac: float = 0.75,
    sample_weight: object = None,
) -> LogLossDecomposition:
    """Split the log loss into a calibration (KL) and refinement (entropy) part.

    The conditional event rate ``c(p)`` is estimated by a LOESS smoother of
    the outcome on the prediction; calibration is the mean
    ``KL(Bernoulli(c) || Bernoulli(p))`` and refinement the mean entropy of
    ``Bernoulli(c)``. Only as good as the plug-in estimate of ``c``.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    frac : float, keyword-only
        LOESS smoothing fraction passed through to the recalibration curve.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    LogLossDecomposition
        Calibration and refinement terms.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    c = np.clip(loess(p_arr, y_arr, frac=frac), EPS, 1.0 - EPS)
    kl = c * (np.log(c) - np.log(p_arr)) + (1.0 - c) * (np.log1p(-c) - np.log1p(-p_arr))
    ent = -(c * np.log(c) + (1.0 - c) * np.log1p(-c))
    return LogLossDecomposition(
        calibration=float(np.average(kl, weights=w)),
        refinement=float(np.average(ent, weights=w)),
    )

binned ¶

Binned calibration-error estimators: ECE family, MCE, Hosmer–Lemeshow.

Pathologies (binning sensitivity, finite-sample bias, HL power issues) are documented in docs/concepts/metrics.md. None of these are selection criteria; the Hosmer–Lemeshow test is report-only. Sample weights follow the package convention in :mod:probcal.metrics._common (relative weights; per-bin variances from the weights themselves wherever a sample size enters).

HosmerLemeshowResult `dataclass` ¶

HosmerLemeshowResult(statistic: float, df: int, p_value: float)

Hosmer–Lemeshow chi-square test (report-only; never a selection criterion).

ATTRIBUTE	DESCRIPTION
`statistic`	Chi-square test statistic. TYPE: `float`
`df`	Degrees of freedom: non-empty groups minus `ddof` (`0` when no degree of freedom is left). TYPE: `int`
`p_value`	Upper-tail p-value of the statistic; `nan` when `df == 0`. TYPE: `float`

ece ¶

ece(y: object, p: object, *, n_bins: int = 15, strategy: str = 'mass', norm: str = 'l1', sample_weight: object = None) -> float

Expected calibration error; norm="max" gives the MCE.

Binning-sensitive and upward-biased in finite samples — report, never select on it (see the metrics chapter's table).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`n_bins`	Requested number of bins (positive integer). TYPE: `(int, keyword - only)` DEFAULT: `15`
`strategy`	`"mass"` (equal-count, default; the "adaptive ECE" of the literature) or `"width"` (equal-width over [0, 1]). TYPE: `(mass, width)` DEFAULT: `"mass"`
`norm`	`"l1"` (default, the usual ECE), `"l2"`, or `"max"` (the MCE). TYPE: `(l1, l2, max)` DEFAULT: `"l1"`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Weighted binned calibration error under the chosen norm.

Source code in src/probcal/metrics/binned.py

def ece(
    y: object,
    p: object,
    *,
    n_bins: int = 15,
    strategy: str = "mass",
    norm: str = "l1",
    sample_weight: object = None,
) -> float:
    """Expected calibration error; ``norm="max"`` gives the MCE.

    Binning-sensitive and upward-biased in finite samples — report, never
    select on it (see the metrics chapter's table).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    n_bins : int, keyword-only
        Requested number of bins (positive integer).
    strategy : {"mass", "width"}, keyword-only
        ``"mass"`` (equal-count, default; the "adaptive ECE" of the
        literature) or ``"width"`` (equal-width over [0, 1]).
    norm : {"l1", "l2", "max"}, keyword-only
        ``"l1"`` (default, the usual ECE), ``"l2"``, or ``"max"`` (the MCE).
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Weighted binned calibration error under the chosen norm.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    bs = bin_stats(y_arr, p_arr, w, n_bins, strategy)
    return _reduce(*_shares_gaps(bs, float(w.sum())), norm)

ece_debiased ¶

ece_debiased(y: object, p: object, *, n_bins: int = 15, strategy: str = 'mass', sample_weight: object = None) -> float

Bias-corrected ECE, floored at zero.

Per-bin squared gaps minus the within-bin variance of the event rate (correction in the spirit of Bröcker 2009 / Ferro & Fricker 2012): sqrt(max(gap_b**2 - rate_b * (1 - rate_b) / (n_b - 1), 0)) for bins with n_b > 1. n_b is the bin's effective count sum_b(w)**2 / sum_b(w**2) -- the size of an unweighted bin whose event rate has the same variance -- so it follows the weights inside the bin, and is the raw count for unit (or within-bin constant) weights.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`n_bins`	Requested number of bins. TYPE: `(int, keyword - only)` DEFAULT: `15`
`strategy`	`"mass"` (equal-count, default) or `"width"` (equal-width over [0, 1]). TYPE: `(mass, width)` DEFAULT: `"mass"`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Bias-corrected calibration error.

Source code in src/probcal/metrics/binned.py

def ece_debiased(
    y: object,
    p: object,
    *,
    n_bins: int = 15,
    strategy: str = "mass",
    sample_weight: object = None,
) -> float:
    """Bias-corrected ECE, floored at zero.

    Per-bin squared gaps minus the within-bin variance of the event rate
    (correction in the spirit of Bröcker 2009 / Ferro & Fricker 2012):
    ``sqrt(max(gap_b**2 - rate_b * (1 - rate_b) / (n_b - 1), 0))`` for bins
    with ``n_b > 1``. ``n_b`` is the bin's effective count
    ``sum_b(w)**2 / sum_b(w**2)`` -- the size of an unweighted bin whose
    event rate has the same variance -- so it follows the weights inside the
    bin, and is the raw count for unit (or within-bin constant) weights.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    n_bins : int, keyword-only
        Requested number of bins.
    strategy : {"mass", "width"}, keyword-only
        ``"mass"`` (equal-count, default) or ``"width"`` (equal-width over [0, 1]).
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Bias-corrected calibration error.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    return _debiased(bin_stats(y_arr, p_arr, w, n_bins, strategy), w)

ece_sweep ¶

ece_sweep(y: object, p: object, *, norm: str = 'l1', sample_weight: object = None, max_bins: int = 100, rule: str = 'first_violation') -> float

Monotonic-sweep calibration error (Roelofs et al., 2022).

Equal-mass bins; the bin count is the largest b such that b and every smaller bin count give non-decreasing bin event rates — the paper's b* = max{b : for all b' <= b, ybar_1 <= ... <= ybar_b'} (its Algorithm 1: sweep b = 2, 3, ... and stop at the first non-monotone binning). The paper sweeps up to b = n; the sweep here stops at max_bins (default 100) as a cost cap.

rule="largest" instead takes the largest monotone b anywhere in 2..max_bins, skipping over non-monotone counts.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`norm`	Norm of the final binned error at the selected bin count. TYPE: `(l1, l2, max)` DEFAULT: `"l1"`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`max_bins`	Upper end of the sweep (capped at `n`). TYPE: `(int, keyword - only)` DEFAULT: `100`
`rule`	`"first_violation"` (default, the paper's rule) or `"largest"` (the largest monotone count, non-monotone counts skipped). TYPE: `(first_violation, largest)` DEFAULT: `"first_violation"`

RETURNS	DESCRIPTION
`float`	Calibration error at the selected bin count.

Source code in src/probcal/metrics/binned.py

def ece_sweep(
    y: object,
    p: object,
    *,
    norm: str = "l1",
    sample_weight: object = None,
    max_bins: int = 100,
    rule: str = "first_violation",
) -> float:
    """Monotonic-sweep calibration error (Roelofs et al., 2022).

    Equal-mass bins; the bin count is the largest ``b`` such that ``b`` *and
    every smaller bin count* give non-decreasing bin event rates — the
    paper's ``b* = max{b : for all b' <= b, ybar_1 <= ... <= ybar_b'}``
    (its Algorithm 1: sweep ``b = 2, 3, ...`` and stop at the first
    non-monotone binning). The paper sweeps up to ``b = n``; the sweep here
    stops at ``max_bins`` (default 100) as a cost cap.

    ``rule="largest"`` instead takes the largest monotone ``b`` anywhere in
    ``2..max_bins``, skipping over non-monotone counts.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    norm : {"l1", "l2", "max"}, keyword-only
        Norm of the final binned error at the selected bin count.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.
    max_bins : int, keyword-only
        Upper end of the sweep (capped at ``n``).
    rule : {"first_violation", "largest"}, keyword-only
        ``"first_violation"`` (default, the paper's rule) or ``"largest"``
        (the largest monotone count, non-monotone counts skipped).

    Returns
    -------
    float
        Calibration error at the selected bin count.
    """
    if rule not in _SWEEP_RULES:
        raise ValueError(f"rule must be one of {_SWEEP_RULES}, got {rule!r}")
    cap = validate_positive_int(max_bins, "max_bins")
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    bs = bin_stats(y_arr, p_arr, w, _sweep_best_b(y_arr, p_arr, w, cap, rule), "mass")
    return _reduce(*_shares_gaps(bs, float(w.sum())), norm)

hosmer_lemeshow ¶

hosmer_lemeshow(y: object, p: object, *, n_bins: int = 10, sample_weight: object = None, ddof: int = 0) -> HosmerLemeshowResult

Hosmer–Lemeshow goodness-of-fit test on n_bins equal-mass risk groups.

The statistic depends on an essentially arbitrary grouping and its power scales with n — see the metrics chapter for why this is report-only.

Per group g, O_g = sum(w y) and E_g = sum(w p) are compared against V_g = S_g * pt_g * (1 - pt_g), where S_g = sum(w**2) and pt_g = sum(w**2 p) / S_g. V_g is the weighted form of the classical n_g pbar_g (1 - pbar_g): it equals the exact null variance sum(w**2 p (1 - p)) of O_g - E_g plus the same within-group spread term sum(w**2 (p - pt_g)**2) the unweighted statistic carries. The statistic is invariant to rescaling the weights, and unit weights give the unweighted test bit for bit.

The reference distribution is chi-square on groups - ddof degrees of freedom. With p fixed in advance (a model or calibrator fitted on other data, the usual validation setting) the per-group standardized residuals are asymptotically independent, so the statistic is chi2(groups): ddof=0, the default (Hosmer, Lemeshow and Sturdivant, 2013, on external validation). The classical groups - 2 is for a logistic model fitted on the same data; pass ddof=2 for it. On out-of-sample p it is anti-conservative (about 12% rejections at a nominal 5% with 10 groups).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`n_bins`	Requested number of equal-mass risk groups. Tied scores can collapse groups; when no more than `ddof` groups are non-empty the test has no degrees of freedom and returns `df=0`, `p_value=nan`. TYPE: `(int, keyword - only)` DEFAULT: `10`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`ddof`	Degrees of freedom subtracted from the number of non-empty groups: `0` (default) for `p` fixed before seeing `y`, `2` for a logistic model fitted on the same `(y, p)` data. TYPE: `(int, keyword - only)` DEFAULT: `0`

RETURNS	DESCRIPTION
`HosmerLemeshowResult`	Chi-square statistic, degrees of freedom, and p-value.

Source code in src/probcal/metrics/binned.py

def hosmer_lemeshow(
    y: object,
    p: object,
    *,
    n_bins: int = 10,
    sample_weight: object = None,
    ddof: int = 0,
) -> HosmerLemeshowResult:
    """Hosmer–Lemeshow goodness-of-fit test on ``n_bins`` equal-mass risk groups.

    The statistic depends on an essentially arbitrary grouping and its power
    scales with n — see the metrics chapter for why this is report-only.

    Per group ``g``, ``O_g = sum(w y)`` and ``E_g = sum(w p)`` are compared
    against ``V_g = S_g * pt_g * (1 - pt_g)``, where ``S_g = sum(w**2)`` and
    ``pt_g = sum(w**2 p) / S_g``. ``V_g`` is the weighted form of the
    classical ``n_g pbar_g (1 - pbar_g)``: it equals the exact null variance
    ``sum(w**2 p (1 - p))`` of ``O_g - E_g`` plus the same within-group
    spread term ``sum(w**2 (p - pt_g)**2)`` the unweighted statistic carries.
    The statistic is invariant to rescaling the weights, and unit weights
    give the unweighted test bit for bit.

    The reference distribution is chi-square on ``groups - ddof`` degrees of
    freedom. With ``p`` fixed in advance (a model or calibrator fitted on
    other data, the usual validation setting) the per-group standardized
    residuals are asymptotically independent, so the statistic is
    ``chi2(groups)``: ``ddof=0``, the default (Hosmer, Lemeshow and
    Sturdivant, 2013, on external validation). The classical ``groups - 2``
    is for a logistic model fitted on the same data; pass ``ddof=2`` for it.
    On out-of-sample ``p`` it is anti-conservative (about 12% rejections at
    a nominal 5% with 10 groups).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    n_bins : int, keyword-only
        Requested number of equal-mass risk groups. Tied scores can collapse
        groups; when no more than ``ddof`` groups are non-empty the test has
        no degrees of freedom and returns ``df=0``, ``p_value=nan``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.
    ddof : int, keyword-only
        Degrees of freedom subtracted from the number of non-empty groups:
        ``0`` (default) for ``p`` fixed before seeing ``y``, ``2`` for a
        logistic model fitted on the same ``(y, p)`` data.

    Returns
    -------
    HosmerLemeshowResult
        Chi-square statistic, degrees of freedom, and p-value.
    """
    ddof = validate_nonneg_int(ddof, "ddof")
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    idx, m = bin_index(p_arr, n_bins, "mass")
    keep = np.bincount(idx, minlength=m) > 0
    w2 = w * w
    obs = np.bincount(idx, weights=w * y_arr, minlength=m)[keep]
    exp = np.bincount(idx, weights=w * p_arr, minlength=m)[keep]
    s2 = np.bincount(idx, weights=w2, minlength=m)[keep]
    s2p = np.bincount(idx, weights=w2 * p_arr, minlength=m)[keep]
    denom = s2p * (1.0 - s2p / s2)
    ok = denom > 0
    stat = float(np.sum((obs[ok] - exp[ok]) ** 2 / denom[ok]))
    df = max(len(s2) - ddof, 0)
    p_value = chi2_sf(stat, df) if df > 0 else math.nan
    return HosmerLemeshowResult(statistic=stat, df=df, p_value=p_value)

smooth ¶

Binning-free calibration metrics: smoothECE, ECCE, ICI family, Spiegelhalter z.

Theory: docs/concepts/metrics.md. Sample weights follow the package convention in :mod:probcal.metrics._common.

Lattice `module-attribute` ¶

Lattice = tuple[float, float, int]

(width, t_lo, bins) of the equal-width logit lattice a smECE was solved on.

EcceResult `dataclass` ¶

EcceResult(stat_max: float, stat_mean: float)

Empirical cumulative calibration error: Kolmogorov-style max and mean of the cumulative deviation over sorted predictions.

ATTRIBUTE	DESCRIPTION
`stat_max`	Maximum absolute cumulative deviation. TYPE: `float`
`stat_mean`	Weighted mean absolute cumulative deviation (each tied-score block weighted by its total sample weight). TYPE: `float`

SpiegelhalterResult `dataclass` ¶

SpiegelhalterResult(z: float, p_value: float)

Spiegelhalter's z test of forecast unbiasedness (two-sided).

ATTRIBUTE	DESCRIPTION
`z`	Standardized test statistic (`nan` when its null variance is 0). TYPE: `float`
`p_value`	Two-sided p-value under the standard normal approximation (`nan` when `z` is). TYPE: `float`

smooth_ece ¶

smooth_ece(y: object, p: object, *, sample_weight: object = None, bins: int | None = 8192) -> float

Kernel-smoothed ECE with a self-consistent bandwidth (Błasiok–Nakkiran).

Residuals are smoothed with a Gaussian kernel on the logit scale (the paper's reflected kernel is a boundary device for [0, 1]; on the unbounded logit scale no reflection is needed), and the reported value is the fixed point smECE(sigma) = sigma found by bisection.

bins aggregates the weighted residual measure onto an equal-width logit lattice and solves the fixed point there by direct convolution (cost independent of n and sigma). If the solved sigma spans fewer than 8 cells, the solve is repeated once on a lattice refined to 8 cells per sigma; if that would exceed 2**20 cells or is still under-resolved -- or with bins=None or a degenerate logit range -- the exact computation (a 257-point grid evaluator, O(n) per bisection step) runs instead. The lattice value differs from the exact one by ~1e-4 on typical portfolios; on wide clipped-logit data the exact path's fixed grid under-resolves small-sigma kernels and the lattice value is the better one.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`bins`	Number of lattice bins for the fast path (default 8192); `None` forces the exact O(n) computation. TYPE: `(int or None, keyword - only)` DEFAULT: `8192`

RETURNS	DESCRIPTION
`float`	smECE: the fixed point `sigma` solving `smECE(sigma) = sigma`.

Source code in src/probcal/metrics/smooth.py

def smooth_ece(
    y: object, p: object, *, sample_weight: object = None, bins: int | None = 8192
) -> float:
    """Kernel-smoothed ECE with a self-consistent bandwidth (Błasiok–Nakkiran).

    Residuals are smoothed with a Gaussian kernel on the logit scale (the
    paper's reflected kernel is a boundary device for [0, 1]; on the
    unbounded logit scale no reflection is needed), and the
    reported value is the fixed point ``smECE(sigma) = sigma`` found by
    bisection.

    ``bins`` aggregates the weighted residual measure onto an equal-width
    logit lattice and solves the fixed point there by direct convolution
    (cost independent of n and sigma). If the solved ``sigma`` spans fewer
    than 8 cells, the solve is repeated once on a lattice refined to 8 cells
    per sigma; if that would exceed ``2**20`` cells or is still
    under-resolved -- or with ``bins=None`` or a degenerate logit range --
    the exact computation (a 257-point grid evaluator, O(n) per bisection
    step) runs instead. The lattice value differs from the exact one by
    ~1e-4 on typical portfolios; on wide clipped-logit data the exact path's
    fixed grid under-resolves small-sigma kernels and the lattice value is
    the better one.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.
    bins : int or None, keyword-only
        Number of lattice bins for the fast path (default 8192); ``None``
        forces the exact O(n) computation.

    Returns
    -------
    float
        smECE: the fixed point ``sigma`` solving ``smECE(sigma) = sigma``.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    t = logit(p_arr)
    mass = (w / w.sum()) * (y_arr - p_arr)
    return _smece_solve(t, mass, bins)[0]

ecce ¶

ecce(y: object, p: object, *, sample_weight: object = None) -> EcceResult

Cumulative-deviation calibration error (Arrieta-Ibarra et al., 2022).

Sort by prediction and walk the cumulative sum of weighted residuals, normalized by the total weight; under calibration the walk hovers near zero, and drift localizes miscalibration without any smoothing parameter. The walk is read only at the end of each block of tied predictions, so the result does not depend on the input order of tied rows; stat_mean weights each block by its total sample weight. Without ties and with unit weights the walk is read after every row.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`EcceResult`	Max and mean absolute cumulative deviation.

Source code in src/probcal/metrics/smooth.py

def ecce(y: object, p: object, *, sample_weight: object = None) -> EcceResult:
    """Cumulative-deviation calibration error (Arrieta-Ibarra et al., 2022).

    Sort by prediction and walk the cumulative sum of weighted residuals,
    normalized by the total weight; under calibration the walk hovers near
    zero, and drift localizes miscalibration without any smoothing
    parameter. The walk is read only at the end of each block of tied
    predictions, so the result does not depend on the input order of tied
    rows; ``stat_mean`` weights each block by its total sample weight.
    Without ties and with unit weights the walk is read after every row.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    EcceResult
        Max and mean absolute cumulative deviation.
    """
    return _ecce_stats(_ecce_walk(*_prep(y, p, sample_weight)))

ici ¶

ici(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None, grid_size: int | None = 512) -> float

Integrated calibration index: weighted mean |LOESS(y|p) - p| (Austin & Steyerberg, 2019).

The LOESS stage itself is unweighted. grid_size=None fits the LOESS exactly at every prediction.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`frac`	LOESS smoothing fraction. TYPE: `(float, keyword - only)` DEFAULT: `0.75`
`sample_weight`	Optional positive weights, same length as `y`; weights only the final averaging step, not the LOESS fit. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`grid_size`	LOESS evaluation grid size; `None` fits exactly at every prediction. TYPE: `(int or None, keyword - only)` DEFAULT: `512`

RETURNS	DESCRIPTION
`float`	Weighted mean absolute LOESS-to-prediction distance.

Source code in src/probcal/metrics/smooth.py

def ici(
    y: object,
    p: object,
    *,
    frac: float = 0.75,
    sample_weight: object = None,
    grid_size: int | None = 512,
) -> float:
    """Integrated calibration index: weighted mean |LOESS(y|p) - p|
    (Austin & Steyerberg, 2019).

    The LOESS stage itself is unweighted.
    ``grid_size=None`` fits the LOESS exactly at every prediction.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    frac : float, keyword-only
        LOESS smoothing fraction.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``; weights only the
        final averaging step, not the LOESS fit.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` fits exactly at every prediction.

    Returns
    -------
    float
        Weighted mean absolute LOESS-to-prediction distance.
    """
    return _ici_public("ici", y, p, frac, sample_weight, grid_size)

e50 ¶

e50(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None, grid_size: int | None = 512) -> float

Median of the |LOESS(y|p) - p| distances.

grid_size=None fits the LOESS exactly at every prediction. The LOESS distances are always unweighted; sample_weight, when given and not uniform, weights only the quantile step (see :func:~probcal._math.weighted_quantile).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`frac`	LOESS smoothing fraction. TYPE: `(float, keyword - only)` DEFAULT: `0.75`
`sample_weight`	Optional positive weights; used only for the quantile step. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`grid_size`	LOESS evaluation grid size; `None` fits exactly at every prediction. TYPE: `(int or None, keyword - only)` DEFAULT: `512`

RETURNS	DESCRIPTION
`float`	Median of the LOESS distances.

Source code in src/probcal/metrics/smooth.py

def e50(
    y: object,
    p: object,
    *,
    frac: float = 0.75,
    sample_weight: object = None,
    grid_size: int | None = 512,
) -> float:
    """Median of the |LOESS(y|p) - p| distances.

    ``grid_size=None`` fits the LOESS exactly at every prediction. The LOESS distances are
    always unweighted; ``sample_weight``, when given and
    not uniform, weights only the quantile step (see
    :func:`~probcal._math.weighted_quantile`).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    frac : float, keyword-only
        LOESS smoothing fraction.
    sample_weight : array_like or None, keyword-only
        Optional positive weights; used only for the quantile step.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` fits exactly at every prediction.

    Returns
    -------
    float
        Median of the LOESS distances.
    """
    return _ici_public("e50", y, p, frac, sample_weight, grid_size)

e90 ¶

e90(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None, grid_size: int | None = 512) -> float

90th percentile of the |LOESS(y|p) - p| distances.

grid_size=None fits the LOESS exactly at every prediction. The LOESS distances are always unweighted; sample_weight, when given and not uniform, weights only the quantile step (see :func:~probcal._math.weighted_quantile).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`frac`	LOESS smoothing fraction. TYPE: `(float, keyword - only)` DEFAULT: `0.75`
`sample_weight`	Optional positive weights; used only for the quantile step. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`grid_size`	LOESS evaluation grid size; `None` fits exactly at every prediction. TYPE: `(int or None, keyword - only)` DEFAULT: `512`

RETURNS	DESCRIPTION
`float`	90th percentile of the LOESS distances.

Source code in src/probcal/metrics/smooth.py

def e90(
    y: object,
    p: object,
    *,
    frac: float = 0.75,
    sample_weight: object = None,
    grid_size: int | None = 512,
) -> float:
    """90th percentile of the |LOESS(y|p) - p| distances.

    ``grid_size=None`` fits the LOESS exactly at every prediction. The LOESS distances are
    always unweighted; ``sample_weight``, when given and
    not uniform, weights only the quantile step (see
    :func:`~probcal._math.weighted_quantile`).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    frac : float, keyword-only
        LOESS smoothing fraction.
    sample_weight : array_like or None, keyword-only
        Optional positive weights; used only for the quantile step.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` fits exactly at every prediction.

    Returns
    -------
    float
        90th percentile of the LOESS distances.
    """
    return _ici_public("e90", y, p, frac, sample_weight, grid_size)

emax ¶

emax(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None, grid_size: int | None = 512) -> float

Maximum of the |LOESS(y|p) - p| distances.

grid_size=None fits the LOESS exactly at every prediction.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`frac`	LOESS smoothing fraction. TYPE: `(float, keyword - only)` DEFAULT: `0.75`
`sample_weight`	Optional positive weights, validated like everywhere else. The weighted maximum over strictly positive weights is the plain maximum, so the value does not depend on them. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`grid_size`	LOESS evaluation grid size; `None` fits exactly at every prediction. TYPE: `(int or None, keyword - only)` DEFAULT: `512`

RETURNS	DESCRIPTION
`float`	Maximum of the LOESS distances.

Source code in src/probcal/metrics/smooth.py

def emax(
    y: object,
    p: object,
    *,
    frac: float = 0.75,
    sample_weight: object = None,
    grid_size: int | None = 512,
) -> float:
    """Maximum of the |LOESS(y|p) - p| distances.

    ``grid_size=None`` fits the LOESS exactly at every prediction.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    frac : float, keyword-only
        LOESS smoothing fraction.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, validated like everywhere else. The
        weighted maximum over strictly positive weights is the plain maximum,
        so the value does not depend on them.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` fits exactly at every prediction.

    Returns
    -------
    float
        Maximum of the LOESS distances.
    """
    return _ici_public("emax", y, p, frac, sample_weight, grid_size)

spiegelhalter_z ¶

spiegelhalter_z(y: object, p: object, *, sample_weight: object = None) -> SpiegelhalterResult

Spiegelhalter (1986) z statistic built on the Brier score.

The numerator has expectation zero under calibration; the statistic is asymptotically standard normal. No binning, no smoothing; aggregates the whole range, so compensating regional errors can cancel.

With weights the numerator is sum(w (y - p) (1 - 2p)) and the null variance its exact variance under calibration, sum(w**2 (1 - 2p)^2 p (1 - p)) (package convention, :mod:probcal.metrics._common): z is invariant to rescaling the weights and unit weights give the unweighted statistic. A Kish-rescaled effective sample size is not used because it misstates the variance when the weights vary with p (1 - p).

The null variance is zero when every prediction is exactly 0.5 (or clipped to 0/1); the statistic is then undefined and z and p_value are nan.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`SpiegelhalterResult`	Z statistic and two-sided p-value.

Source code in src/probcal/metrics/smooth.py

def spiegelhalter_z(y: object, p: object, *, sample_weight: object = None) -> SpiegelhalterResult:
    """Spiegelhalter (1986) z statistic built on the Brier score.

    The numerator has expectation zero under calibration; the statistic is
    asymptotically standard normal. No binning, no smoothing; aggregates the
    whole range, so compensating regional errors can cancel.

    With weights the numerator is ``sum(w (y - p) (1 - 2p))`` and the null
    variance its exact variance under calibration, ``sum(w**2 (1 - 2p)^2 p
    (1 - p))`` (package convention, :mod:`probcal.metrics._common`): ``z``
    is invariant to rescaling the weights and unit weights give the
    unweighted statistic. A Kish-rescaled effective sample size is not used
    because it misstates the variance when the weights vary with ``p (1 - p)``.

    The null variance is zero when every prediction is exactly 0.5 (or
    clipped to 0/1); the statistic is then undefined and ``z`` and
    ``p_value`` are ``nan``.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    SpiegelhalterResult
        Z statistic and two-sided p-value.
    """
    return _spiegelhalter(*_prep(y, p, sample_weight))

regression ¶

Recalibration-regression framework: calibration intercept, slope, and joint test.

The Cox (1958) framework; lineage through Miller, Hui & Tierney (1991). Theory: docs/concepts/metrics.md. Sample weights follow the package convention in :mod:probcal.metrics._common.

CalibrationTestResult `dataclass` ¶

CalibrationTestResult(statistic: float, p_value: float, alpha: float, beta: float, df: float = 2.0)

Likelihood-ratio test of (intercept, slope) = (0, 1) — the Cox-framed weak calibration test.

ATTRIBUTE	DESCRIPTION
`statistic`	Likelihood-ratio test statistic, referred to `chi2(df)` (the Rao-Scott corrected statistic for non-uniform weights). TYPE: `float`
`p_value`	Upper-tail p-value of the statistic. TYPE: `float`
`alpha`	Fitted intercept (also available as :attr:`intercept`). TYPE: `float`
`beta`	Fitted slope (also available as :attr:`slope`); `nan` when `logit(p)` is constant. TYPE: `float`
`df`	Degrees of freedom of the reference chi-square: `2` (`1` when only the intercept is testable), smaller and fractional after the Rao-Scott correction for non-uniform weights. TYPE: `float`

intercept `property` ¶

intercept: float

Fitted intercept; the name used by the rest of the package for alpha.

slope `property` ¶

slope: float

Fitted slope; the name used by the rest of the package for beta.

GuardrailReport `dataclass` ¶

GuardrailReport(slope: float, intercept: float, spiegelhalter_p: float, slope_ok: bool, intercept_ok: bool, spiegelhalter_ok: bool, all_ok: bool)

Three-flag calibration health summary used across the package.

Thresholds are conventions, not theorems: slope within [0.9, 1.1], intercept within +/-0.1 log-odds, Spiegelhalter p above 0.05.

ATTRIBUTE	DESCRIPTION
`slope`	Fitted Cox calibration slope. TYPE: `float`
`intercept`	Fitted calibration-in-the-large intercept (log-odds). TYPE: `float`
`spiegelhalter_p`	Spiegelhalter test p-value. TYPE: `float`
`slope_ok`	Whether `slope` lies in `[0.9, 1.1]`. TYPE: `bool`
`intercept_ok`	Whether `abs(intercept) <= 0.1`. TYPE: `bool`
`spiegelhalter_ok`	Whether `spiegelhalter_p > 0.05`. TYPE: `bool`
`all_ok`	Conjunction of the three flags above. TYPE: `bool`

calibration_intercept ¶

calibration_intercept(y: object, p: object, *, sample_weight: object = None) -> float

Calibration-in-the-large in log-odds.

Logistic intercept with the slope fixed at 1 (offset regression on logit(p)).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Fitted intercept in log-odds units.

Source code in src/probcal/metrics/regression.py

def calibration_intercept(y: object, p: object, *, sample_weight: object = None) -> float:
    """Calibration-in-the-large in log-odds.

    Logistic intercept with the slope fixed at 1 (offset regression on
    logit(p)).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Fitted intercept in log-odds units.
    """
    return _intercept(*_prep(y, p, sample_weight))

calibration_slope ¶

calibration_slope(y: object, p: object, *, sample_weight: object = None) -> float

Cox calibration slope.

< 1 means overfitting/overconfident spread, > 1 underfitting. When every prediction has the same logit(p) (a constant p, or a single PD level) the slope is not identified and nan is returned with a UserWarning.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`float`	Fitted slope on the logit scale (`nan` when `logit(p)` is constant).

Source code in src/probcal/metrics/regression.py

def calibration_slope(y: object, p: object, *, sample_weight: object = None) -> float:
    """Cox calibration slope.

    ``< 1`` means overfitting/overconfident spread, ``> 1`` underfitting.
    When every prediction has the same ``logit(p)`` (a constant ``p``, or a
    single PD level) the slope is not identified and ``nan`` is returned
    with a ``UserWarning``.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    float
        Fitted slope on the logit scale (``nan`` when ``logit(p)`` is
        constant).
    """
    return float(_slope_fit(*_prep(y, p, sample_weight))[1][1])

calibration_test ¶

calibration_test(y: object, p: object, *, sample_weight: object = None) -> CalibrationTestResult

Likelihood-ratio test of joint calibration (alpha, beta) = (0, 1).

With weights, the fit and the likelihood ratio use sample_weight as given, and the weighted LR statistic -- distributed under calibration as a weighted sum of chi2_1 variables, not as chi2_2 -- is referred to a second-order Rao-Scott corrected chi-square (:func:probcal.metrics._common.rao_scott; the design effects use the exact null variance sum(w**2 p (1 - p) x x') of the weighted score). The test is invariant to rescaling the weights, and unit (or constant) weights give the plain 2-df likelihood-ratio test.

When every prediction has the same logit(p) the slope is not identified: only calibration-in-the-large is tested (slope fixed at 1, df = 1), beta is nan, and a UserWarning is emitted.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`CalibrationTestResult`	Test statistic, p-value, and fitted intercept/slope.

Source code in src/probcal/metrics/regression.py

def calibration_test(
    y: object, p: object, *, sample_weight: object = None
) -> CalibrationTestResult:
    """Likelihood-ratio test of joint calibration (alpha, beta) = (0, 1).

    With weights, the fit and the likelihood ratio use ``sample_weight`` as
    given, and the weighted LR statistic -- distributed under calibration as
    a weighted sum of ``chi2_1`` variables, not as ``chi2_2`` -- is referred
    to a second-order Rao-Scott corrected chi-square
    (:func:`probcal.metrics._common.rao_scott`; the design effects use the
    exact null variance ``sum(w**2 p (1 - p) x x')`` of the weighted score).
    The test is invariant to rescaling the weights, and unit (or constant)
    weights give the plain 2-df likelihood-ratio test.

    When every prediction has the same ``logit(p)`` the slope is not
    identified: only calibration-in-the-large is tested (slope fixed at 1,
    ``df = 1``), ``beta`` is ``nan``, and a ``UserWarning`` is emitted.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    CalibrationTestResult
        Test statistic, p-value, and fitted intercept/slope.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    X, beta_hat = _slope_fit(y_arr, p_arr, w, "calibration_test")
    alpha, beta = float(beta_hat[0]), float(beta_hat[1])

    if math.isnan(beta):  # constant logit(p): intercept-only (offset) test
        X = X[:, :1]
        fitted = expit(alpha + logit(p_arr))
    else:
        fitted = expit(X @ beta_hat)
    lr = max(2.0 * (bern_loglik(y_arr, fitted, w) - bern_loglik(y_arr, p_arr, w)), 0.0)
    stat, df = rao_scott(lr, X, p_arr, w)
    p_value = chi2_sf(stat, df)
    return CalibrationTestResult(statistic=stat, p_value=p_value, alpha=alpha, beta=beta, df=df)

calibration_guardrails ¶

calibration_guardrails(y: object, p: object, *, sample_weight: object = None) -> GuardrailReport

Evaluate the three guardrail flags.

Printed in selection reports and offset audit reports.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`sample_weight`	Optional positive weights, same length as `y`. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`

RETURNS	DESCRIPTION
`GuardrailReport`	Slope, intercept, and Spiegelhalter-p values with pass/fail flags.

Source code in src/probcal/metrics/regression.py

def calibration_guardrails(
    y: object, p: object, *, sample_weight: object = None
) -> GuardrailReport:
    """Evaluate the three guardrail flags.

    Printed in selection reports and offset audit reports.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    sample_weight : array_like or None, keyword-only
        Optional positive weights, same length as ``y``.

    Returns
    -------
    GuardrailReport
        Slope, intercept, and Spiegelhalter-p values with pass/fail flags.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    slope = float(_slope_fit(y_arr, p_arr, w, "calibration_guardrails")[1][1])
    intercept = _intercept(y_arr, p_arr, w)
    sp = _spiegelhalter(y_arr, p_arr, w)
    slope_ok = 0.9 <= slope <= 1.1
    intercept_ok = abs(intercept) <= 0.1
    sp_ok = sp.p_value > 0.05
    return GuardrailReport(
        slope=slope,
        intercept=intercept,
        spiegelhalter_p=sp.p_value,
        slope_ok=slope_ok,
        intercept_ok=intercept_ok,
        spiegelhalter_ok=sp_ok,
        all_ok=slope_ok and intercept_ok and sp_ok,
    )

grade ¶

Per-grade binomial and Jeffreys backtests (credit-risk rating grades).

Supervisory context (BCBS WP14; ECB 2019 instructions): each rating grade's realized default count is tested against its assigned PD. Theory: docs/concepts/metrics.md.

BinomialGradeResult `dataclass` ¶

BinomialGradeResult(grades: tuple, n: ndarray, k: ndarray, pd: ndarray, p_exact: ndarray, p_normal: ndarray, light: tuple, ci_low: ndarray, ci_high: ndarray)

Bases: _ResultBase

Exact and approximate one-sided binomial backtest per rating grade.

ATTRIBUTE	DESCRIPTION
`grades`	Grade labels: the `order` passed, else best to worst for a `Masterscale`, else sorted. TYPE: `tuple of str`
`n`	Observation count per grade. TYPE: `ndarray`
`k`	Default count per grade. TYPE: `ndarray`
`pd`	Assigned PD per grade (mean of `p` within the grade). TYPE: `ndarray`
`p_exact`	Exact binomial tail p-value per grade (also available as :attr:`p_value`, the name the other tests use). TYPE: `ndarray`
`p_normal`	Normal-approximation p-value per grade. TYPE: `ndarray`
`light`	Traffic light per grade (`"green"`, `"amber"`, or `"red"`), derived from `p_exact`. TYPE: `tuple of str`
`ci_low, ci_high`	90% Clopper-Pearson display interval for the observed rate. TYPE: `ndarray`

p_value `property` ¶

p_value: ndarray

Alias of :attr:p_exact (the exact test decides the traffic light).

JeffreysGradeResult `dataclass` ¶

JeffreysGradeResult(grades: tuple, n: ndarray, k: ndarray, pd: ndarray, p_value: ndarray, light: tuple, ci_low: ndarray, ci_high: ndarray)

Bases: _ResultBase

Jeffreys-posterior backtest per rating grade (ECB IRB practice).

ATTRIBUTE	DESCRIPTION
`grades`	Grade labels: the `order` passed, else best to worst for a `Masterscale`, else sorted. TYPE: `tuple of str`
`n`	Observation count per grade. TYPE: `ndarray`
`k`	Default count per grade. TYPE: `ndarray`
`pd`	Assigned PD per grade (mean of `p` within the grade). TYPE: `ndarray`
`p_value`	Posterior `P(theta <= PD \| k, n)` per grade. TYPE: `ndarray`
`light`	Traffic light per grade (`"green"`, `"amber"`, or `"red"`), derived from `p_value`. TYPE: `tuple of str`
`ci_low, ci_high`	Central 90% Jeffreys posterior display interval. TYPE: `ndarray`

binomial_grade_test ¶

binomial_grade_test(y: object, p: object, grades: object, *, sample_weight: object = None, order: object = None) -> BinomialGradeResult

Exact binomial tail test per grade: P(X >= k | n, PD).

Small p-values flag grades with more defaults than the assigned PD supports. The exact tail uses the incomplete-beta identity P(X >= k) = I_PD(k, n - k + 1); the normal approximation is reported alongside. Traffic lights: green > 0.05, amber > 0.01, red <= 0.01. ci_low/ci_high are 90% Clopper-Pearson display intervals for the observed rate; the traffic light itself remains the one-sided exact test, unchanged.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities (assigned PDs) in `[0, 1]`. TYPE: `array_like`
`grades`	Rating grade label per observation, or a :class:`probcal.Masterscale` that assigns them from `p` (results then come out best to worst). TYPE: `array_like or Masterscale`
`sample_weight`	Not used: grade tests use raw integer counts. A `UserWarning` is emitted if the weights are non-uniform. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`order`	Explicit grade order for the result (a permutation of the grade labels present); `None` keeps the masterscale order or sorts the labels. TYPE: `sequence of str or None, keyword-only` DEFAULT: `None`

RETURNS	DESCRIPTION
`BinomialGradeResult`	Per-grade counts, p-values, traffic lights, and display intervals.

Source code in src/probcal/metrics/grade.py

def binomial_grade_test(
    y: object,
    p: object,
    grades: object,
    *,
    sample_weight: object = None,
    order: object = None,
) -> BinomialGradeResult:
    """Exact binomial tail test per grade: P(X >= k | n, PD).

    Small p-values flag grades with more defaults than the assigned PD
    supports. The exact tail uses the incomplete-beta identity
    ``P(X >= k) = I_PD(k, n - k + 1)``; the normal approximation is reported
    alongside. Traffic lights: green > 0.05, amber > 0.01, red <= 0.01.
    ``ci_low``/``ci_high`` are 90% Clopper-Pearson display intervals for the
    observed rate; the traffic light itself remains the one-sided exact test,
    unchanged.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities (assigned PDs) in ``[0, 1]``.
    grades : array_like or Masterscale
        Rating grade label per observation, or a :class:`probcal.Masterscale`
        that assigns them from ``p`` (results then come out best to worst).
    sample_weight : array_like or None, keyword-only
        Not used: grade tests use raw integer counts. A ``UserWarning`` is
        emitted if the weights are non-uniform.
    order : sequence of str or None, keyword-only
        Explicit grade order for the result (a permutation of the grade
        labels present); ``None`` keeps the masterscale order or sorts the
        labels.

    Returns
    -------
    BinomialGradeResult
        Per-grade counts, p-values, traffic lights, and display intervals.
    """
    labels, n, k, pd = _grade_counts(y, p, grades, sample_weight, order)
    p_exact, p_normal, ci_low, ci_high = (np.empty(len(labels)) for _ in range(4))
    for i in range(len(labels)):
        ki, ni = int(k[i]), int(n[i])
        p_exact[i] = 1.0 if ki == 0 else float(betainc(float(ki), float(ni - ki + 1), pd[i]))
        se = np.sqrt(n[i] * pd[i] * (1.0 - pd[i]))
        p_normal[i] = norm_sf((k[i] - n[i] * pd[i]) / se if se > 0 else 0.0)
        ci_low[i] = 0.0 if ki == 0 else beta_ppf(0.05, float(ki), float(ni - ki + 1))
        ci_high[i] = 1.0 if ki == ni else beta_ppf(0.95, float(ki + 1), float(ni - ki))
    light = tuple(_traffic_light(v) for v in p_exact)
    return BinomialGradeResult(
        grades=labels,
        n=n,
        k=k,
        pd=pd,
        p_exact=p_exact,
        p_normal=p_normal,
        light=light,
        ci_low=ci_low,
        ci_high=ci_high,
    )

jeffreys_grade_test ¶

jeffreys_grade_test(y: object, p: object, grades: object, *, sample_weight: object = None, order: object = None) -> JeffreysGradeResult

Jeffreys test per grade: posterior P(theta <= PD | k, n) under Beta(k+1/2, n-k+1/2).

One-sided and conservative by design: a small value flags a grade whose PD is likely understated. Do not read it two-sided (a recurring validation error — see the metrics chapter). ci_low/ci_high are the central 90% Jeffreys posterior display intervals; the traffic light itself remains the one-sided posterior test, unchanged.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities (assigned PDs) in `[0, 1]`. TYPE: `array_like`
`grades`	Rating grade label per observation, or a :class:`probcal.Masterscale` that assigns them from `p` (results then come out best to worst). TYPE: `array_like or Masterscale`
`sample_weight`	Not used: grade tests use raw integer counts. A `UserWarning` is emitted if the weights are non-uniform. TYPE: `(array_like or None, keyword - only)` DEFAULT: `None`
`order`	Explicit grade order for the result (a permutation of the grade labels present); `None` keeps the masterscale order or sorts the labels. TYPE: `sequence of str or None, keyword-only` DEFAULT: `None`

RETURNS	DESCRIPTION
`JeffreysGradeResult`	Per-grade counts, p-values, traffic lights, and display intervals.

Source code in src/probcal/metrics/grade.py

def jeffreys_grade_test(
    y: object,
    p: object,
    grades: object,
    *,
    sample_weight: object = None,
    order: object = None,
) -> JeffreysGradeResult:
    """Jeffreys test per grade: posterior P(theta <= PD | k, n) under Beta(k+1/2, n-k+1/2).

    One-sided and conservative by design: a small value flags a grade whose
    PD is likely understated. Do not read it two-sided (a recurring
    validation error — see the metrics chapter). ``ci_low``/``ci_high`` are
    the central 90% Jeffreys posterior display intervals; the traffic light
    itself remains the one-sided posterior test, unchanged.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities (assigned PDs) in ``[0, 1]``.
    grades : array_like or Masterscale
        Rating grade label per observation, or a :class:`probcal.Masterscale`
        that assigns them from ``p`` (results then come out best to worst).
    sample_weight : array_like or None, keyword-only
        Not used: grade tests use raw integer counts. A ``UserWarning`` is
        emitted if the weights are non-uniform.
    order : sequence of str or None, keyword-only
        Explicit grade order for the result (a permutation of the grade
        labels present); ``None`` keeps the masterscale order or sorts the
        labels.

    Returns
    -------
    JeffreysGradeResult
        Per-grade counts, p-values, traffic lights, and display intervals.
    """
    labels, n, k, pd = _grade_counts(y, p, grades, sample_weight, order)
    p_value, ci_low, ci_high = (np.empty(len(labels)) for _ in range(3))
    for i in range(len(labels)):
        a, b = k[i] + 0.5, n[i] - k[i] + 0.5
        p_value[i] = float(betainc(a, b, pd[i]))
        ci_low[i] = beta_ppf(0.05, a, b)
        ci_high[i] = beta_ppf(0.95, a, b)
    light = tuple(_traffic_light(v) for v in p_value)
    return JeffreysGradeResult(
        grades=labels,
        n=n,
        k=k,
        pd=pd,
        p_value=p_value,
        light=light,
        ci_low=ci_low,
        ci_high=ci_high,
    )

kernel ¶

Kernel calibration error (SKCE) and its calibration tests.

Theory: docs/concepts/metrics.md. Widmann, Lindsten & Zachariah (2019), "Calibration tests in multi-class classification: A unifying framework", NeurIPS 32 (arXiv:1910.11385).

Binary specialization: with the identity-matrix kernel construction and a prediction represented as the 2-vector (1 - p, p), the paper's kernel term reduces to h_ij = 2 * k(s_i, s_j) * (y_i - p_i) * (y_j - p_j) — the factor 2 keeps values comparable with the paper's framework. Residuals always stay on the probability scale; only the kernel input s may be logit-transformed.

No sample_weight: the cited U-statistic theory (unbiasedness, the degenerate limit, the distribution-free bounds) is stated for unweighted i.i.d. samples. Refusing the argument is honest; improvising weighted inference is not.

Complexity: "uq", "biased", and the bootstrap test are O(n^2) memory and O(n_boot * n^2) time; prefer method="asymptotic" for n >~ 20 000.

SkceTestResult `dataclass` ¶

SkceTestResult(statistic: float, estimator: str, method: str, p_value: float, p_value_bound: float, bandwidth: float, n_boot: int | None)

One-sided SKCE calibration test (H0: calibrated; large positive rejects).

ATTRIBUTE	DESCRIPTION
`statistic`	SKCE point estimate (`"ul"` for the asymptotic method, `"uq"` for the bootstrap method). TYPE: `float`
`estimator`	Estimator used for `statistic` (`"ul"` or `"uq"`). TYPE: `str`
`method`	Test method used (`"asymptotic"` or `"bootstrap"`). TYPE: `str`
`p_value`	Test p-value. TYPE: `float`
`p_value_bound`	Distribution-free worst-case p-value bound (valid without asymptotics). TYPE: `float`
`bandwidth`	Kernel bandwidth used (resolved from `bandwidth=None` if applicable). TYPE: `float`
`n_boot`	Bootstrap replicate count for the bootstrap method; `None` for the asymptotic method. TYPE: `int or None`

skce ¶

skce(y: object, p: object, *, estimator: str = 'uq', kernel: str = 'laplace', bandwidth: float | None = None, scale: str = 'probability', random_state: int = 42) -> float

Squared kernel calibration error (Widmann et al., 2019, Table 1).

"uq" (default) is the unbiased quadratic estimator (may be negative); "ul" the unbiased linear O(n) estimator over seeded disjoint pairs (random_state controls the pairing); "biased" the nonnegative V-statistic. bandwidth=None uses the deterministic median heuristic; scale="logit" transforms the kernel input only (the low-PD option).

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`estimator`	Estimator variant; see above. TYPE: `(uq, ul, biased)` DEFAULT: `"uq"`
`kernel`	Kernel family applied to the (scaled) score distance. TYPE: `(laplace, gaussian)` DEFAULT: `"laplace"`
`bandwidth`	Kernel bandwidth; `None` (default) uses the median-heuristic distance (mean fallback if the median is 0). TYPE: `(float or None, keyword - only)` DEFAULT: `None`
`scale`	Scale on which the kernel input `s` is computed; residuals stay on the probability scale regardless. TYPE: `(probability, logit)` DEFAULT: `"probability"`
`random_state`	Seed for the `"ul"` estimator's disjoint-pairing permutation (unused by `"uq"`/`"biased"`). TYPE: `(int, keyword - only)` DEFAULT: `42`

RETURNS	DESCRIPTION
`float`	SKCE point estimate.

RAISES	DESCRIPTION
`ValueError`	If `estimator` is not one of `"uq"`, `"ul"`, `"biased"`, or if fewer than 2 observations are given.

Source code in src/probcal/metrics/kernel.py

def skce(
    y: object,
    p: object,
    *,
    estimator: str = "uq",
    kernel: str = "laplace",
    bandwidth: float | None = None,
    scale: str = "probability",
    random_state: int = 42,
) -> float:
    """Squared kernel calibration error (Widmann et al., 2019, Table 1).

    ``"uq"`` (default) is the unbiased quadratic estimator (may be negative);
    ``"ul"`` the unbiased linear O(n) estimator over seeded disjoint pairs
    (``random_state`` controls the pairing); ``"biased"`` the nonnegative
    V-statistic. ``bandwidth=None`` uses the deterministic median heuristic;
    ``scale="logit"`` transforms the kernel input only (the low-PD option).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    estimator : {"uq", "ul", "biased"}, keyword-only
        Estimator variant; see above.
    kernel : {"laplace", "gaussian"}, keyword-only
        Kernel family applied to the (scaled) score distance.
    bandwidth : float or None, keyword-only
        Kernel bandwidth; ``None`` (default) uses the median-heuristic
        distance (mean fallback if the median is 0).
    scale : {"probability", "logit"}, keyword-only
        Scale on which the kernel input ``s`` is computed; residuals stay on
        the probability scale regardless.
    random_state : int, keyword-only
        Seed for the ``"ul"`` estimator's disjoint-pairing permutation
        (unused by ``"uq"``/``"biased"``).

    Returns
    -------
    float
        SKCE point estimate.

    Raises
    ------
    ValueError
        If ``estimator`` is not one of ``"uq"``, ``"ul"``, ``"biased"``, or if
        fewer than 2 observations are given.
    """
    if estimator not in ("uq", "ul", "biased"):
        raise ValueError(f"estimator must be 'uq', 'ul', or 'biased', got {estimator!r}")
    y_arr, p_arr, s, bw = _setup(y, p, scale, bandwidth, 2, "skce")
    if estimator == "ul":
        return float(np.mean(_ul_terms(y_arr, p_arr, s, kernel, bw, random_state)))
    h = _h_full(y_arr, p_arr, s, kernel, bw)
    if estimator == "biased":
        return float(h.sum() / len(p_arr) ** 2)
    return _uq(h)

skce_test ¶

skce_test(y: object, p: object, *, method: str = 'bootstrap', n_boot: int = 999, kernel: str = 'laplace', bandwidth: float | None = None, scale: str = 'probability', random_state: int = 42) -> SkceTestResult

Calibration test on the SKCE (Widmann et al., 2019, Sec. 6 / App. G).

"bootstrap" (default): quadratic statistic with Arcones–Giné centered resampling; O(n_boot * n^2) — the more powerful choice. "asymptotic": linear statistic, normal approximation (Corollary G.3); O(n), preferred for n >~ 20 000, but a single random pairing can miss slope-type miscalibration that the bootstrap test rejects (the paper's documented power gap). p_value_bound is the distribution-free worst case.

PARAMETER	DESCRIPTION
`y`	Binary outcomes in `{0, 1}`. TYPE: `array_like`
`p`	Predicted probabilities in `[0, 1]`. TYPE: `array_like`
`method`	Test method; see above. TYPE: `(bootstrap, asymptotic)` DEFAULT: `"bootstrap"`
`n_boot`	Bootstrap replicate count (`"bootstrap"` method only). TYPE: `(int, keyword - only)` DEFAULT: `999`
`kernel`	Kernel family applied to the (scaled) score distance. TYPE: `(laplace, gaussian)` DEFAULT: `"laplace"`
`bandwidth`	Kernel bandwidth; `None` (default) uses the median-heuristic distance (mean fallback if the median is 0). TYPE: `(float or None, keyword - only)` DEFAULT: `None`
`scale`	Scale on which the kernel input `s` is computed; residuals stay on the probability scale regardless. TYPE: `(probability, logit)` DEFAULT: `"probability"`
`random_state`	Seed for the resampling (`"bootstrap"`) or disjoint-pairing (`"asymptotic"`) randomness. TYPE: `(int, keyword - only)` DEFAULT: `42`

RETURNS	DESCRIPTION
`SkceTestResult`	Test statistic, method, p-value, and worst-case bound.

Source code in src/probcal/metrics/kernel.py

def skce_test(
    y: object,
    p: object,
    *,
    method: str = "bootstrap",
    n_boot: int = 999,
    kernel: str = "laplace",
    bandwidth: float | None = None,
    scale: str = "probability",
    random_state: int = 42,
) -> SkceTestResult:
    """Calibration test on the SKCE (Widmann et al., 2019, Sec. 6 / App. G).

    ``"bootstrap"`` (default): quadratic statistic with Arcones–Giné centered
    resampling; O(n_boot * n^2) — the more powerful choice. ``"asymptotic"``:
    linear statistic, normal approximation (Corollary G.3); O(n), preferred
    for n >~ 20 000, but a single random pairing can miss slope-type
    miscalibration that the bootstrap test rejects (the paper's documented
    power gap). ``p_value_bound`` is the distribution-free worst case.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    method : {"bootstrap", "asymptotic"}, keyword-only
        Test method; see above.
    n_boot : int, keyword-only
        Bootstrap replicate count (``"bootstrap"`` method only).
    kernel : {"laplace", "gaussian"}, keyword-only
        Kernel family applied to the (scaled) score distance.
    bandwidth : float or None, keyword-only
        Kernel bandwidth; ``None`` (default) uses the median-heuristic
        distance (mean fallback if the median is 0).
    scale : {"probability", "logit"}, keyword-only
        Scale on which the kernel input ``s`` is computed; residuals stay on
        the probability scale regardless.
    random_state : int, keyword-only
        Seed for the resampling (``"bootstrap"``) or disjoint-pairing
        (``"asymptotic"``) randomness.

    Returns
    -------
    SkceTestResult
        Test statistic, method, p-value, and worst-case bound.
    """
    if method not in ("bootstrap", "asymptotic"):
        raise ValueError(f"method must be 'bootstrap' or 'asymptotic', got {method!r}")
    n_boot = validate_positive_int(n_boot, "n_boot")
    y_arr, p_arr, s, bw = _setup(y, p, scale, bandwidth, 4, "skce_test")
    n = len(p_arr)

    if method == "asymptotic":
        terms = _ul_terms(y_arr, p_arr, s, kernel, bw, random_state)
        stat = float(np.mean(terms))
        sd = float(np.std(terms, ddof=1))
        if sd == 0.0:
            p_value = 1.0 if stat <= 0.0 else 0.0
        else:
            z = math.sqrt(len(terms)) * stat / sd
            p_value = norm_sf(z)
        return SkceTestResult(
            statistic=stat,
            estimator="ul",
            method="asymptotic",
            p_value=p_value,
            p_value_bound=_p_value_bound(stat, n),
            bandwidth=bw,
            n_boot=None,
        )

    h = _h_full(y_arr, p_arr, s, kernel, bw)
    stat = _uq(h)
    t_obs = n * stat
    c = h.mean(axis=1)
    h_tilde = h - c[:, None] - c[None, :] + h.mean()
    rng = np.random.default_rng(random_state)
    counts = rng.multinomial(n, np.full(n, 1.0 / n), size=n_boot).astype(np.float64)
    quad = np.einsum("bi,ij,bj->b", counts, h_tilde, counts)
    # Same scale as t_obs = n * U_n = S / (n - 1), with S the off-diagonal sum:
    # the resampled off-diagonal sum over (n - 1). Dividing by n instead would
    # shrink every replicate by (n - 1) / n and make p-values anti-conservative.
    t_b = (quad - counts @ np.diag(h_tilde)) / (n - 1)
    p_value = float((1 + int(np.sum(t_b >= t_obs))) / (n_boot + 1))
    return SkceTestResult(
        statistic=stat,
        estimator="uq",
        method="bootstrap",
        p_value=p_value,
        p_value_bound=_p_value_bound(stat, n),
        bandwidth=bw,
        n_boot=n_boot,
    )

API: metrics and tests¶

metrics ¶

GroupedMetricReport dataclass ¶

to_frame ¶

ReliabilitySummary dataclass ¶

HosmerLemeshowResult dataclass ¶

BinomialGradeResult dataclass ¶

p_value property ¶

HlEResult dataclass ¶

interpret ¶

JeffreysGradeResult dataclass ¶

PlutoTascheResult dataclass ¶

interpret ¶

SkceTestResult dataclass ¶

CalibrationTestResult dataclass ¶

intercept property ¶

slope property ¶

GuardrailReport dataclass ¶

LogLossDecomposition dataclass ¶

MurphyCurve dataclass ¶

MurphyDecomposition dataclass ¶

EcceResult dataclass ¶

SpiegelhalterResult dataclass ¶

evaluate ¶

reliability_summary ¶

ece ¶

ece_debiased ¶

ece_sweep ¶

hosmer_lemeshow ¶

binomial_grade_test ¶

hl_e_test ¶

jeffreys_grade_test ¶

jeffreys_upper_bands ¶

pluto_tasche ¶

pluto_tasche_from_arrays ¶

skce ¶

skce_test ¶

calibration_guardrails ¶

calibration_intercept ¶

calibration_slope ¶

calibration_test ¶

brier_score ¶

brier_skill_score ¶

log_loss ¶

logloss_calibration_refinement ¶

murphy_curve ¶

murphy_decomposition ¶

e50 ¶

e90 ¶

ecce ¶

emax ¶

ici ¶

smooth_ece ¶

spiegelhalter_z ¶

scores ¶

MurphyDecomposition dataclass ¶

MurphyCurve dataclass ¶

LogLossDecomposition dataclass ¶

log_loss ¶

brier_score ¶

brier_skill_score ¶

murphy_decomposition ¶

murphy_curve ¶

logloss_calibration_refinement ¶

binned ¶

HosmerLemeshowResult dataclass ¶

ece ¶

ece_debiased ¶

ece_sweep ¶

hosmer_lemeshow ¶

smooth ¶

Lattice module-attribute ¶

EcceResult dataclass ¶

SpiegelhalterResult dataclass ¶

smooth_ece ¶

ecce ¶

ici ¶

e50 ¶

e90 ¶

emax ¶

GroupedMetricReport `dataclass` ¶

ReliabilitySummary `dataclass` ¶

HosmerLemeshowResult `dataclass` ¶

BinomialGradeResult `dataclass` ¶

p_value `property` ¶

HlEResult `dataclass` ¶

JeffreysGradeResult `dataclass` ¶

PlutoTascheResult `dataclass` ¶

SkceTestResult `dataclass` ¶

CalibrationTestResult `dataclass` ¶

intercept `property` ¶

slope `property` ¶

GuardrailReport `dataclass` ¶

LogLossDecomposition `dataclass` ¶

MurphyCurve `dataclass` ¶

MurphyDecomposition `dataclass` ¶

EcceResult `dataclass` ¶

SpiegelhalterResult `dataclass` ¶

MurphyDecomposition `dataclass` ¶

MurphyCurve `dataclass` ¶

LogLossDecomposition `dataclass` ¶

HosmerLemeshowResult `dataclass` ¶

Lattice `module-attribute` ¶

EcceResult `dataclass` ¶

SpiegelhalterResult `dataclass` ¶

CalibrationTestResult `dataclass` ¶

intercept `property` ¶

slope `property` ¶

GuardrailReport `dataclass` ¶

BinomialGradeResult `dataclass` ¶

p_value `property` ¶

JeffreysGradeResult `dataclass` ¶

SkceTestResult `dataclass` ¶