Merge #969: Reject undefined Hosmer-Lemeshow calibration results

specialized/specialized-model-qa.md

@@ -167,6 +167,8 @@ def discrimination_report(y_true: pd.Series, y_score: pd.Series) -> dict:
 ### Calibration Test (Hosmer-Lemeshow)
 
 ```python
+import numpy as np
+import pandas as pd
 from scipy.stats import chi2
 
 def hosmer_lemeshow_test(
@@ -175,8 +177,18 @@ def hosmer_lemeshow_test(
     """
     Hosmer-Lemeshow goodness-of-fit test for calibration.
     p-value < 0.05 suggests significant miscalibration.
+    A non-significant result is not proof of calibration. Undefined tests
+    raise ValueError rather than returning NaN and a misleading verdict.
     """
+    if not isinstance(groups, int) or isinstance(groups, bool) or groups < 3:
+        raise ValueError("HL requires an integer group count >= 3")
+    if y_true.empty or not y_true.index.equals(y_pred.index):
+        raise ValueError("HL requires nonempty, aligned observations")
     data = pd.DataFrame({"y": y_true, "p": y_pred})
+    if not np.isfinite(data.to_numpy(dtype=float)).all():
+        raise ValueError("HL observations and probabilities must be finite")
+    if not data["y"].isin([0, 1]).all() or not data["p"].between(0, 1).all():
+        raise ValueError("HL requires binary outcomes and probabilities in [0, 1]")
     data["bucket"] = pd.qcut(data["p"], groups, duplicates="drop")
 
     agg = data.groupby("bucket", observed=True).agg(
@@ -185,13 +197,15 @@ def hosmer_lemeshow_test(
         expected=("p", "sum"),
     )
 
-    hl_stat = (
-        ((agg["observed"] - agg["expected"]) ** 2)
-        / (agg["expected"] * (1 - agg["expected"] / agg["n"]))
-    ).sum()
+    # Ties can collapse qcut groups; zero/one probability groups have zero
+    # expected variance. Neither case permits an ordinary chi-square result.
+    variance = agg["expected"] * (1 - agg["expected"] / agg["n"])
+    if len(agg) < 3 or not (variance > 0).all():
+        raise ValueError("HL needs at least three groups with positive expected variance")
+    hl_stat = (((agg["observed"] - agg["expected"]) ** 2) / variance).sum()
 
     dof = len(agg) - 2
-    p_value = 1 - chi2.cdf(hl_stat, dof)
+    p_value = chi2.sf(hl_stat, dof)
 
     return {
         "HL_statistic": round(hl_stat, 4),