Skip to content

API: metrics and tests

metrics

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

evaluate lives here because it aggregates across every submodule. 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 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
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.
    """
    rows = [
        {"group": group, "metric": n, "value": v, "ci_low": lo, "ci_high": hi}
        for group, n, v, lo, hi in self._rows()
    ]
    try:
        import pandas as pd
    except ImportError:
        return rows
    return pd.DataFrame(rows)

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 (used groups minus 2, floored at 1).

TYPE: int

p_value

Upper-tail p-value of the chi-square statistic.

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)

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

ATTRIBUTE DESCRIPTION
grades

Grade labels: sorted when a label array was given, best to worst when a Masterscale was given.

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.

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

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, 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 
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
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)

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

ATTRIBUTE DESCRIPTION
grades

Grade labels: sorted when a label array was given, best to worst when a Masterscale was given.

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 only if pd_upper needed the running-maximum touch-up (see below) to stay non-decreasing best to worst. Pooled sets are nested (grade i's pooled set contains grade i + 1's), so for a portfolio whose observed per-grade default rates already respect rating order, pd_upper comes out non-decreasing on its own and this flag is False; a noisy grade whose own rate exceeds the worse-grade pool it joins can still produce a real (not merely floating-point) local dip. pd_upper is always non-decreasing on return either way: the raw bound is replaced by its cumulative maximum best to worst (a prudent hull), which never lowers any grade's bound -- only raises a grade whose raw bound fell below a better grade's, never the reverse.

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 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
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)

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

ATTRIBUTE DESCRIPTION
statistic

Likelihood-ratio test statistic (chi-square, 2 df).

TYPE: float

p_value

Upper-tail p-value of the statistic.

TYPE: float

alpha

Fitted intercept.

TYPE: float

beta

Fitted slope.

TYPE: float

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

Mean absolute cumulative deviation.

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.

TYPE: float

p_value

Two-sided p-value under the standard normal approximation.

TYPE: float

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

adaptive_ece 

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

Adaptive ECE: an explicit alias for equal-mass ece.

The literature uses both names for the same estimator.

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

norm

Norm passed through to :func:ece.

TYPE: (l1, l2, max) DEFAULT: "l1"

sample_weight

Optional non-negative weights, same length as y.

TYPE: (array_like or None, keyword - only) DEFAULT: None

RETURNS DESCRIPTION
float

Equal-mass binned calibration error under the chosen norm.
Source code in src/probcal/metrics/binned.py 
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
def adaptive_ece(
    y: object,
    p: object,
    *,
    n_bins: int = 15,
    norm: str = "l1",
    sample_weight: object = None,
) -> float:
    """Adaptive ECE: an explicit alias for equal-mass ``ece``.

    The literature uses both names for the same estimator.

    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.
    norm : {"l1", "l2", "max"}, keyword-only
        Norm passed through to :func:`ece`.
    sample_weight : array_like or None, keyword-only
        Optional non-negative weights, same length as ``y``.

    Returns
    -------
    float
        Equal-mass binned calibration error under the chosen norm.
    """
    return ece(y, p, n_bins=n_bins, strategy="mass", norm=norm, sample_weight=sample_weight)

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.

TYPE: (int, keyword - only) DEFAULT: 15

strategy

"mass" (equal-count, default) 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 non-negative 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 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
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.
    strategy : {"mass", "width"}, keyword-only
        ``"mass"`` (equal-count, default) 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 non-negative 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)
    shares, gaps, _, _ = _bin_gaps(y_arr, p_arr, w, n_bins, strategy)
    return _ece_from_gaps(shares, gaps, 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).

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 non-negative 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 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
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).

    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 non-negative weights, same length as ``y``.

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

ece_sweep 

ece_sweep(y: object, p: object, *, norm: str = 'l1', sample_weight: object = None) -> float

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

Uses equal-mass bins with the largest B whose bin event rates remain monotone non-decreasing (scan 2..min(n, 100)).

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

norm

Norm passed to the final :func:ece call at the selected bin count.

TYPE: (l1, l2, max) DEFAULT: "l1"

sample_weight

Optional non-negative weights, same length as y.

TYPE: (array_like or None, keyword - only) DEFAULT: None

RETURNS DESCRIPTION
float

Calibration error at the largest monotone bin count.
Source code in src/probcal/metrics/binned.py 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def ece_sweep(
    y: object,
    p: object,
    *,
    norm: str = "l1",
    sample_weight: object = None,
) -> float:
    """Monotonic-sweep calibration error (Roelofs et al., 2022).

    Uses equal-mass bins with the largest ``B`` whose bin event rates remain
    monotone non-decreasing (scan 2..min(n, 100)).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    norm : {"l1", "l2", "max"}, keyword-only
        Norm passed to the final :func:`ece` call at the selected bin count.
    sample_weight : array_like or None, keyword-only
        Optional non-negative weights, same length as ``y``.

    Returns
    -------
    float
        Calibration error at the largest monotone bin count.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    best_b = 1
    for b in range(2, min(len(p_arr), 100) + 1):
        _, _, rates, _ = _bin_gaps(y_arr, p_arr, w, b, "mass")
        if np.all(np.diff(rates) >= 0.0):
            best_b = b
    return _ece_at_best_b(y_arr, p_arr, w, best_b, norm)

hosmer_lemeshow 

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

Hosmer–Lemeshow goodness-of-fit test on g 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.

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

g

Requested number of equal-mass risk groups.

TYPE: (int, keyword - only) DEFAULT: 10

sample_weight

Optional non-negative weights, same length as y.

TYPE: (array_like or None, keyword - only) DEFAULT: None

RETURNS DESCRIPTION
HosmerLemeshowResult

Chi-square statistic, degrees of freedom, and p-value.
Source code in src/probcal/metrics/binned.py 
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
def hosmer_lemeshow(
    y: object,
    p: object,
    *,
    g: int = 10,
    sample_weight: object = None,
) -> HosmerLemeshowResult:
    """Hosmer–Lemeshow goodness-of-fit test on ``g`` 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.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    g : int, keyword-only
        Requested number of equal-mass risk groups.
    sample_weight : array_like or None, keyword-only
        Optional non-negative weights, same length as ``y``.

    Returns
    -------
    HosmerLemeshowResult
        Chi-square statistic, degrees of freedom, and p-value.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    idx, m = _bin_index(p_arr, g, "mass")
    stat = 0.0
    used = 0
    for b in range(m):
        mask = idx == b
        if not np.any(mask):
            continue
        nb = float(w[mask].sum())
        obs = float(np.sum(w[mask] * y_arr[mask]))
        exp = float(np.sum(w[mask] * p_arr[mask]))
        denom = exp * (1.0 - exp / nb)
        if denom > 0:
            stat += (obs - exp) ** 2 / denom
        used += 1
    df = max(used - 2, 1)
    p_value = 1.0 - float(gammainc_lower(df / 2.0, stat / 2.0))
    return HosmerLemeshowResult(statistic=float(stat), df=df, p_value=p_value)

binomial_grade_test 

binomial_grade_test(y: object, p: object, grades: object, *, sample_weight: 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

RETURNS DESCRIPTION
BinomialGradeResult

Per-grade counts, p-values, traffic lights, and display intervals.
Source code in src/probcal/metrics/grade.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def binomial_grade_test(
    y: object, p: object, grades: object, *, sample_weight: 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.

    Returns
    -------
    BinomialGradeResult
        Per-grade counts, p-values, traffic lights, and display intervals.
    """
    y_arr, p_arr, _ = _prep(y, p, None)
    _check_weights(sample_weight, len(y_arr))
    g_arr, order = _resolve_grades(grades, p_arr)
    labels, n, k, pd = _per_grade(y_arr, p_arr, g_arr, order)
    p_exact = np.empty(len(labels))
    p_normal = np.empty(len(labels))
    for i in range(len(labels)):
        if k[i] == 0:
            p_exact[i] = 1.0
        else:
            p_exact[i] = float(betainc(float(k[i]), float(n[i] - k[i] + 1), pd[i]))
        se = np.sqrt(n[i] * pd[i] * (1.0 - pd[i]))
        z = (k[i] - n[i] * pd[i]) / se if se > 0 else 0.0
        p_normal[i] = float(1.0 - norm_cdf(np.array([z]))[0])
    light = tuple(_traffic_light(v) for v in p_exact)
    ci_low = np.empty(len(labels))
    ci_high = np.empty(len(labels))
    for i in range(len(labels)):
        ki, ni = int(k[i]), int(n[i])
        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))
    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 (monitor._processes .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 on the Bernoulli factors (passed straight into bern_log_lr) -- consistent with how CalibrationMonitor and the rest of probcal.metrics treat weights, but note that non-integer weights break the interpretation of LR as a genuine likelihood ratio of independent Bernoulli draws (the same caveat docs/concepts/monitoring.md records for the monitor).

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 non-negative weights, same length as y.

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 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
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 (``monitor._processes
    .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 on the Bernoulli factors
    (passed straight into ``bern_log_lr``) -- consistent with how
    ``CalibrationMonitor`` and the rest of ``probcal.metrics`` treat
    weights, but note that non-integer weights break the interpretation of
    ``LR`` as a genuine likelihood ratio of independent Bernoulli draws
    (the same caveat ``docs/concepts/monitoring.md`` records for the
    monitor).

    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 non-negative weights, same length as ``y``.

    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)
    from .grade import _resolve_grades

    g_arr, order = _resolve_grades(grades, p_arr)
    if g_arr.ndim != 1 or len(g_arr) != len(y_arr):
        raise ValueError("grades must be a 1-D array matching y and p in length")
    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)

    g_str = g_arr
    if order is not None:
        labels = order
    else:
        labels = tuple(str(label) for label in sorted(np.unique(g_str)))
    log_e_grade = np.empty(len(labels))
    for i, label in enumerate(labels):
        mask = g_str == label
        y_g, p_g, z_g, w_g = y_arr[mask], p_arr[mask], z_arr[mask], 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) -> 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

RETURNS DESCRIPTION
JeffreysGradeResult

Per-grade counts, p-values, traffic lights, and display intervals.
Source code in src/probcal/metrics/grade.py 
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
def jeffreys_grade_test(
    y: object, p: object, grades: object, *, sample_weight: 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.

    Returns
    -------
    JeffreysGradeResult
        Per-grade counts, p-values, traffic lights, and display intervals.
    """
    y_arr, p_arr, _ = _prep(y, p, None)
    _check_weights(sample_weight, len(y_arr))
    g_arr, order = _resolve_grades(grades, p_arr)
    labels, n, k, pd = _per_grade(y_arr, p_arr, g_arr, order)
    p_value = np.empty(len(labels))
    for i in range(len(labels)):
        p_value[i] = float(betainc(k[i] + 0.5, n[i] - k[i] + 0.5, pd[i]))
    light = tuple(_traffic_light(v) for v in p_value)
    ci_low = np.empty(len(labels))
    ci_high = np.empty(len(labels))
    for i in range(len(labels)):
        a, b = k[i] + 0.5, n[i] - k[i] + 0.5
        ci_low[i] = beta_ppf(0.05, a, b)
        ci_high[i] = beta_ppf(0.95, a, b)
    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, *, level: 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(level, 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 non-decreasing (a noisy grade can have a smaller posterior upper bound than a better grade), which would make the bands overlap or invert. It is monotonized by :func:probcal._math.pava (weighted isotonic regression, weights = grade size n) in the given order — the minimum-adjustment non-decreasing fit, not a running maximum — with a UserWarning emitted only when that adjustment actually changed a value.

The resulting {grade: (lo, hi)} table is exactly the shape :func:probcal.thresholds.calibrated_bands_to_raw consumes 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

level

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) orders 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 to (lo, hi) calibrated-probability bounds, contiguous and non-decreasing best to worst.
RAISES DESCRIPTION
ValueError

If y/p/grades are not equal-length 1-D arrays, level 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, level=0.9)
>>> bands["A"][0]
0.0
>>> bands["A"][1] < bands["B"][1]
True
Source code in src/probcal/metrics/_conservative.py 
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
def jeffreys_upper_bands(
    y: object,
    p: object,
    grades: object,
    *,
    level: 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(level, 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 non-decreasing (a noisy grade
    can have a smaller posterior upper bound than a better grade), which
    would make the bands overlap or invert. It is monotonized by
    :func:`probcal._math.pava` (weighted isotonic regression, weights = grade
    size ``n``) in the given ``order`` — the minimum-adjustment non-decreasing
    fit, not a running maximum — with a ``UserWarning`` emitted only when
    that adjustment actually changed a value.

    The resulting ``{grade: (lo, hi)}`` table is exactly the shape
    :func:`probcal.thresholds.calibrated_bands_to_raw` consumes 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``.
    level : 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) orders grades by their mean
        ``p``, ascending (lowest predicted PD first).

    Returns
    -------
    dict[str, tuple[float, float]]
        Mapping of grade label to ``(lo, hi)`` calibrated-probability bounds,
        contiguous and non-decreasing best to worst.

    Raises
    ------
    ValueError
        If ``y``/``p``/``grades`` are not equal-length 1-D arrays, ``level``
        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, level=0.9)
    >>> bands["A"][0]
    0.0
    >>> bands["A"][1] < bands["B"][1]
    True
    """
    y_arr = _validate_binary_y(y)
    p_arr = validate_scores(p, name="p")
    if len(p_arr) != len(y_arr):
        raise ValueError("y and p must have equal length")
    from .grade import _resolve_grades

    g_arr, default_order = _resolve_grades(grades, p_arr)
    if g_arr.ndim != 1 or len(g_arr) != len(y_arr):
        raise ValueError("grades must be a 1-D array matching y and p in length")
    if not 0.0 < level < 1.0:
        raise ValueError("level must lie in (0, 1)")

    g_str = np.array([str(g) for g in g_arr])
    unique_labels = tuple(sorted(str(u) for u in np.unique(g_str)))

    if order is None and default_order is not None:
        order_t = default_order
    elif order is None:
        mean_p = {lab: float(np.mean(p_arr[g_str == lab])) for lab in unique_labels}
        order_t = tuple(sorted(unique_labels, key=lambda lab: mean_p[lab]))
    else:
        order_t = tuple(str(g) for g in np.asarray(order).reshape(-1))
        if tuple(sorted(order_t)) != unique_labels:
            raise ValueError(
                f"order {order_t} does not match the unique grade labels {unique_labels}"
            )

    n = np.empty(len(order_t), dtype=np.float64)
    k = np.empty(len(order_t), dtype=np.float64)
    for i, label in enumerate(order_t):
        mask = g_str == label
        n[i] = float(np.sum(mask))
        k[i] = float(np.sum(y_arr[mask]))

    hi_raw = np.array([beta_ppf(level, k[i] + 0.5, n[i] - k[i] + 0.5) for i in range(len(order_t))])
    hi = pava(hi_raw, n).fitted
    if np.any(hi != hi_raw):
        warnings.warn(
            "jeffreys_upper_bands: PAVA adjusted the upper-band sequence to keep "
            "it non-decreasing best to worst",
            UserWarning,
            stacklevel=2,
        )

    lo = np.concatenate(([0.0], hi[:-1]))
    return {str(label): (float(lo[i]), float(hi[i])) for i, label in enumerate(order_t)}

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 smallest PD under which observing at most d*_i defaults in n*_i obligors 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 
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
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 smallest PD
    under which observing at most ``d*_i`` defaults in ``n*_i`` obligors 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")
    if not 0.0 < confidence < 1.0:
        raise ValueError("confidence must lie in (0, 1)")
    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)")

    pd_upper = np.empty(k, dtype=np.float64)
    for i in range(k):
        ns, ds = n_pooled[i], d_pooled[i]
        # n_pooled - d_pooled cannot go negative here: d <= n was validated
        # per grade above, and both are cumulative sums over the same
        # (worse-grades-first) order, so the inequality is preserved term
        # by term -- beta_ppf's second shape parameter stays > 0 (b == 0
        # only in the d* == n* case handled explicitly).
        pd_upper[i] = 1.0 if ds == ns else beta_ppf(confidence, ds + 1.0, ns - ds)

    # Pooled sets are nested (grade i's pooled set contains grade i + 1's),
    # so pd_upper comes out non-decreasing already whenever the observed
    # per-grade default rates respect rating order. A noisy grade whose own
    # rate exceeds the worse-grade pool it joins can still produce a real
    # local dip; the most-prudent reading of such a dip is to raise the
    # worse grade's bound up to the better grade's, never to lower the
    # better grade's bound to match -- so the touch-up is the cumulative
    # maximum best to worst (a prudent hull), which never reduces any
    # grade's bound.
    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 non-negative weights, same length as y; per-grade counts become weighted sums.

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 
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
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 non-negative weights, same length as ``y``; per-grade
        counts become weighted sums.

    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.])
    """
    from .grade import _resolve_grades

    y_arr = _validate_binary_y(y)
    if hasattr(grades, "assign") and hasattr(grades, "names"):
        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")
        g_str, default_order = _resolve_grades(grades, p_arr)
    else:
        if order is None:
            raise ValueError("order is required when grades is a label array")
        g_str, default_order = _resolve_grades(grades, np.empty(0))
    if g_str.ndim != 1 or len(g_str) != len(y_arr):
        raise ValueError("grades and y must be 1-D arrays of equal length")
    w_arr = validate_weights(sample_weight, len(y_arr))

    unique_labels = tuple(sorted(np.unique(g_str)))
    if order is None:
        order_t = default_order if default_order is not None else ()
    else:
        order_t = tuple(str(g) for g in np.asarray(order).reshape(-1))
        if tuple(sorted(order_t)) != unique_labels:
            raise ValueError(
                f"order {order_t} does not match the unique grade labels {unique_labels}"
            )

    n = np.empty(len(order_t), dtype=np.float64)
    d = np.empty(len(order_t), dtype=np.float64)
    for i, label in enumerate(order_t):
        mask = g_str == label
        n[i] = float(np.sum(w_arr[mask]))
        d[i] = float(np.sum(w_arr[mask] * y_arr[mask]))

    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 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
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, _ = _prep(y, p, None)
    n = len(p_arr)
    if n < 2:
        raise ValueError(f"skce needs at least 2 observations, got {n}")
    s = _kernel_input(p_arr, scale)
    bw = _resolve_bandwidth(s, bandwidth)
    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() / n**2)
    return float((h.sum() - np.trace(h)) / (n * (n - 1)))

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 
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
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}")
    if n_boot < 1:
        raise ValueError(f"n_boot must be at least 1, got {n_boot}")
    y_arr, p_arr, _ = _prep(y, p, None)
    n = len(p_arr)
    if n < 4:
        raise ValueError(f"skce_test needs at least 4 observations, got {n}")
    s = _kernel_input(p_arr, scale)
    bw = _resolve_bandwidth(s, bandwidth)

    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 = float(1.0 - norm_cdf(np.array([z]))[0])
        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 = float((h.sum() - np.trace(h)) / (n * (n - 1)))
    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)
    t_b = (quad - counts @ np.diag(h_tilde)) / n
    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 non-negative 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 
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
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 non-negative weights, same length as ``y``.

    Returns
    -------
    GuardrailReport
        Slope, intercept, and Spiegelhalter-p values with pass/fail flags.
    """
    slope = calibration_slope(y, p, sample_weight=sample_weight)
    intercept = calibration_intercept(y, p, sample_weight=sample_weight)
    sp = spiegelhalter_z(y, p, sample_weight=sample_weight)
    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 non-negative 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 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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 non-negative weights, same length as ``y``.

    Returns
    -------
    float
        Fitted intercept in log-odds units.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    z = logit(p_arr)
    res = irls_logistic(np.ones((len(z), 1)), y_arr, w=w, offset=z)
    return float(res.beta[0])

calibration_slope 

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

Cox calibration slope.

< 1 means overfitting/overconfident spread, > 1 underfitting.

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

sample_weight

Optional non-negative weights, same length as y.

TYPE: (array_like or None, keyword - only) DEFAULT: None

RETURNS DESCRIPTION
float

Fitted slope on the logit scale.
Source code in src/probcal/metrics/regression.py 
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
def calibration_slope(y: object, p: object, *, sample_weight: object = None) -> float:
    """Cox calibration slope.

    ``< 1`` means overfitting/overconfident spread, ``> 1`` underfitting.

    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 non-negative weights, same length as ``y``.

    Returns
    -------
    float
        Fitted slope on the logit scale.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    z = logit(p_arr)
    X = np.column_stack([np.ones_like(z), z])
    res = irls_logistic(X, y_arr, w=w)
    return float(res.beta[1])

calibration_test 

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

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

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

sample_weight

Optional non-negative 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 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def calibration_test(
    y: object, p: object, *, sample_weight: object = None
) -> CalibrationTestResult:
    """Likelihood-ratio test of joint calibration (alpha, beta) = (0, 1).

    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 non-negative 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)
    z = logit(p_arr)
    X = np.column_stack([np.ones_like(z), z])
    fit = irls_logistic(X, y_arr, w=w)
    alpha, beta = float(fit.beta[0]), float(fit.beta[1])

    def _ll(prob: np.ndarray) -> float:
        prob = np.clip(prob, 1e-12, 1.0 - 1e-12)
        return float(np.sum(w * (y_arr * np.log(prob) + (1.0 - y_arr) * np.log1p(-prob))))

    ll_fit = _ll(expit(X @ fit.beta))
    ll_null = _ll(p_arr)
    lr = max(2.0 * (ll_fit - ll_null), 0.0)
    p_value = 1.0 - float(gammainc_lower(1.0, lr / 2.0))  # chi-square, df = 2
    return CalibrationTestResult(statistic=lr, p_value=p_value, alpha=alpha, beta=beta)

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 non-negative 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 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
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 non-negative 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 non-negative 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 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
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 non-negative 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 non-negative 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 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
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 non-negative 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 non-negative 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 
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
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 non-negative 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), 1e-12, 1.0 - 1e-12)
    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 exactly (a per-observation calculation: integrating a single observation's elementary score over theta in [0, 1] gives p**2 / 2 when y=0 and (1 - p)**2 / 2 when y=1, whose doubled weighted mean is exactly E[(1-y)*p**2 + y*(1-p)**2] == E[(p - y)**2], the Brier score). S_theta is piecewise linear in theta on each open interval between consecutive breakpoints u = sorted(unique(p) | {0, 1}) and jumps only exactly at those breakpoints (the observation at that p crosses sides), so: evaluated at the midpoint of each interval — an interior point, never a breakpoint, so the jump ambiguity never arises — the midpoint rule is exact for a linear function on an interval, and 2 * sum(np.diff(u) * murphy_curve(y, p, thresholds=mid).score) (mid = (u[1:] + u[:-1]) / 2) reproduces the Brier score to machine precision; evaluated directly at the breakpoints instead (e.g. via plain np.trapezoid over u itself), each sampled value is a one-sided limit of the jump there, so that discretization converges to the Brier identity only at a rate that shrinks with the sample size — already far inside the default grid's 1e-3 budget at realistic n, but not exact at any finite n. Isotonic (PAV) recalibration of p never increases the score at any threshold (Ehm et al., 2016), so the two curves' relative position diagnoses the value of recalibration without collapsing to one scalar. Computed by sorting p once and accumulating weighted class-conditional sums via searchsorted — O(n log n + T log n) for T (arbitrary, not necessarily uniformly spaced) thresholds, never the naive O(n * T) mask.

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 non-negative 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 
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
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
    exactly (a per-observation calculation: integrating a single
    observation's elementary score over ``theta in [0, 1]`` gives
    ``p**2 / 2`` when ``y=0`` and ``(1 - p)**2 / 2`` when ``y=1``, whose
    doubled weighted mean is exactly ``E[(1-y)*p**2 + y*(1-p)**2] ==
    E[(p - y)**2]``, the Brier score). ``S_theta`` is piecewise linear in
    ``theta`` on each *open* interval between consecutive breakpoints
    ``u = sorted(unique(p) | {0, 1})`` and jumps only exactly at those
    breakpoints (the observation at that ``p`` crosses sides), so:
    evaluated at the *midpoint* of each interval — an interior point,
    never a breakpoint, so the jump ambiguity never arises — the midpoint
    rule is exact for a linear function on an interval, and
    ``2 * sum(np.diff(u) * murphy_curve(y, p, thresholds=mid).score)``
    (``mid = (u[1:] + u[:-1]) / 2``) reproduces the Brier score to machine
    precision; evaluated directly *at* the breakpoints instead (e.g. via
    plain ``np.trapezoid`` over ``u`` itself), each sampled value is a
    one-sided limit of the jump there, so that discretization converges to
    the Brier identity only at a rate that shrinks with the sample size —
    already far inside the default grid's 1e-3 budget at realistic ``n``,
    but not exact at any finite ``n``. Isotonic (PAV) recalibration of
    ``p`` never increases the score at any threshold (Ehm et al., 2016),
    so the two curves' relative position diagnoses the value of
    recalibration without collapsing to one scalar. Computed by sorting
    ``p`` once and accumulating weighted class-conditional sums via
    ``searchsorted`` — ``O(n log n + T log n)`` for ``T`` (arbitrary,
    not necessarily uniformly spaced) thresholds, never the naive
    ``O(n * T)`` mask.

    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 non-negative 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.

TYPE: (bool, keyword - only) DEFAULT: False

sample_weight

Optional non-negative 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 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
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.
    sample_weight : array_like or None, keyword-only
        Optional non-negative weights, same length as ``y``.

    Returns
    -------
    MurphyDecomposition
        Reliability, resolution, and uncertainty terms.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    from .binned import _bin_index

    idx, m = _bin_index(p_arr, n_bins, strategy)
    w_tot = float(w.sum())
    y_bar = float(np.average(y_arr, weights=w))
    rel = res = 0.0
    for b in range(m):
        mask = idx == b
        if not np.any(mask):
            continue
        wb = float(w[mask].sum())
        pb = float(np.average(p_arr[mask], weights=w[mask]))
        yb = float(np.average(y_arr[mask], weights=w[mask]))
        nb = int(np.sum(mask))
        rel_term = (pb - yb) ** 2
        res_term = (yb - y_bar) ** 2
        if bias_corrected and nb > 1:
            var_yb = yb * (1.0 - yb) / (nb - 1)
            rel_term = max(rel_term - var_yb, 0.0)
            res_term = max(res_term - var_yb, 0.0)
        rel += (wb / w_tot) * rel_term
        res += (wb / w_tot) * 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 recovers 0.1.2 values exactly. The LOESS distances are always unweighted; sample_weight, when given and not uniform, weights only the quantile step (see :func: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 non-negative weights; used only for the quantile step.

TYPE: (array_like or None, keyword - only) DEFAULT: None

grid_size

LOESS evaluation grid size; None recovers 0.1.2 values exactly.

TYPE: (int or None, keyword - only) DEFAULT: 512

RETURNS DESCRIPTION
float

Median of the LOESS distances.
Source code in src/probcal/metrics/smooth.py 
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
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`` recovers 0.1.2 values exactly. The LOESS distances are
    always unweighted; ``sample_weight``, when given and
    not uniform, weights only the quantile step (see
    :func:`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 non-negative weights; used only for the quantile step.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` recovers 0.1.2 values exactly.

    Returns
    -------
    float
        Median of the LOESS distances.
    """
    d = _ici_distances(y, p, frac, grid_size)
    return _ici_quantile(d, 0.5, y, p, sample_weight)

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 recovers 0.1.2 values exactly. The LOESS distances are always unweighted; sample_weight, when given and not uniform, weights only the quantile step (see :func: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 non-negative weights; used only for the quantile step.

TYPE: (array_like or None, keyword - only) DEFAULT: None

grid_size

LOESS evaluation grid size; None recovers 0.1.2 values exactly.

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 
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
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`` recovers 0.1.2 values exactly. The LOESS distances are
    always unweighted; ``sample_weight``, when given and
    not uniform, weights only the quantile step (see
    :func:`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 non-negative weights; used only for the quantile step.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` recovers 0.1.2 values exactly.

    Returns
    -------
    float
        90th percentile of the LOESS distances.
    """
    d = _ici_distances(y, p, frac, grid_size)
    return _ici_quantile(d, 0.9, y, p, sample_weight)

ecce 

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

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

Sort by prediction and walk the cumulative sum of weighted residuals; under calibration the walk hovers near zero, and drift localizes miscalibration without any smoothing parameter.

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

sample_weight

Optional non-negative weights, same length as y.

TYPE: (array_like or None, keyword - only) DEFAULT: None

presorted

Declare that p is already sorted ascending, so the internal argsort can be skipped. Purely a throughput switch for callers that already hold a sorted copy (evaluate's bootstrap sorts each replicate once and shares that order across metrics); the result is unchanged when the declaration holds and meaningless when it does not, and nothing checks it.

TYPE: (bool, keyword - only) DEFAULT: False

RETURNS DESCRIPTION
EcceResult

Max and mean absolute cumulative deviation.
Source code in src/probcal/metrics/smooth.py 
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
def ecce(
    y: object, p: object, *, sample_weight: object = None, presorted: bool = False
) -> EcceResult:
    """Cumulative-deviation calibration error (Arrieta-Ibarra et al., 2022).

    Sort by prediction and walk the cumulative sum of weighted residuals;
    under calibration the walk hovers near zero, and drift localizes
    miscalibration without any smoothing parameter.

    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 non-negative weights, same length as ``y``.
    presorted : bool, keyword-only
        Declare that ``p`` is already sorted ascending, so the internal
        ``argsort`` can be skipped. Purely a throughput switch for callers that
        already hold a sorted copy (``evaluate``'s bootstrap sorts each
        replicate once and shares that order across metrics); the result is
        unchanged when the declaration holds and meaningless when it does not,
        and nothing checks it.

    Returns
    -------
    EcceResult
        Max and mean absolute cumulative deviation.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    if presorted:
        c = np.cumsum(w * (y_arr - p_arr)) / w.sum()
    else:
        order = np.argsort(p_arr, kind="stable")
        c = np.cumsum(w[order] * (y_arr[order] - p_arr[order])) / w.sum()
    return EcceResult(stat_max=float(np.max(np.abs(c))), stat_mean=float(np.mean(np.abs(c))))

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 recovers 0.1.2 values exactly.

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

Accepted for signature parity with the other ICI-family metrics but not used: the maximum is a weight-independent order statistic.

TYPE: (array_like or None, keyword - only) DEFAULT: None

grid_size

LOESS evaluation grid size; None recovers 0.1.2 values exactly.

TYPE: (int or None, keyword - only) DEFAULT: 512

RETURNS DESCRIPTION
float

Maximum of the LOESS distances.
Source code in src/probcal/metrics/smooth.py 
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
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`` recovers 0.1.2 values exactly.

    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
        Accepted for signature parity with the other ICI-family metrics but
        not used: the maximum is a weight-independent order statistic.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` recovers 0.1.2 values exactly.

    Returns
    -------
    float
        Maximum of the LOESS distances.
    """
    return float(np.max(_ici_distances(y, p, frac, 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 recovers 0.1.2 values exactly.

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 non-negative 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 recovers 0.1.2 values exactly.

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 
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
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`` recovers 0.1.2 values exactly.

    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 non-negative 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`` recovers 0.1.2 values exactly.

    Returns
    -------
    float
        Weighted mean absolute LOESS-to-prediction distance.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    c = loess(p_arr, y_arr, frac=frac, grid_size=grid_size)
    return float(np.average(np.abs(c - p_arr), weights=w))

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 pre-aggregates the weighted residual measure onto a regular grid over the logit range before solving the fixed point; the binned measure is then evaluated in closed form on its own lattice by direct Gaussian convolution, at a cost independent of n and of sigma. The lattice path engages for every call with a non-degenerate logit range (0.1.3 engaged it only for n > bins, leaving typical calibration-set sizes on the exact O(n)-per-step path — the "size cliff"). With bins=None, or a degenerate range (t.max() == t.min()), the exact 0.1.2 computation runs bit-for-bit. Otherwise, if the found sigma is smaller than 8 bin widths (the kernel would be under-resolved by the bins), the solve is repeated once on an adaptively refined binning (bins <- ceil(range / (sigma/8))); the exact computation is used only when that refinement is infeasible (refined bin count above 2**20) or still under-resolved — reachable for near-perfectly-calibrated data spread over a wide logit range (e.g. extreme/clipped scores), so the worst case matches the pre-0.1.3 O(n) cost. For n <= bins the lattice value may differ from the exact grid at the ~1e-4 level on typical portfolios (measured <= 2.4e-4 on make_pd_portfolio); on wide clipped-logit-range data the gap can be much larger because the exact path's fixed 257-point grid under-resolves small-sigma kernels there — in that regime the lattice value is the better one (>= 8 samples per sigma). bins=None recovers the old values.

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

sample_weight

Optional non-negative 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 
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
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`` pre-aggregates the weighted residual measure onto a regular grid
    over the logit range before solving the fixed point; the binned measure
    is then evaluated in closed form on its own lattice by direct Gaussian
    convolution, at a cost independent of n and of sigma. The lattice path
    engages for every call with a non-degenerate logit range
    (0.1.3 engaged it only for ``n > bins``, leaving typical calibration-set
    sizes on the exact O(n)-per-step path — the "size cliff").
    With ``bins=None``, or a degenerate range
    (``t.max() == t.min()``), the exact 0.1.2 computation runs bit-for-bit.
    Otherwise, if the found ``sigma`` is smaller than 8 bin widths (the
    kernel would be under-resolved by the bins), the solve is repeated once
    on an adaptively refined binning (``bins <- ceil(range / (sigma/8))``);
    the exact computation is used only when that refinement is infeasible
    (refined bin count above ``2**20``) or still under-resolved — reachable
    for near-perfectly-calibrated data spread over a wide logit range (e.g.
    extreme/clipped scores), so the worst case matches the pre-0.1.3 O(n)
    cost. For ``n <= bins`` the lattice value may differ from the exact
    grid at the ~1e-4 level on typical portfolios (measured <= 2.4e-4 on
    ``make_pd_portfolio``); on wide clipped-logit-range data the gap can be
    much larger because the exact path's fixed 257-point grid under-resolves
    small-sigma kernels there — in that regime the lattice value is the
    better one (>= 8 samples per sigma). ``bins=None`` recovers the old
    values.

    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 non-negative 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.

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

sample_weight

Optional non-negative 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 
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
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.

    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 non-negative weights, same length as ``y``.

    Returns
    -------
    SpiegelhalterResult
        Z statistic and two-sided p-value.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    num = float(np.sum(w * (y_arr - p_arr) * (1.0 - 2.0 * p_arr)))
    var = float(np.sum(w**2 * (1.0 - 2.0 * p_arr) ** 2 * p_arr * (1.0 - p_arr)))
    z = num / math.sqrt(var)
    p_value = float(2.0 * (1.0 - norm_cdf(np.array([abs(z)]))[0]))
    return SpiegelhalterResult(z=z, p_value=p_value)

evaluate 

evaluate(y: object, p: object, *, sample_weight: object = None, n_boot: int = 1000, seed: 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, seed: 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, seed: 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

Observation weights (resampled together with the observations).

TYPE: array_like or None DEFAULT: None

n_boot

Case-resampling bootstrap replicates (percentile CIs at 2.5/97.5).

TYPE: int DEFAULT: 1000

seed

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 each label is stringified and a separate report is computed per sorted label — group i (in sorted-label order) is evaluated with seed + 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 seed 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 by is given with a length that does not match y; 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

Cost model per replicate: scores and regression metrics and ECCE are O(n); binned ECEs are O(n log n); the ICI family (ici/e50/e90/emax) shares one LOESS fit at O(grid_size * frac * n); smECE is O(n + 257 * bins) per bisection step. All of the above are paid n_boot times — for n > 1e6, reduce n_boot or pass a metrics= subset. With by given, the whole cost model above is paid once per group plus once for the pooled report.

Each replicate is sorted by prediction once and that order is shared: the LOESS fit and ECCE skip their own sorts, ece/ece_debiased/ mce share one 15-bin equal-mass binning pass, ece_sweep's ~99-candidate scan reads per-bin sums off prefix-sum differences at searchsorted cut positions, and the LOESS anchor fits are solved in vectorized blocks rather than one Python iteration per anchor. The reported point estimates are computed on the unsorted, scalar path and are bit-for-bit what 0.2.x produced; only the replicates take the fast path, whose reordered sums move percentile CI bounds in their last bits (measured <= 4e-11 relative) and whose tricube weight cubes by multiplication rather than ** 3 (<= 2.3e-16 relative on a well-conditioned window; on a rank-deficient one the abs(det) < _FPMIN guard in the local-linear solve can select a different branch than the scalar loop, where the swy / sw branch is the well-defined answer — see _math._loess_fit_sorted_vec. Anchors are data quantiles, so this has not been observed to reach a reported value). On the dev host at n=1e4 a full-catalog replicate costs 0.089s — 58% of it the ICI family's LOESS fit, 27% the ece_sweep scan, 10% intercept/slope, 0.5% the whole binned ECE family — and the full run (n_boot=1000) takes 87s against 304s in 0.2.x. Excluding the ICI family via metrics= remains the single largest lever on cost. See docs/concepts/metrics.md for the measured table.

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/__init__.py 
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
def evaluate(
    y: object,
    p: object,
    *,
    sample_weight: object = None,
    n_boot: int = 1000,
    seed: 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
        Observation weights (resampled together with the observations).
    n_boot : int
        Case-resampling bootstrap replicates (percentile CIs at 2.5/97.5).
    seed : 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
        each label is stringified and a separate report is computed per
        sorted label — group ``i`` (in sorted-label order) is evaluated
        with ``seed + 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 ``seed``
        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
        ``by`` is given with a length that does not match ``y``; 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
    -----
    Cost model per replicate: scores and regression metrics and ECCE are
    O(n); binned ECEs are O(n log n); the ICI family (ici/e50/e90/emax)
    shares one LOESS fit at O(grid_size * frac * n); smECE is
    O(n + 257 * bins) per bisection step. All of the above are paid
    ``n_boot`` times — for n > 1e6, reduce ``n_boot`` or pass a ``metrics=``
    subset. With ``by`` given, the whole cost model above is paid once per
    group plus once for the pooled report.

    Each replicate is sorted by prediction once and that order is shared:
    the LOESS fit and ECCE skip their own sorts, ``ece``/``ece_debiased``/
    ``mce`` share one 15-bin equal-mass binning pass, ``ece_sweep``'s
    ~99-candidate scan reads per-bin sums off prefix-sum differences at
    ``searchsorted`` cut positions, and the LOESS anchor fits are solved in
    vectorized blocks rather than one Python iteration per anchor. The
    reported *point* estimates are computed on the unsorted, scalar path and
    are bit-for-bit what 0.2.x produced; only the replicates take the fast
    path, whose reordered sums move percentile CI bounds in their last bits
    (measured <= 4e-11 relative) and whose tricube weight cubes by
    multiplication rather than ``** 3`` (<= 2.3e-16 relative on a
    well-conditioned window; on a rank-deficient one the
    ``abs(det) < _FPMIN`` guard in the local-linear solve can select a
    different branch than the scalar loop, where the ``swy / sw`` branch is
    the well-defined answer — see ``_math._loess_fit_sorted_vec``. Anchors
    are data quantiles, so this has not been observed to reach a reported
    value). On the dev
    host at n=1e4 a full-catalog replicate costs 0.089s — 58% of it the ICI
    family's LOESS fit, 27% the ``ece_sweep`` scan, 10% intercept/slope,
    0.5% the whole binned ECE family — and the full run
    (``n_boot=1000``) takes 87s against 304s in 0.2.x. Excluding the ICI
    family via ``metrics=`` remains the single largest lever on cost. See
    ``docs/concepts/metrics.md`` for the measured table.

    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
    """
    if by is not None:
        return _evaluate_grouped(
            y,
            p,
            by,
            sample_weight=sample_weight,
            n_boot=n_boot,
            seed=seed,
            metrics=metrics,
            stratify=stratify,
        )

    from .scores import _prep

    if metrics is not None:
        unknown = sorted(set(metrics) - set(_METRIC_CATALOG))
        if unknown:
            raise ValueError(
                f"unknown metric names {unknown}; valid names: {list(_METRIC_CATALOG)}"
            )
        names = tuple(k for k in _METRIC_CATALOG if k in set(metrics))
    else:
        names = _METRIC_CATALOG

    y_arr, p_arr, w_arr = _prep(y, p, sample_weight)
    point = _point_metrics(y_arr, p_arr, w_arr, names)
    values = np.array([point[k] for k in names])

    rng = np.random.default_rng(seed)
    n = len(y_arr)
    # _prep -> validate_binary_y already rejects single-class y unconditionally
    # (both idx0 and idx1 are therefore guaranteed non-empty here).
    idx0 = np.flatnonzero(y_arr == 0)
    idx1 = np.flatnonzero(y_arr == 1)

    boot = np.empty((n_boot, len(names)))
    for b in range(n_boot):
        if stratify:
            idx = np.concatenate(
                [
                    idx0[rng.integers(0, len(idx0), len(idx0))],
                    idx1[rng.integers(0, len(idx1), len(idx1))],
                ]
            )
        else:
            for _attempt in range(100):
                idx = rng.integers(0, n, n)
                if y_arr[idx].min() != y_arr[idx].max():
                    break
            else:
                raise RuntimeError(
                    "100 consecutive degenerate (single-class) bootstrap draws; "
                    "pass stratify=True or supply more data"
                )
        # One stable sort per replicate, shared by every metric that would
        # otherwise sort (or re-bin) on its own; see ``_point_metrics``.
        idx = idx[np.argsort(p_arr[idx], kind="stable")]
        yb, pb, wb = y_arr[idx], p_arr[idx], w_arr[idx]
        pm = _point_metrics(yb, pb, wb, names, presorted=True)
        boot[b] = [pm[k] for k in names]
    ci_low = np.percentile(boot, 2.5, axis=0)
    ci_high = np.percentile(boot, 97.5, axis=0)
    return MetricReport(names=names, values=values, ci_low=ci_low, ci_high=ci_high)

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 recovers 0.1.2 values exactly.

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

sample_weight

Optional non-negative 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 recovers 0.1.2 values exactly.

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/__init__.py 
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
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`` recovers
    0.1.2 values exactly.

    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 non-negative weights, same length as ``y``.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size for the ICI/E90 terms; ``None`` recovers
        0.1.2 values exactly.

    Returns
    -------
    ReliabilitySummary
        Stats-box fields for the annotated reliability diagram.
    """
    from .scores import _prep

    y_arr, p_arr, w = _prep(y, p, sample_weight)
    wq = None if sample_weight is None else w
    return ReliabilitySummary(
        n=len(y_arr),
        events=int(y_arr.sum()),
        intercept=calibration_intercept(y_arr, p_arr, sample_weight=wq),
        slope=calibration_slope(y_arr, p_arr, sample_weight=wq),
        ici=ici(y_arr, p_arr, sample_weight=wq, grid_size=grid_size),
        e90=e90(y_arr, p_arr, sample_weight=wq, grid_size=grid_size),
        spiegelhalter_p=spiegelhalter_z(y_arr, p_arr, sample_weight=wq).p_value,
    )

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 non-negative 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 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
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 non-negative 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 non-negative 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 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
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 non-negative 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 non-negative 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 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
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 non-negative 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.

TYPE: (bool, keyword - only) DEFAULT: False

sample_weight

Optional non-negative 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 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
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.
    sample_weight : array_like or None, keyword-only
        Optional non-negative weights, same length as ``y``.

    Returns
    -------
    MurphyDecomposition
        Reliability, resolution, and uncertainty terms.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    from .binned import _bin_index

    idx, m = _bin_index(p_arr, n_bins, strategy)
    w_tot = float(w.sum())
    y_bar = float(np.average(y_arr, weights=w))
    rel = res = 0.0
    for b in range(m):
        mask = idx == b
        if not np.any(mask):
            continue
        wb = float(w[mask].sum())
        pb = float(np.average(p_arr[mask], weights=w[mask]))
        yb = float(np.average(y_arr[mask], weights=w[mask]))
        nb = int(np.sum(mask))
        rel_term = (pb - yb) ** 2
        res_term = (yb - y_bar) ** 2
        if bias_corrected and nb > 1:
            var_yb = yb * (1.0 - yb) / (nb - 1)
            rel_term = max(rel_term - var_yb, 0.0)
            res_term = max(res_term - var_yb, 0.0)
        rel += (wb / w_tot) * rel_term
        res += (wb / w_tot) * 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 exactly (a per-observation calculation: integrating a single observation's elementary score over theta in [0, 1] gives p**2 / 2 when y=0 and (1 - p)**2 / 2 when y=1, whose doubled weighted mean is exactly E[(1-y)*p**2 + y*(1-p)**2] == E[(p - y)**2], the Brier score). S_theta is piecewise linear in theta on each open interval between consecutive breakpoints u = sorted(unique(p) | {0, 1}) and jumps only exactly at those breakpoints (the observation at that p crosses sides), so: evaluated at the midpoint of each interval — an interior point, never a breakpoint, so the jump ambiguity never arises — the midpoint rule is exact for a linear function on an interval, and 2 * sum(np.diff(u) * murphy_curve(y, p, thresholds=mid).score) (mid = (u[1:] + u[:-1]) / 2) reproduces the Brier score to machine precision; evaluated directly at the breakpoints instead (e.g. via plain np.trapezoid over u itself), each sampled value is a one-sided limit of the jump there, so that discretization converges to the Brier identity only at a rate that shrinks with the sample size — already far inside the default grid's 1e-3 budget at realistic n, but not exact at any finite n. Isotonic (PAV) recalibration of p never increases the score at any threshold (Ehm et al., 2016), so the two curves' relative position diagnoses the value of recalibration without collapsing to one scalar. Computed by sorting p once and accumulating weighted class-conditional sums via searchsorted — O(n log n + T log n) for T (arbitrary, not necessarily uniformly spaced) thresholds, never the naive O(n * T) mask.

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 non-negative 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 
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
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
    exactly (a per-observation calculation: integrating a single
    observation's elementary score over ``theta in [0, 1]`` gives
    ``p**2 / 2`` when ``y=0`` and ``(1 - p)**2 / 2`` when ``y=1``, whose
    doubled weighted mean is exactly ``E[(1-y)*p**2 + y*(1-p)**2] ==
    E[(p - y)**2]``, the Brier score). ``S_theta`` is piecewise linear in
    ``theta`` on each *open* interval between consecutive breakpoints
    ``u = sorted(unique(p) | {0, 1})`` and jumps only exactly at those
    breakpoints (the observation at that ``p`` crosses sides), so:
    evaluated at the *midpoint* of each interval — an interior point,
    never a breakpoint, so the jump ambiguity never arises — the midpoint
    rule is exact for a linear function on an interval, and
    ``2 * sum(np.diff(u) * murphy_curve(y, p, thresholds=mid).score)``
    (``mid = (u[1:] + u[:-1]) / 2``) reproduces the Brier score to machine
    precision; evaluated directly *at* the breakpoints instead (e.g. via
    plain ``np.trapezoid`` over ``u`` itself), each sampled value is a
    one-sided limit of the jump there, so that discretization converges to
    the Brier identity only at a rate that shrinks with the sample size —
    already far inside the default grid's 1e-3 budget at realistic ``n``,
    but not exact at any finite ``n``. Isotonic (PAV) recalibration of
    ``p`` never increases the score at any threshold (Ehm et al., 2016),
    so the two curves' relative position diagnoses the value of
    recalibration without collapsing to one scalar. Computed by sorting
    ``p`` once and accumulating weighted class-conditional sums via
    ``searchsorted`` — ``O(n log n + T log n)`` for ``T`` (arbitrary,
    not necessarily uniformly spaced) thresholds, never the naive
    ``O(n * T)`` mask.

    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 non-negative 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 non-negative 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 
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
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 non-negative 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), 1e-12, 1.0 - 1e-12)
    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.

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 (used groups minus 2, floored at 1).

TYPE: int

p_value

Upper-tail p-value of the chi-square statistic.

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.

TYPE: (int, keyword - only) DEFAULT: 15

strategy

"mass" (equal-count, default) 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 non-negative 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 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
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.
    strategy : {"mass", "width"}, keyword-only
        ``"mass"`` (equal-count, default) 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 non-negative 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)
    shares, gaps, _, _ = _bin_gaps(y_arr, p_arr, w, n_bins, strategy)
    return _ece_from_gaps(shares, gaps, 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).

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 non-negative 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 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
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).

    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 non-negative weights, same length as ``y``.

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

ece_sweep 

ece_sweep(y: object, p: object, *, norm: str = 'l1', sample_weight: object = None) -> float

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

Uses equal-mass bins with the largest B whose bin event rates remain monotone non-decreasing (scan 2..min(n, 100)).

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

norm

Norm passed to the final :func:ece call at the selected bin count.

TYPE: (l1, l2, max) DEFAULT: "l1"

sample_weight

Optional non-negative weights, same length as y.

TYPE: (array_like or None, keyword - only) DEFAULT: None

RETURNS DESCRIPTION
float

Calibration error at the largest monotone bin count.
Source code in src/probcal/metrics/binned.py 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def ece_sweep(
    y: object,
    p: object,
    *,
    norm: str = "l1",
    sample_weight: object = None,
) -> float:
    """Monotonic-sweep calibration error (Roelofs et al., 2022).

    Uses equal-mass bins with the largest ``B`` whose bin event rates remain
    monotone non-decreasing (scan 2..min(n, 100)).

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    norm : {"l1", "l2", "max"}, keyword-only
        Norm passed to the final :func:`ece` call at the selected bin count.
    sample_weight : array_like or None, keyword-only
        Optional non-negative weights, same length as ``y``.

    Returns
    -------
    float
        Calibration error at the largest monotone bin count.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    best_b = 1
    for b in range(2, min(len(p_arr), 100) + 1):
        _, _, rates, _ = _bin_gaps(y_arr, p_arr, w, b, "mass")
        if np.all(np.diff(rates) >= 0.0):
            best_b = b
    return _ece_at_best_b(y_arr, p_arr, w, best_b, norm)

adaptive_ece 

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

Adaptive ECE: an explicit alias for equal-mass ece.

The literature uses both names for the same estimator.

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

norm

Norm passed through to :func:ece.

TYPE: (l1, l2, max) DEFAULT: "l1"

sample_weight

Optional non-negative weights, same length as y.

TYPE: (array_like or None, keyword - only) DEFAULT: None

RETURNS DESCRIPTION
float

Equal-mass binned calibration error under the chosen norm.
Source code in src/probcal/metrics/binned.py 
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
def adaptive_ece(
    y: object,
    p: object,
    *,
    n_bins: int = 15,
    norm: str = "l1",
    sample_weight: object = None,
) -> float:
    """Adaptive ECE: an explicit alias for equal-mass ``ece``.

    The literature uses both names for the same estimator.

    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.
    norm : {"l1", "l2", "max"}, keyword-only
        Norm passed through to :func:`ece`.
    sample_weight : array_like or None, keyword-only
        Optional non-negative weights, same length as ``y``.

    Returns
    -------
    float
        Equal-mass binned calibration error under the chosen norm.
    """
    return ece(y, p, n_bins=n_bins, strategy="mass", norm=norm, sample_weight=sample_weight)

hosmer_lemeshow 

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

Hosmer–Lemeshow goodness-of-fit test on g 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.

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

g

Requested number of equal-mass risk groups.

TYPE: (int, keyword - only) DEFAULT: 10

sample_weight

Optional non-negative weights, same length as y.

TYPE: (array_like or None, keyword - only) DEFAULT: None

RETURNS DESCRIPTION
HosmerLemeshowResult

Chi-square statistic, degrees of freedom, and p-value.
Source code in src/probcal/metrics/binned.py 
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
def hosmer_lemeshow(
    y: object,
    p: object,
    *,
    g: int = 10,
    sample_weight: object = None,
) -> HosmerLemeshowResult:
    """Hosmer–Lemeshow goodness-of-fit test on ``g`` 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.

    Parameters
    ----------
    y : array_like
        Binary outcomes in ``{0, 1}``.
    p : array_like
        Predicted probabilities in ``[0, 1]``.
    g : int, keyword-only
        Requested number of equal-mass risk groups.
    sample_weight : array_like or None, keyword-only
        Optional non-negative weights, same length as ``y``.

    Returns
    -------
    HosmerLemeshowResult
        Chi-square statistic, degrees of freedom, and p-value.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    idx, m = _bin_index(p_arr, g, "mass")
    stat = 0.0
    used = 0
    for b in range(m):
        mask = idx == b
        if not np.any(mask):
            continue
        nb = float(w[mask].sum())
        obs = float(np.sum(w[mask] * y_arr[mask]))
        exp = float(np.sum(w[mask] * p_arr[mask]))
        denom = exp * (1.0 - exp / nb)
        if denom > 0:
            stat += (obs - exp) ** 2 / denom
        used += 1
    df = max(used - 2, 1)
    p_value = 1.0 - float(gammainc_lower(df / 2.0, stat / 2.0))
    return HosmerLemeshowResult(statistic=float(stat), df=df, p_value=p_value)

smooth

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

Theory: docs/concepts/metrics.md.

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

Mean absolute cumulative deviation.

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.

TYPE: float

p_value

Two-sided p-value under the standard normal approximation.

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 pre-aggregates the weighted residual measure onto a regular grid over the logit range before solving the fixed point; the binned measure is then evaluated in closed form on its own lattice by direct Gaussian convolution, at a cost independent of n and of sigma. The lattice path engages for every call with a non-degenerate logit range (0.1.3 engaged it only for n > bins, leaving typical calibration-set sizes on the exact O(n)-per-step path — the "size cliff"). With bins=None, or a degenerate range (t.max() == t.min()), the exact 0.1.2 computation runs bit-for-bit. Otherwise, if the found sigma is smaller than 8 bin widths (the kernel would be under-resolved by the bins), the solve is repeated once on an adaptively refined binning (bins <- ceil(range / (sigma/8))); the exact computation is used only when that refinement is infeasible (refined bin count above 2**20) or still under-resolved — reachable for near-perfectly-calibrated data spread over a wide logit range (e.g. extreme/clipped scores), so the worst case matches the pre-0.1.3 O(n) cost. For n <= bins the lattice value may differ from the exact grid at the ~1e-4 level on typical portfolios (measured <= 2.4e-4 on make_pd_portfolio); on wide clipped-logit-range data the gap can be much larger because the exact path's fixed 257-point grid under-resolves small-sigma kernels there — in that regime the lattice value is the better one (>= 8 samples per sigma). bins=None recovers the old values.

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

sample_weight

Optional non-negative 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 
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
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`` pre-aggregates the weighted residual measure onto a regular grid
    over the logit range before solving the fixed point; the binned measure
    is then evaluated in closed form on its own lattice by direct Gaussian
    convolution, at a cost independent of n and of sigma. The lattice path
    engages for every call with a non-degenerate logit range
    (0.1.3 engaged it only for ``n > bins``, leaving typical calibration-set
    sizes on the exact O(n)-per-step path — the "size cliff").
    With ``bins=None``, or a degenerate range
    (``t.max() == t.min()``), the exact 0.1.2 computation runs bit-for-bit.
    Otherwise, if the found ``sigma`` is smaller than 8 bin widths (the
    kernel would be under-resolved by the bins), the solve is repeated once
    on an adaptively refined binning (``bins <- ceil(range / (sigma/8))``);
    the exact computation is used only when that refinement is infeasible
    (refined bin count above ``2**20``) or still under-resolved — reachable
    for near-perfectly-calibrated data spread over a wide logit range (e.g.
    extreme/clipped scores), so the worst case matches the pre-0.1.3 O(n)
    cost. For ``n <= bins`` the lattice value may differ from the exact
    grid at the ~1e-4 level on typical portfolios (measured <= 2.4e-4 on
    ``make_pd_portfolio``); on wide clipped-logit-range data the gap can be
    much larger because the exact path's fixed 257-point grid under-resolves
    small-sigma kernels there — in that regime the lattice value is the
    better one (>= 8 samples per sigma). ``bins=None`` recovers the old
    values.

    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 non-negative 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, presorted: bool = False) -> EcceResult

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

Sort by prediction and walk the cumulative sum of weighted residuals; under calibration the walk hovers near zero, and drift localizes miscalibration without any smoothing parameter.

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

sample_weight

Optional non-negative weights, same length as y.

TYPE: (array_like or None, keyword - only) DEFAULT: None

presorted

Declare that p is already sorted ascending, so the internal argsort can be skipped. Purely a throughput switch for callers that already hold a sorted copy (evaluate's bootstrap sorts each replicate once and shares that order across metrics); the result is unchanged when the declaration holds and meaningless when it does not, and nothing checks it.

TYPE: (bool, keyword - only) DEFAULT: False

RETURNS DESCRIPTION
EcceResult

Max and mean absolute cumulative deviation.
Source code in src/probcal/metrics/smooth.py 
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
def ecce(
    y: object, p: object, *, sample_weight: object = None, presorted: bool = False
) -> EcceResult:
    """Cumulative-deviation calibration error (Arrieta-Ibarra et al., 2022).

    Sort by prediction and walk the cumulative sum of weighted residuals;
    under calibration the walk hovers near zero, and drift localizes
    miscalibration without any smoothing parameter.

    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 non-negative weights, same length as ``y``.
    presorted : bool, keyword-only
        Declare that ``p`` is already sorted ascending, so the internal
        ``argsort`` can be skipped. Purely a throughput switch for callers that
        already hold a sorted copy (``evaluate``'s bootstrap sorts each
        replicate once and shares that order across metrics); the result is
        unchanged when the declaration holds and meaningless when it does not,
        and nothing checks it.

    Returns
    -------
    EcceResult
        Max and mean absolute cumulative deviation.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    if presorted:
        c = np.cumsum(w * (y_arr - p_arr)) / w.sum()
    else:
        order = np.argsort(p_arr, kind="stable")
        c = np.cumsum(w[order] * (y_arr[order] - p_arr[order])) / w.sum()
    return EcceResult(stat_max=float(np.max(np.abs(c))), stat_mean=float(np.mean(np.abs(c))))

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 recovers 0.1.2 values exactly.

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 non-negative 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 recovers 0.1.2 values exactly.

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 
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
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`` recovers 0.1.2 values exactly.

    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 non-negative 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`` recovers 0.1.2 values exactly.

    Returns
    -------
    float
        Weighted mean absolute LOESS-to-prediction distance.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    c = loess(p_arr, y_arr, frac=frac, grid_size=grid_size)
    return float(np.average(np.abs(c - p_arr), weights=w))

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 recovers 0.1.2 values exactly. The LOESS distances are always unweighted; sample_weight, when given and not uniform, weights only the quantile step (see :func: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 non-negative weights; used only for the quantile step.

TYPE: (array_like or None, keyword - only) DEFAULT: None

grid_size

LOESS evaluation grid size; None recovers 0.1.2 values exactly.

TYPE: (int or None, keyword - only) DEFAULT: 512

RETURNS DESCRIPTION
float

Median of the LOESS distances.
Source code in src/probcal/metrics/smooth.py 
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
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`` recovers 0.1.2 values exactly. The LOESS distances are
    always unweighted; ``sample_weight``, when given and
    not uniform, weights only the quantile step (see
    :func:`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 non-negative weights; used only for the quantile step.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` recovers 0.1.2 values exactly.

    Returns
    -------
    float
        Median of the LOESS distances.
    """
    d = _ici_distances(y, p, frac, grid_size)
    return _ici_quantile(d, 0.5, y, p, sample_weight)

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 recovers 0.1.2 values exactly. The LOESS distances are always unweighted; sample_weight, when given and not uniform, weights only the quantile step (see :func: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 non-negative weights; used only for the quantile step.

TYPE: (array_like or None, keyword - only) DEFAULT: None

grid_size

LOESS evaluation grid size; None recovers 0.1.2 values exactly.

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 
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
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`` recovers 0.1.2 values exactly. The LOESS distances are
    always unweighted; ``sample_weight``, when given and
    not uniform, weights only the quantile step (see
    :func:`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 non-negative weights; used only for the quantile step.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` recovers 0.1.2 values exactly.

    Returns
    -------
    float
        90th percentile of the LOESS distances.
    """
    d = _ici_distances(y, p, frac, grid_size)
    return _ici_quantile(d, 0.9, 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 recovers 0.1.2 values exactly.

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

Accepted for signature parity with the other ICI-family metrics but not used: the maximum is a weight-independent order statistic.

TYPE: (array_like or None, keyword - only) DEFAULT: None

grid_size

LOESS evaluation grid size; None recovers 0.1.2 values exactly.

TYPE: (int or None, keyword - only) DEFAULT: 512

RETURNS DESCRIPTION
float

Maximum of the LOESS distances.
Source code in src/probcal/metrics/smooth.py 
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
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`` recovers 0.1.2 values exactly.

    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
        Accepted for signature parity with the other ICI-family metrics but
        not used: the maximum is a weight-independent order statistic.
    grid_size : int or None, keyword-only
        LOESS evaluation grid size; ``None`` recovers 0.1.2 values exactly.

    Returns
    -------
    float
        Maximum of the LOESS distances.
    """
    return float(np.max(_ici_distances(y, p, frac, 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.

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

sample_weight

Optional non-negative 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 
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
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.

    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 non-negative weights, same length as ``y``.

    Returns
    -------
    SpiegelhalterResult
        Z statistic and two-sided p-value.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    num = float(np.sum(w * (y_arr - p_arr) * (1.0 - 2.0 * p_arr)))
    var = float(np.sum(w**2 * (1.0 - 2.0 * p_arr) ** 2 * p_arr * (1.0 - p_arr)))
    z = num / math.sqrt(var)
    p_value = float(2.0 * (1.0 - norm_cdf(np.array([abs(z)]))[0]))
    return SpiegelhalterResult(z=z, p_value=p_value)

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.

CalibrationTestResult dataclass 

CalibrationTestResult(statistic: float, p_value: float, alpha: float, beta: float)

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

ATTRIBUTE DESCRIPTION
statistic

Likelihood-ratio test statistic (chi-square, 2 df).

TYPE: float

p_value

Upper-tail p-value of the statistic.

TYPE: float

alpha

Fitted intercept.

TYPE: float

beta

Fitted slope.

TYPE: float

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 non-negative 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 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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 non-negative weights, same length as ``y``.

    Returns
    -------
    float
        Fitted intercept in log-odds units.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    z = logit(p_arr)
    res = irls_logistic(np.ones((len(z), 1)), y_arr, w=w, offset=z)
    return float(res.beta[0])

calibration_slope 

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

Cox calibration slope.

< 1 means overfitting/overconfident spread, > 1 underfitting.

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

sample_weight

Optional non-negative weights, same length as y.

TYPE: (array_like or None, keyword - only) DEFAULT: None

RETURNS DESCRIPTION
float

Fitted slope on the logit scale.
Source code in src/probcal/metrics/regression.py 
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
def calibration_slope(y: object, p: object, *, sample_weight: object = None) -> float:
    """Cox calibration slope.

    ``< 1`` means overfitting/overconfident spread, ``> 1`` underfitting.

    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 non-negative weights, same length as ``y``.

    Returns
    -------
    float
        Fitted slope on the logit scale.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    z = logit(p_arr)
    X = np.column_stack([np.ones_like(z), z])
    res = irls_logistic(X, y_arr, w=w)
    return float(res.beta[1])

calibration_test 

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

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

PARAMETER DESCRIPTION
y

Binary outcomes in {0, 1}.

TYPE: array_like

p

Predicted probabilities in [0, 1].

TYPE: array_like

sample_weight

Optional non-negative 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 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def calibration_test(
    y: object, p: object, *, sample_weight: object = None
) -> CalibrationTestResult:
    """Likelihood-ratio test of joint calibration (alpha, beta) = (0, 1).

    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 non-negative 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)
    z = logit(p_arr)
    X = np.column_stack([np.ones_like(z), z])
    fit = irls_logistic(X, y_arr, w=w)
    alpha, beta = float(fit.beta[0]), float(fit.beta[1])

    def _ll(prob: np.ndarray) -> float:
        prob = np.clip(prob, 1e-12, 1.0 - 1e-12)
        return float(np.sum(w * (y_arr * np.log(prob) + (1.0 - y_arr) * np.log1p(-prob))))

    ll_fit = _ll(expit(X @ fit.beta))
    ll_null = _ll(p_arr)
    lr = max(2.0 * (ll_fit - ll_null), 0.0)
    p_value = 1.0 - float(gammainc_lower(1.0, lr / 2.0))  # chi-square, df = 2
    return CalibrationTestResult(statistic=lr, p_value=p_value, alpha=alpha, beta=beta)

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 non-negative 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 
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
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 non-negative weights, same length as ``y``.

    Returns
    -------
    GuardrailReport
        Slope, intercept, and Spiegelhalter-p values with pass/fail flags.
    """
    slope = calibration_slope(y, p, sample_weight=sample_weight)
    intercept = calibration_intercept(y, p, sample_weight=sample_weight)
    sp = spiegelhalter_z(y, p, sample_weight=sample_weight)
    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)

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

ATTRIBUTE DESCRIPTION
grades

Grade labels: sorted when a label array was given, best to worst when a Masterscale was given.

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.

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

JeffreysGradeResult dataclass 

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

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

ATTRIBUTE DESCRIPTION
grades

Grade labels: sorted when a label array was given, best to worst when a Masterscale was given.

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) -> 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

RETURNS DESCRIPTION
BinomialGradeResult

Per-grade counts, p-values, traffic lights, and display intervals.
Source code in src/probcal/metrics/grade.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def binomial_grade_test(
    y: object, p: object, grades: object, *, sample_weight: 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.

    Returns
    -------
    BinomialGradeResult
        Per-grade counts, p-values, traffic lights, and display intervals.
    """
    y_arr, p_arr, _ = _prep(y, p, None)
    _check_weights(sample_weight, len(y_arr))
    g_arr, order = _resolve_grades(grades, p_arr)
    labels, n, k, pd = _per_grade(y_arr, p_arr, g_arr, order)
    p_exact = np.empty(len(labels))
    p_normal = np.empty(len(labels))
    for i in range(len(labels)):
        if k[i] == 0:
            p_exact[i] = 1.0
        else:
            p_exact[i] = float(betainc(float(k[i]), float(n[i] - k[i] + 1), pd[i]))
        se = np.sqrt(n[i] * pd[i] * (1.0 - pd[i]))
        z = (k[i] - n[i] * pd[i]) / se if se > 0 else 0.0
        p_normal[i] = float(1.0 - norm_cdf(np.array([z]))[0])
    light = tuple(_traffic_light(v) for v in p_exact)
    ci_low = np.empty(len(labels))
    ci_high = np.empty(len(labels))
    for i in range(len(labels)):
        ki, ni = int(k[i]), int(n[i])
        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))
    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) -> 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

RETURNS DESCRIPTION
JeffreysGradeResult

Per-grade counts, p-values, traffic lights, and display intervals.
Source code in src/probcal/metrics/grade.py 
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
def jeffreys_grade_test(
    y: object, p: object, grades: object, *, sample_weight: 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.

    Returns
    -------
    JeffreysGradeResult
        Per-grade counts, p-values, traffic lights, and display intervals.
    """
    y_arr, p_arr, _ = _prep(y, p, None)
    _check_weights(sample_weight, len(y_arr))
    g_arr, order = _resolve_grades(grades, p_arr)
    labels, n, k, pd = _per_grade(y_arr, p_arr, g_arr, order)
    p_value = np.empty(len(labels))
    for i in range(len(labels)):
        p_value[i] = float(betainc(k[i] + 0.5, n[i] - k[i] + 0.5, pd[i]))
    light = tuple(_traffic_light(v) for v in p_value)
    ci_low = np.empty(len(labels))
    ci_high = np.empty(len(labels))
    for i in range(len(labels)):
        a, b = k[i] + 0.5, n[i] - k[i] + 0.5
        ci_low[i] = beta_ppf(0.05, a, b)
        ci_high[i] = beta_ppf(0.95, a, b)
    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 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
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, _ = _prep(y, p, None)
    n = len(p_arr)
    if n < 2:
        raise ValueError(f"skce needs at least 2 observations, got {n}")
    s = _kernel_input(p_arr, scale)
    bw = _resolve_bandwidth(s, bandwidth)
    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() / n**2)
    return float((h.sum() - np.trace(h)) / (n * (n - 1)))

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 
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
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}")
    if n_boot < 1:
        raise ValueError(f"n_boot must be at least 1, got {n_boot}")
    y_arr, p_arr, _ = _prep(y, p, None)
    n = len(p_arr)
    if n < 4:
        raise ValueError(f"skce_test needs at least 4 observations, got {n}")
    s = _kernel_input(p_arr, scale)
    bw = _resolve_bandwidth(s, bandwidth)

    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 = float(1.0 - norm_cdf(np.array([z]))[0])
        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 = float((h.sum() - np.trace(h)) / (n * (n - 1)))
    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)
    t_b = (quad - counts @ np.diag(h_tilde)) / n
    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,
    )