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 (DECISIONS entry). Selection guidance — what may be optimized and what is report-only — is the table in docs/concepts/metrics.md.

HosmerLemeshowResult dataclass 

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

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

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.

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

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

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.

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.

LogLossDecomposition dataclass 

LogLossDecomposition(calibration: float, refinement: float)

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

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

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.

SpiegelhalterResult dataclass 

SpiegelhalterResult(z: float, p_value: float)

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

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.

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

Source code in src/probcal/metrics/binned.py 
122
123
124
125
126
127
128
129
130
131
132
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)."""
    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).

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
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).
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    shares, gaps, _, _ = _bin_gaps(y_arr, p_arr, w, n_bins, strategy)
    if norm == "l1":
        return float(np.sum(shares * gaps))
    if norm == "l2":
        return float(np.sqrt(np.sum(shares * gaps**2)))
    if norm == "max":
        return float(gaps.max())
    raise ValueError(f"norm must be 'l1', 'l2', or 'max', got {norm!r}")

ece_debiased 

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

Bias-corrected ECE: per-bin squared gaps minus the within-bin variance of the event rate, floored at zero (correction in the spirit of Bröcker 2009 / Ferro & Fricker 2012; exact estimator in the DECISIONS log).

Source code in src/probcal/metrics/binned.py 
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
def ece_debiased(
    y: object,
    p: object,
    *,
    n_bins: int = 15,
    strategy: str = "mass",
    sample_weight: object = None,
) -> float:
    """Bias-corrected ECE: per-bin squared gaps minus the within-bin variance
    of the event rate, floored at zero (correction in the spirit of Bröcker
    2009 / Ferro & Fricker 2012; exact estimator in the DECISIONS log)."""
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    shares, gaps, rates, counts = _bin_gaps(y_arr, p_arr, w, n_bins, strategy)
    corrected = np.empty_like(gaps)
    for i in range(len(gaps)):
        if counts[i] > 1:
            var_b = rates[i] * (1.0 - rates[i]) / (counts[i] - 1)
            corrected[i] = np.sqrt(max(gaps[i] ** 2 - var_b, 0.0))
        else:
            corrected[i] = gaps[i]
    return float(np.sum(shares * corrected))

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); DECISIONS entry).

Source code in src/probcal/metrics/binned.py 
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
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); DECISIONS entry).
    """
    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
    if best_b == 1:
        pb = float(np.average(p_arr, weights=w))
        yb = float(np.average(y_arr, weights=w))
        return abs(pb - yb)
    return ece(y_arr, p_arr, n_bins=best_b, strategy="mass", norm=norm, sample_weight=w)

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.

Source code in src/probcal/metrics/binned.py 
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
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.
    """
    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.

Source code in src/probcal/metrics/grade.py 
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 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
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.
    """
    y_arr, p_arr, _ = _prep(y, p, None)
    _check_weights(sample_weight, len(y_arr))
    g_arr = np.asarray(grades)
    labels, n, k, pd = _per_grade(y_arr, p_arr, g_arr)
    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.

Source code in src/probcal/metrics/grade.py 
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
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.
    """
    y_arr, p_arr, _ = _prep(y, p, None)
    _check_weights(sample_weight, len(y_arr))
    g_arr = np.asarray(grades)
    labels, n, k, pd = _per_grade(y_arr, p_arr, g_arr)
    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,
    )

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

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

Source code in src/probcal/metrics/kernel.py 
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
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.
    """
    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).

Source code in src/probcal/metrics/regression.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
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)."""
    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)).

Source code in src/probcal/metrics/regression.py 
16
17
18
19
20
21
22
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))."""
    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.

Source code in src/probcal/metrics/regression.py 
25
26
27
28
29
30
31
32
def calibration_slope(y: object, p: object, *, sample_weight: object = None) -> float:
    """Cox calibration slope: < 1 means overfitting/overconfident spread,
    > 1 underfitting."""
    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).

Source code in src/probcal/metrics/regression.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
def calibration_test(
    y: object, p: object, *, sample_weight: object = None
) -> CalibrationTestResult:
    """Likelihood-ratio test of joint calibration (alpha, beta) = (0, 1)."""
    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).

Source code in src/probcal/metrics/scores.py 
31
32
33
34
def brier_score(y: object, p: object, *, sample_weight: object = None) -> float:
    """Weighted mean squared error of the probability forecast (strictly proper)."""
    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.

Source code in src/probcal/metrics/scores.py 
37
38
39
40
41
42
43
44
45
46
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.
    """
    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).

Source code in src/probcal/metrics/scores.py 
23
24
25
26
27
28
def log_loss(y: object, p: object, *, sample_weight: object = None) -> float:
    """Weighted mean negative log-likelihood (strictly proper; the default
    selection criterion)."""
    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.

Source code in src/probcal/metrics/scores.py 
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
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``.
    """
    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_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.

Source code in src/probcal/metrics/scores.py 
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
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.
    """
    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) -> float

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

Source code in src/probcal/metrics/smooth.py 
89
90
91
def e50(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None) -> float:
    """Median of the |LOESS(y|p) - p| distances."""
    return float(np.quantile(_ici_distances(y, p, frac), 0.5))

e90 

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

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

Source code in src/probcal/metrics/smooth.py 
94
95
96
def e90(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None) -> float:
    """90th percentile of the |LOESS(y|p) - p| distances."""
    return float(np.quantile(_ici_distances(y, p, frac), 0.9))

ecce 

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

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

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

Source code in src/probcal/metrics/smooth.py 
62
63
64
65
66
67
68
69
70
71
72
def ecce(y: object, p: object, *, sample_weight: object = None) -> EcceResult:
    """Cumulative-deviation calibration error (Arrieta-Ibarra et al., 2022).

    Sort by prediction and walk the cumulative sum of weighted residuals;
    under calibration the walk hovers near zero, and drift localizes
    miscalibration without any smoothing parameter.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    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) -> float

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

Source code in src/probcal/metrics/smooth.py 
 99
100
101
def emax(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None) -> float:
    """Maximum of the |LOESS(y|p) - p| distances."""
    return float(np.max(_ici_distances(y, p, frac)))

ici 

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

Integrated calibration index: weighted mean |LOESS(y|p) - p| (Austin & Steyerberg, 2019). The LOESS stage itself is unweighted (DECISIONS entry).

Source code in src/probcal/metrics/smooth.py 
75
76
77
78
79
80
81
def ici(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None) -> float:
    """Integrated calibration index: weighted mean |LOESS(y|p) - p|
    (Austin & Steyerberg, 2019). The LOESS stage itself is unweighted
    (DECISIONS entry)."""
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    c = loess(p_arr, y_arr, frac=frac)
    return float(np.average(np.abs(c - p_arr), weights=w))

smooth_ece 

smooth_ece(y: object, p: object, *, sample_weight: object = None) -> 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 — DECISIONS entry), and the reported value is the fixed point smECE(sigma) = sigma found by bisection.

Source code in src/probcal/metrics/smooth.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
def smooth_ece(y: object, p: object, *, sample_weight: object = None) -> 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 — DECISIONS entry), and the
    reported value is the fixed point ``smECE(sigma) = sigma`` found by
    bisection.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    t = logit(p_arr)
    r = y_arr - p_arr
    wn = w / w.sum()
    lo, hi = 1e-4, 2.0
    f_lo = _smece_at_sigma(t, r, wn, lo) - lo
    if f_lo <= 0.0:  # essentially perfectly calibrated at the finest scale
        return _smece_at_sigma(t, r, wn, lo)
    for _ in range(40):
        mid = 0.5 * (lo + hi)
        f_mid = _smece_at_sigma(t, r, wn, mid) - mid
        if abs(hi - lo) < 1e-4:
            break
        if f_mid > 0.0:
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)

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.

Source code in src/probcal/metrics/smooth.py 
112
113
114
115
116
117
118
119
120
121
122
123
124
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.
    """
    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) -> MetricReport

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

RETURNS DESCRIPTION
MetricReport

Point estimates and CI bounds for the full catalog. Note the caveat from the metrics chapter: a bootstrap CI around a biased estimator (plain ECE) quantifies its variance, not its bias.
Source code in src/probcal/metrics/__init__.py 
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
def evaluate(
    y: object,
    p: object,
    *,
    sample_weight: object = None,
    n_boot: int = 1000,
    seed: int = 42,
) -> MetricReport:
    """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.

    Returns
    -------
    MetricReport
        Point estimates and CI bounds for the full catalog. Note the caveat
        from the metrics chapter: a bootstrap CI around a *biased* estimator
        (plain ECE) quantifies its variance, not its bias.
    """
    from .scores import _prep

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

    rng = np.random.default_rng(seed)
    n = len(y_arr)
    boot = np.empty((n_boot, len(names)))
    for b in range(n_boot):
        idx = rng.integers(0, n, n)
        yb, pb, wb = y_arr[idx], p_arr[idx], w_arr[idx]
        if yb.min() == yb.max():  # degenerate resample: keep the point estimate
            boot[b] = values
            continue
        pm = _point_metrics(yb, pb, wb)
        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) -> 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.

Source code in src/probcal/metrics/__init__.py 
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
def reliability_summary(
    y: object, p: object, *, sample_weight: object = None
) -> 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.
    """
    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),
        e90=e90(y_arr, p_arr),
        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).

LogLossDecomposition dataclass 

LogLossDecomposition(calibration: float, refinement: float)

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

log_loss 

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

Weighted mean negative log-likelihood (strictly proper; the default selection criterion).

Source code in src/probcal/metrics/scores.py 
23
24
25
26
27
28
def log_loss(y: object, p: object, *, sample_weight: object = None) -> float:
    """Weighted mean negative log-likelihood (strictly proper; the default
    selection criterion)."""
    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).

Source code in src/probcal/metrics/scores.py 
31
32
33
34
def brier_score(y: object, p: object, *, sample_weight: object = None) -> float:
    """Weighted mean squared error of the probability forecast (strictly proper)."""
    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.

Source code in src/probcal/metrics/scores.py 
37
38
39
40
41
42
43
44
45
46
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.
    """
    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.

Source code in src/probcal/metrics/scores.py 
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
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.
    """
    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)

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.

Source code in src/probcal/metrics/scores.py 
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
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``.
    """
    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).

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

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
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).
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    shares, gaps, _, _ = _bin_gaps(y_arr, p_arr, w, n_bins, strategy)
    if norm == "l1":
        return float(np.sum(shares * gaps))
    if norm == "l2":
        return float(np.sqrt(np.sum(shares * gaps**2)))
    if norm == "max":
        return float(gaps.max())
    raise ValueError(f"norm must be 'l1', 'l2', or 'max', got {norm!r}")

ece_debiased 

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

Bias-corrected ECE: per-bin squared gaps minus the within-bin variance of the event rate, floored at zero (correction in the spirit of Bröcker 2009 / Ferro & Fricker 2012; exact estimator in the DECISIONS log).

Source code in src/probcal/metrics/binned.py 
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
def ece_debiased(
    y: object,
    p: object,
    *,
    n_bins: int = 15,
    strategy: str = "mass",
    sample_weight: object = None,
) -> float:
    """Bias-corrected ECE: per-bin squared gaps minus the within-bin variance
    of the event rate, floored at zero (correction in the spirit of Bröcker
    2009 / Ferro & Fricker 2012; exact estimator in the DECISIONS log)."""
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    shares, gaps, rates, counts = _bin_gaps(y_arr, p_arr, w, n_bins, strategy)
    corrected = np.empty_like(gaps)
    for i in range(len(gaps)):
        if counts[i] > 1:
            var_b = rates[i] * (1.0 - rates[i]) / (counts[i] - 1)
            corrected[i] = np.sqrt(max(gaps[i] ** 2 - var_b, 0.0))
        else:
            corrected[i] = gaps[i]
    return float(np.sum(shares * corrected))

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); DECISIONS entry).

Source code in src/probcal/metrics/binned.py 
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
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); DECISIONS entry).
    """
    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
    if best_b == 1:
        pb = float(np.average(p_arr, weights=w))
        yb = float(np.average(y_arr, weights=w))
        return abs(pb - yb)
    return ece(y_arr, p_arr, n_bins=best_b, strategy="mass", norm=norm, sample_weight=w)

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

Source code in src/probcal/metrics/binned.py 
122
123
124
125
126
127
128
129
130
131
132
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)."""
    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.

Source code in src/probcal/metrics/binned.py 
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
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.
    """
    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.

SpiegelhalterResult dataclass 

SpiegelhalterResult(z: float, p_value: float)

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

smooth_ece 

smooth_ece(y: object, p: object, *, sample_weight: object = None) -> 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 — DECISIONS entry), and the reported value is the fixed point smECE(sigma) = sigma found by bisection.

Source code in src/probcal/metrics/smooth.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
def smooth_ece(y: object, p: object, *, sample_weight: object = None) -> 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 — DECISIONS entry), and the
    reported value is the fixed point ``smECE(sigma) = sigma`` found by
    bisection.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    t = logit(p_arr)
    r = y_arr - p_arr
    wn = w / w.sum()
    lo, hi = 1e-4, 2.0
    f_lo = _smece_at_sigma(t, r, wn, lo) - lo
    if f_lo <= 0.0:  # essentially perfectly calibrated at the finest scale
        return _smece_at_sigma(t, r, wn, lo)
    for _ in range(40):
        mid = 0.5 * (lo + hi)
        f_mid = _smece_at_sigma(t, r, wn, mid) - mid
        if abs(hi - lo) < 1e-4:
            break
        if f_mid > 0.0:
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)

ecce 

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

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

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

Source code in src/probcal/metrics/smooth.py 
62
63
64
65
66
67
68
69
70
71
72
def ecce(y: object, p: object, *, sample_weight: object = None) -> EcceResult:
    """Cumulative-deviation calibration error (Arrieta-Ibarra et al., 2022).

    Sort by prediction and walk the cumulative sum of weighted residuals;
    under calibration the walk hovers near zero, and drift localizes
    miscalibration without any smoothing parameter.
    """
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    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) -> float

Integrated calibration index: weighted mean |LOESS(y|p) - p| (Austin & Steyerberg, 2019). The LOESS stage itself is unweighted (DECISIONS entry).

Source code in src/probcal/metrics/smooth.py 
75
76
77
78
79
80
81
def ici(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None) -> float:
    """Integrated calibration index: weighted mean |LOESS(y|p) - p|
    (Austin & Steyerberg, 2019). The LOESS stage itself is unweighted
    (DECISIONS entry)."""
    y_arr, p_arr, w = _prep(y, p, sample_weight)
    c = loess(p_arr, y_arr, frac=frac)
    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) -> float

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

Source code in src/probcal/metrics/smooth.py 
89
90
91
def e50(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None) -> float:
    """Median of the |LOESS(y|p) - p| distances."""
    return float(np.quantile(_ici_distances(y, p, frac), 0.5))

e90 

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

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

Source code in src/probcal/metrics/smooth.py 
94
95
96
def e90(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None) -> float:
    """90th percentile of the |LOESS(y|p) - p| distances."""
    return float(np.quantile(_ici_distances(y, p, frac), 0.9))

emax 

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

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

Source code in src/probcal/metrics/smooth.py 
 99
100
101
def emax(y: object, p: object, *, frac: float = 0.75, sample_weight: object = None) -> float:
    """Maximum of the |LOESS(y|p) - p| distances."""
    return float(np.max(_ici_distances(y, p, frac)))

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.

Source code in src/probcal/metrics/smooth.py 
112
113
114
115
116
117
118
119
120
121
122
123
124
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.
    """
    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.

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.

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

Source code in src/probcal/metrics/regression.py 
16
17
18
19
20
21
22
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))."""
    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.

Source code in src/probcal/metrics/regression.py 
25
26
27
28
29
30
31
32
def calibration_slope(y: object, p: object, *, sample_weight: object = None) -> float:
    """Cox calibration slope: < 1 means overfitting/overconfident spread,
    > 1 underfitting."""
    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).

Source code in src/probcal/metrics/regression.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
def calibration_test(
    y: object, p: object, *, sample_weight: object = None
) -> CalibrationTestResult:
    """Likelihood-ratio test of joint calibration (alpha, beta) = (0, 1)."""
    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).

Source code in src/probcal/metrics/regression.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
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)."""
    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.

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

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.

Source code in src/probcal/metrics/grade.py 
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 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
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.
    """
    y_arr, p_arr, _ = _prep(y, p, None)
    _check_weights(sample_weight, len(y_arr))
    g_arr = np.asarray(grades)
    labels, n, k, pd = _per_grade(y_arr, p_arr, g_arr)
    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.

Source code in src/probcal/metrics/grade.py 
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
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.
    """
    y_arr, p_arr, _ = _prep(y, p, None)
    _check_weights(sample_weight, len(y_arr))
    g_arr = np.asarray(grades)
    labels, n, k, pd = _per_grade(y_arr, p_arr, g_arr)
    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).

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

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

Source code in src/probcal/metrics/kernel.py 
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
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.
    """
    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,
    )