import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingRegressor
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
h=pd.read_csv("cali_housing.csv"); hf=[c for c in h.columns if c!="MedHouseVal"]
X=h[hf].values; y=h["MedHouseVal"].values
Xtr,Xtmp,ytr,ytmp=train_test_split(X,y,test_size=0.4,random_state=0); Xcal,Xte,ycal,yte=train_test_split(Xtmp,ytmp,test_size=0.5,random_state=0)
m=GradientBoostingRegressor(n_estimators=300,max_depth=3,learning_rate=0.05).fit(Xtr,ytr)
scores=np.abs(ycal-m.predict(Xcal))                              # conformity = |residual|
def qhat(s,al):
    """The conformal quantile. method='higher' matters: the guarantee is stated for the EMPIRICAL
    quantile, and numpy's default linear interpolation returns something slightly smaller."""
    return np.quantile(s, np.ceil((len(s)+1)*(1-al))/len(s), method="higher")
alphas=np.array([0.02,0.05,0.1,0.15,0.2,0.3]); cov=[]; wid=[]
for al in alphas:
    q=qhat(scores,al); pr=m.predict(Xte); cov.append(np.mean((yte>=pr-q)&(yte<=pr+q))); wid.append(2*q)
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].plot(1-alphas,cov,"o-",color=BLUE,lw=2,label="empirical coverage"); ax[0].plot([0.6,1],[0.6,1],"k--",lw=1,label="target = 1−α")
ax[0].set_xlabel("target coverage 1−α"); ax[0].set_ylabel("empirical coverage"); ax[0].set_title("The guarantee holds — coverage tracks the target"); ax[0].legend(fontsize=8)
pr=m.predict(Xte); q=qhat(scores,0.1); o=np.argsort(pr)[::30]
ax[1].errorbar(np.arange(len(o)),pr[o],yerr=q,fmt="o",ms=3,color=BLUE,ecolor=GREY,capsize=2,label="90% conformal interval")
ax[1].scatter(np.arange(len(o)),yte[o],color=RED,s=14,zorder=5,label="actual"); ax[1].set_xlabel("test point (sorted)"); ax[1].set_ylabel("median house value ($100k)"); ax[1].set_title("90% intervals (constant width)"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"target 0.90 -> empirical coverage {cov[2]:.3f}, interval width {wid[2]:.2f}. Coverage matched at every level, from a")
print("plain gradient booster with no distributional assumptions -- the finite-sample conformal guarantee.")

target 0.90 -> empirical coverage 0.897, interval width 1.57. Coverage matched at every level, from a
plain gradient booster with no distributional assumptions -- the finite-sample conformal guarantee.

lo=GradientBoostingRegressor(loss="quantile",alpha=0.05,n_estimators=300,max_depth=3,learning_rate=0.05).fit(Xtr,ytr)
hi=GradientBoostingRegressor(loss="quantile",alpha=0.95,n_estimators=300,max_depth=3,learning_rate=0.05).fit(Xtr,ytr)
E=np.maximum(lo.predict(Xcal)-ycal, ycal-hi.predict(Xcal))       # CQR conformity score
qc=qhat(E,0.1); loT,hiT=lo.predict(Xte)-qc, hi.predict(Xte)+qc
covC=np.mean((yte>=loT)&(yte<=hiT)); widC=hiT-loT
covS=np.abs(yte-pr)<=q; covCQ=(yte>=loT)&(yte<=hiT)
dec=np.digitize(pr,np.quantile(pr,np.linspace(0,1,11)[1:-1]))     # deciles of the model's own prediction
cs=[covS[dec==k].mean() for k in range(10)]; cc=[covCQ[dec==k].mean() for k in range(10)]
wc=[widC[dec==k].mean() for k in range(10)]

fig,ax=plt.subplots(1,3,figsize=(16,4.2))
ax[0].hist(widC,bins=30,color=GREEN,edgecolor="white"); ax[0].axvline(2*q,color=GREY,lw=2,label="split conformal (constant)")
ax[0].set_xlabel("interval width"); ax[0].set_ylabel("# test points"); ax[0].set_title("CQR widths vary"); ax[0].legend(fontsize=8)
xd=np.arange(10)
ax[1].plot(xd,cs,"o-",color=GREY,lw=2,label=f"split conformal (marginal {covS.mean():.3f})")
ax[1].plot(xd,cc,"o-",color=GREEN,lw=2,label=f"CQR (marginal {covCQ.mean():.3f})")
ax[1].axhline(0.9,color=RED,ls="--",lw=1.5,label="target 0.90")
ax[1].set_xlabel("decile of predicted value"); ax[1].set_ylabel("coverage within decile")
ax[1].set_title("Conditional coverage — the marginal guarantee hides this"); ax[1].legend(fontsize=7); ax[1].set_ylim(0.6,1.02)
ax[2].plot(xd,[2*q]*10,"o-",color=GREY,lw=2,label="split conformal"); ax[2].plot(xd,wc,"o-",color=GREEN,lw=2,label="CQR")
ax[2].set_xlabel("decile of predicted value"); ax[2].set_ylabel("mean interval width"); ax[2].set_title("Where the width goes"); ax[2].legend(fontsize=8)
plt.tight_layout(); plt.show()

print(f"Both methods hit the marginal target: split conformal {covS.mean():.3f}, CQR {covCQ.mean():.3f}, against 0.90.")
print(f"On average width, split conformal is the NARROWER of the two: {2*q:.3f} against CQR's {widC.mean():.3f}.")
print("On that evidence alone CQR looks like a step backwards. Coverage by decile says otherwise:\n")
print(f"  {'decile':>7} {'mean prediction':>16} {'split coverage':>16} {'CQR coverage':>14} {'CQR width':>11}")
for k in range(10):
    print(f"  {k:>7} {pr[dec==k].mean():>16.2f} {cs[k]:>16.3f} {cc[k]:>14.3f} {wc[k]:>11.2f}")
print(f"\nSplit conformal's constant interval covers {max(cs):.0%} of the cheapest decile and {min(cs):.0%} of the eighth -- a spread of")
print(f"{max(cs)-min(cs):.3f}. It is a 90% interval only on average; for expensive houses it is closer to a {min(cs):.0%} interval, and")
print("nobody pricing an individual property would be served by that. The marginal guarantee is real and it is exactly")
print("as narrow as it sounds.")
print(f"\nCQR spends its width where the uncertainty is -- {min(wc):.2f} in the cheapest decile rising to {max(wc):.2f} in the dearest -- and")
print(f"cuts the coverage spread from {max(cs)-min(cs):.3f} to {max(cc)-min(cc):.3f}. It is {widC.mean()/(2*q)-1:+.0%} wider on average and far more honest per")
print("stratum. That trade, not the varying widths themselves, is the reason to prefer it.")
print(f"\nNeither achieves exact conditional coverage, and no distribution-free method can: that is a theorem")
print("(Foygel Barber et al., 2021), not a limitation of this implementation. Approximate conditional validity is the")
print("most that is available, and it has to be checked rather than assumed.")

Both methods hit the marginal target: split conformal 0.897, CQR 0.895, against 0.90.
On average width, split conformal is the NARROWER of the two: 1.574 against CQR's 1.710.
On that evidence alone CQR looks like a step backwards. Coverage by decile says otherwise:

   decile  mean prediction   split coverage   CQR coverage   CQR width
        0             0.72            0.981          0.937        1.06
        1             1.06            0.973          0.860        1.27
        2             1.31            0.969          0.864        1.30
        3             1.53            0.956          0.903        1.35
        4             1.75            0.944          0.901        1.43
        5             1.96            0.930          0.920        1.60
        6             2.23            0.862          0.881        1.87
        7             2.62            0.821          0.898        2.17
        8             3.18            0.719          0.855        2.50
        9             4.20            0.811          0.935        2.56

Split conformal's constant interval covers 98% of the cheapest decile and 72% of the eighth -- a spread of
0.262. It is a 90% interval only on average; for expensive houses it is closer to a 72% interval, and
nobody pricing an individual property would be served by that. The marginal guarantee is real and it is exactly
as narrow as it sounds.

CQR spends its width where the uncertainty is -- 1.06 in the cheapest decile rising to 2.56 in the dearest -- and
cuts the coverage spread from 0.262 to 0.082. It is +9% wider on average and far more honest per
stratum. That trade, not the varying widths themselves, is the reason to prefer it.

Neither achieves exact conditional coverage, and no distribution-free method can: that is a theorem
(Foygel Barber et al., 2021), not a limitation of this implementation. Approximate conditional validity is the
most that is available, and it has to be checked rather than assumed.

from sklearn.datasets import load_digits
from sklearn.ensemble import RandomForestClassifier
dig=load_digits(); rng=np.random.default_rng(0); Xd=dig.data+rng.normal(0,6,dig.data.shape); yd=dig.target
Xtr,Xtmp,ytr,ytmp=train_test_split(Xd,yd,train_size=250,random_state=0,stratify=yd); Xcal,Xte,ycal,yte=train_test_split(Xtmp,ytmp,test_size=0.5,random_state=0)
clf=RandomForestClassifier(200,random_state=0).fit(Xtr,ytr)
def aps_score(P,yv):
    o=np.argsort(-P,1); s=np.zeros(len(yv))
    for i in range(len(yv)):
        c=0.0
        for k in o[i]:
            c+=P[i,k]
            if k==yv[i]: s[i]=c; break
    return s
sc=aps_score(clf.predict_proba(Xcal),ycal); qh=qhat(sc,0.1)
Pte=clf.predict_proba(Xte); o=np.argsort(-Pte,1); sets=np.zeros_like(Pte,bool); sets_loose=np.zeros_like(Pte,bool)
for i in range(len(Pte)):
    c=0.0
    for j,k in enumerate(o[i]):
        c+=Pte[i,k]                                              # cumulative INCLUDING this label
        if c<=qh or j==0: sets[i,k]=True                         # keep it only if the score still fits; never empty
        else: break
    c=0.0
    for k in o[i]:                                               # the tempting variant: add, then test
        sets_loose[i,k]=True; c+=Pte[i,k]
        if c>=qh: break
cov=sets[np.arange(len(yte)),yte].mean(); size=sets.sum(1); correct=clf.predict(Xte)==yte
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
ax[0].hist(size,bins=range(1,int(size.max())+2),color=PURP,edgecolor="white",align="left"); ax[0].set_xlabel("prediction-set size"); ax[0].set_ylabel("# test digits"); ax[0].set_title(f"Set sizes (coverage {cov:.3f}, target 0.90)")
ax[1].bar(["model correct\n(easy)","model wrong\n(hard)"],[size[correct].mean(),size[~correct].mean()],color=[GREEN,RED]); ax[1].set_ylabel("avg set size"); ax[1].set_title("Set size adapts to difficulty")
for i,v in enumerate([size[correct].mean(),size[~correct].mean()]): ax[1].text(i,v+0.05,f"{v:.2f}",ha="center")
plt.tight_layout(); plt.show()
cov_l=sets_loose[np.arange(len(yte)),yte].mean(); size_l=sets_loose.sum(1)
print(f"APS conformal on noisy digits (base accuracy {correct.mean():.2f}, calibration n={len(ycal)}), target coverage 0.90:")
print(f"   set = {{k : cumulative_<=k <= qhat}}          coverage {cov:.4f}   mean set size {size.mean():.2f}")
print(f"   the 'add, then test' variant                 coverage {cov_l:.4f}   mean set size {size_l.mean():.2f}")
print(f"\nThe second construction carries one extra label in every set. It costs {size_l.mean()/size.mean()-1:+.0%} in set size and pushes coverage")
print(f"to {cov_l:.3f} against a {0.90:.2f} target -- over-covering is not free, because an interval or a set is only useful to the")
print("extent that it is small. Matching the set construction to the score definition is the whole of the fix.")
print(f"\nEven done correctly the coverage is {cov:.3f} rather than 0.900. That is expected: without the randomised")
print("tie-breaking term this form of APS over-covers by construction, and the guarantee is one-sided anyway.")
print(f"\nSet size carries the uncertainty: {size[correct].mean():.2f} labels when the model gets it right against {size[~correct].mean():.2f} when it does not.")
print("A conformal classifier does not become more accurate -- it becomes honest about when to hedge.")

APS conformal on noisy digits (base accuracy 0.68, calibration n=773), target coverage 0.90:
   set = {k : cumulative_<=k <= qhat}          coverage 0.9160   mean set size 2.83
   the 'add, then test' variant                 coverage 0.9548   mean set size 3.79

The second construction carries one extra label in every set. It costs +34% in set size and pushes coverage
to 0.955 against a 0.90 target -- over-covering is not free, because an interval or a set is only useful to the
extent that it is small. Matching the set construction to the score definition is the whole of the fix.

Even done correctly the coverage is 0.916 rather than 0.900. That is expected: without the randomised
tie-breaking term this form of APS over-covers by construction, and the guarantee is one-sided anyway.

Set size carries the uncertainty: 2.65 labels when the model gets it right against 3.22 when it does not.
A conformal classifier does not become more accurate -- it becomes honest about when to hedge.

Model Evaluation IV — Conformal Prediction (finale of the ML arc)¶

Distribution-free prediction sets with a guaranteed coverage rate (Vovk; Angelopoulos & Bates)¶

1. Split conformal for regression — the guarantee¶

2. Adaptive intervals — conformalized quantile regression¶

3. Classification — prediction sets with adaptive size¶

4. Conformal vs Bayesian — and the end of the ML arc¶