Conformal prediction — নিশ্চয়তাসহ অনুমানের পরিসর

পর্ব ২৪-এ quantile regression দিয়ে ১০–৯০ পরিসর বানিয়ে held-out ডেটায় coverage মেপেছেন, আর দেখেছেন যে ৮০% আসে শুধু model ঠিক থাকলে। পর্ব ২১-এ bootstrap দিয়ে অনিশ্চয়তা মেপেছেন, পর্ব ৪৬-এ Bayesian পরিসর দেখেছেন, আর পর্ব ২৬-এ সম্ভাবনা calibrate করেছেন। সবগুলোর একটি দুর্বলতা আছে: model-এর অনুমান ভুল হলে পরিসরও ভুল। এই পর্বের পদ্ধতি সেই নির্ভরতা কাটায়। ডেটা ভাগ করার নিয়ম (AI-এর ভিত্তি কোর্সের পর্ব ৫৬) আর মূল্যায়নের মাপ (একই কোর্সের পর্ব ৬০) ধরে নিচ্ছি। চারটি চালানো যায় এমন ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, pandas, scikit-learn আর matplotlib); প্রথমবার scikit-learn নামতে প্রায় দশ সেকেন্ড লাগে।

১ · একটি সংখ্যা নয়, একটি প্রতিশ্রুতি — conformal-এর মূল ধারণা

Conformal prediction ঠিক এই কাজ করে। যেকোনো model নিন: linear, boosting, neural network. আলাদা একটি calibration set-এ মাপুন, model-এর ভুল কত বড় হয়। সেই ভুলগুলোর একটি নির্দিষ্ট ক্রমের মান দিয়ে নতুন উদাহরণের চারপাশে পরিসর বানান। ফলাফল একটি prediction interval (regression-এ) বা একটি prediction set (classification-এ), সঙ্গে একটি গাণিতিক প্রতিশ্রুতি:

P( Y_নতুন ∈ C(X_নতুন) ) ≥ ১ − α

এখানে C(x) হলো x-এর জন্য বানানো পরিসর বা সেট, আর α হলো আপনি কতটা ভুল মেনে নিতে রাজি (α = ০.১ মানে লক্ষ্য ৯০%). সম্ভাবনাটি নেওয়া হয় calibration ডেটা আর নতুন উদাহরণ, দুটির এলোমেলোতা মিলিয়ে। একমাত্র বড় শর্ত exchangeability: calibration-এর উদাহরণগুলো আর নতুন উদাহরণটির ক্রম বদলালে যৌথ বণ্টন বদলায় না। i.i.d. ডেটা exchangeable, কিন্তু শর্তটি তার চেয়েও একটু দুর্বল। Model ঠিক, ভুল normal, বা নমুনা বড়, এর কোনোটিই লাগে না। Model খারাপ হলে প্রতিশ্রুতি ভাঙে না; শুধু পরিসর চওড়া হয়, অর্থাৎ কম কাজের হয়।

ভুলের মাপকে বলে nonconformity score: একটি উদাহরণ model-এর চোখে কতটা “বেমানান”। Regression-এ সবচেয়ে সরল score ∣y − ŷ∣; classification-এ ১ − p̂(আসল শ্রেণি)। Score যত বড়, উদাহরণ তত অস্বাভাবিক। পুরো পদ্ধতির ধার কোন score বাছছেন তার ওপর: ভালো score মানে সরু পরিসর, খারাপ score মানে চওড়া পরিসর; কিন্তু দুই ক্ষেত্রেই coverage-এর প্রতিশ্রুতি একই।

২ · Split conformal হাতে — score, ⌈(n + ১)(১ − α)⌉ আর interval

Split conformal (inductive conformal নামেও পরিচিত) সবচেয়ে প্রচলিত রূপ। ধাপগুলো ছবির মতোই: train ভাগে model fit, calibration ভাগে n টি score s_১, …, s_n, তারপর

k = ⌈(n + ১)(১ − α)⌉,   q̂ = সাজানো score-এর k-তম মান,   C(x) = [ŷ(x) − q̂, ŷ(x) + q̂]

⌈ ⌉ মানে ওপরের দিকে পূর্ণসংখ্যায় তোলা (ceiling). কেন n নয়, n + ১? নতুন উদাহরণের score-টিকে n টি calibration score-এর সঙ্গে মিলিয়ে মোট n + ১টি সংখ্যা ভাবুন। Exchangeable হলে নতুন score-টি এই n + ১টির মধ্যে যেকোনো স্থানে বসার সম্ভাবনা সমান, প্রতিটি ১/(n + ১)। তাই সে প্রথম k টি স্থানের মধ্যে থাকার সম্ভাবনা অন্তত k/(n + ১) ≥ ১ − α. নিচের stepper নয়টি score-এ পুরো হিসাবটি করে।

এখন একই নিয়ম আসল আকারে। নিচের module বাজারঘরের (কাল্পনিক দোকান) তিন রকম ডেটা seed দিয়ে বানায়, কোনো আসল দোকানের নয়। delivery-তে ডেলিভারির সময় দূরত্ব, বৃষ্টি আর চর-হাওরের মতো দূরবর্তী এলাকার (remote, প্রায় ২০% অর্ডার) ওপর নির্ভর করে, আর ছড়ানো দূরত্বের সঙ্গে বাড়ে; দূরবর্তী এলাকায় ছড়ানো আরও আড়াই গুণ। অর্থাৎ ডেটা heteroscedastic. returns ফেরত দেওয়া পণ্যের চারটি কারণ বানায়, আর monthly আট মাসের অর্ডার, যার ষষ্ঠ মাস থেকে বর্ষা।

# বাজারঘর (কাল্পনিক দোকান): seed দিয়ে বানানো ডেলিভারি আর ফেরতের ডেটা
import numpy as np
import pandas as pd

def delivery(n=3000, seed=58):
    rng = np.random.default_rng(seed)
    km = rng.uniform(2, 80, n).round(1)
    rain = rng.binomial(1, 0.3, n)
    remote = rng.binomial(1, 0.2, n)        # চর/হাওর এলাকা: সময়ের ছড়ানো বেশি
    sd = (0.5 + 0.06 * km) * (1 + 1.5 * remote)
    hours = 6 + 0.25 * km + 3 * rain + 2 * remote + sd * rng.standard_normal(n)
    X = pd.DataFrame({"km": km, "rain": rain, "remote": remote})
    return X, np.maximum(hours, 0.5).round(2)

REASONS = ["size", "damaged", "late", "wrong_item"]   # ফেরতের চারটি কারণ

def returns(n=4000, seed=58):
    rng = np.random.default_rng(seed)
    cat = rng.choice(3, n, p=[0.5, 0.3, 0.2])          # fashion, electronics, home
    hours = rng.gamma(3, 8, n).round(1)                 # ডেলিভারিতে কত ঘণ্টা
    items = rng.integers(1, 6, n)                       # অর্ডারে কয়টি পণ্য
    logit = np.column_stack([
        1.0 + 2.5 * (cat == 0) - 1.5 * (cat == 1),
        2.8 * (cat == 1) + 1.2 * (cat == 2),
        -2.2 + 0.09 * hours,
        -2.6 + 0.8 * items])
    p = np.exp(logit) / np.exp(logit).sum(1, keepdims=True)
    y = (p.cumsum(1) > rng.random((n, 1))).argmax(1)
    X = np.column_stack([cat == 0, cat == 1, hours / 24, items]).astype(float)
    return X, y

def monthly(per_month=300, seed=58, monsoon_from=6):
    rng = np.random.default_rng(seed)
    month = np.repeat(np.arange(1, 9), per_month)
    km = rng.uniform(2, 80, month.size).round(1)
    rain = rng.binomial(1, 0.3, month.size)
    wet = month >= monsoon_from                         # বর্ষা: দূরের রাস্তায় দেরি
    sd = (0.5 + 0.06 * km) * np.where(wet, 1.6, 1.0)
    hours = 6 + 0.25 * km + 3 * rain + 3 * wet * (km > 40) + sd * rng.standard_normal(month.size)
    return pd.DataFrame({"month": month, "km": km, "rain": rain,
                         "hours": np.maximum(hours, 0.5).round(2)})

প্রথম ব্লক ২,০০০টি অর্ডার থেকে ১,০০০টি train, ৫০০টি calibration আর ৫০০টি test নেয়, LinearRegression শেখায়, আর উপরের নিয়মে পরিসর বানায়। তারপর একই কাজ ২০০টি আলাদা এলোমেলো ভাগে করে প্রতিবারের test coverage গোনে। ছবির লেখা ইংরেজিতে, কারণ পাতার matplotlib বাংলা অক্ষর ঠিকমতো আঁকে না।

import numpy as np
import pandas as pd                  # bazar_conformal module-এর জন্য লাগে
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from bazar_conformal import delivery

X, y = delivery(n=2000)
X, alpha = X.to_numpy(), 0.1

def split_conformal(seed, n_cal=500):
    idx = np.random.default_rng(seed).permutation(len(y))
    tr, cal, te = idx[:1000], idx[1000:1000 + n_cal], idx[1000 + n_cal:]
    model = LinearRegression().fit(X[tr], y[tr])
    scores = np.abs(y[cal] - model.predict(X[cal]))    # residual score
    k = int(np.ceil((n_cal + 1) * (1 - alpha)))         # কততম ছোট score
    q = np.sort(scores)[k - 1]
    pred = model.predict(X[te])                         # interval: pred ± q
    return k, q, np.mean(np.abs(y[te] - pred) <= q)

k, q, cov = split_conformal(seed=0)
print(f"n_cal=500  k={k}  q_hat={q:.2f} hours  test coverage={cov:.3f}")
covs = np.array([split_conformal(s)[2] for s in range(200)])
print(f"200 splits: mean={covs.mean():.4f}  sd={covs.std():.4f}  "
      f"min={covs.min():.3f}  max={covs.max():.3f}")
print(f"theory: mean in [{1 - alpha:.4f}, {k / 501:.4f}]")
print("splits below 0.88:", int((covs < 0.88).sum()), "of 200")

fig, ax = plt.subplots(figsize=(6, 3.4))
fig.subplots_adjust(bottom=0.17)
ax.hist(covs, bins=20, color="tab:blue", alpha=0.8)
ax.axvline(0.9, color="r", ls="--", label="target 1 - alpha = 0.9")
ax.set_xlabel("test coverage of one split"); ax.set_ylabel("number of splits")
ax.legend(fontsize=8); plt.show()

n_cal=500  k=451  q_hat=7.09 hours  test coverage=0.904
200 splits: mean=0.8991  sd=0.0181  min=0.840  max=0.948
theory: mean in [0.9000, 0.9002]
splits below 0.88: 29 of 200

৩ · “৯০% coverage” ঠিকভাবে পড়া — marginal, প্রতিটি ক্রেতার নয়

Conformal-এর প্রতিশ্রুতি marginal coverage: সব x আর সব calibration ভাগ মিলিয়ে গড়ে অন্তত ১ − α. এটি conditional coverage নয়, অর্থাৎ “এই নির্দিষ্ট দূরত্বের, এই এলাকার অর্ডারের জন্য ৯০%” নয়। দুটি প্রশ্ন আলাদা, আর ভুল পড়লে সবচেয়ে বেশি ঠকেন সংখ্যালঘু দলের ক্রেতারা। নিচের stepper একটি সরল মিশ্রণ দিয়ে কথাটি পরিষ্কার করে।

কোনো model-নির্ভরতা ছাড়া প্রতিটি x-এর জন্য conditional coverage দেওয়া সম্ভব কি না, এর উত্তর গবেষণায় মোটামুটি স্থির: সাধারণভাবে না। Vovk (২০১২), Lei আর Wasserman (২০১৪), আর পরে Foygel Barber ও সহকর্মীরা (২০২১) দেখান, বণ্টন-মুক্ত পদ্ধতিতে সঠিক conditional coverage চাইলে সাধারণ পরিস্থিতিতে পরিসর অসীম চওড়া হয়ে যায়। বাস্তবের পথ দুটি: score এমনভাবে বাছা যাতে প্রস্থ ঝুঁকির সঙ্গে বদলায় (অংশ ৪), আর জানা দলগুলোর জন্য আলাদা calibration (অংশ ৫-এর Mondrian).

৪ · Conformalized quantile regression — প্রস্থ যেখানে দরকার

∣y − ŷ∣ score-এর পরিসর সব অর্ডারে একই চওড়া। Heteroscedastic ডেটায় তাই কাছের অর্ডারে অযথা চওড়া, দূরের অর্ডারে সরু। পর্ব ২৪-এর quantile regression প্রস্থকে ঝুঁকির সঙ্গে বদলায়, কিন্তু coverage-এর নিশ্চয়তা দেয় না। Conformalized quantile regression (CQR; Romano, Patterson আর Candès, ২০১৯) দুটিকে মেলায়:

1. Train ভাগে দুটি quantile model শেখান: q̂_lo(x) (৫তম percentile) আর q̂_hi(x) (৯৫তম)। পর্ব ৩৫-এর মতো HistGradientBoostingRegressor(loss="quantile") এ কাজ করে, pinball loss কমিয়ে।
2. Calibration-এ score E_i = max(q̂_lo(x_i) − y_i, y_i − q̂_hi(x_i)). y ব্যান্ডের ভেতরে থাকলে E ঋণাত্মক (নিকটতম সীমা থেকে কত ভেতরে), বাইরে থাকলে ধনাত্মক (কত বাইরে)।
3. আগের মতো q̂ = E-এর k-তম ছোট মান, আর পরিসর [q̂_lo(x) − q̂, q̂_hi(x) + q̂]। q̂ ধনাত্মক হলে ব্যান্ড চওড়া হয়, ঋণাত্মক হলে সরু।

Quantile model ভালো হলে q̂ ছোট, আর প্রস্থ মূলত model-এর নিজের। Model খারাপ হলেও marginal coverage রক্ষা পায়, কারণ যুক্তিটি অংশ ২-এর ক্রমের যুক্তিই। নিচের ব্লক একই ডেটায় সাধারণ split conformal (HistGradientBoosting-এর গড়-model আর ∣y − ŷ∣) আর CQR পাশাপাশি চালায়, তারপর তিনটি দলে coverage আর গড় প্রস্থ মাপে: কাছে (৪০ km পর্যন্ত), দূরে, আর দূরবর্তী এলাকা (remote).

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import HistGradientBoostingRegressor as HGB
from bazar_conformal import delivery
X, y = delivery(n=3000)
tr, cal, te = np.split(np.random.default_rng(1).permutation(3000), [1500, 2250])
k = int(np.ceil((len(cal) + 1) * 0.9))                  # alpha = 0.1
kth = lambda s: np.sort(s)[k - 1]
fit = lambda **kw: HGB(random_state=0, **kw).fit(X.iloc[tr], y[tr])
mean = fit()                                            # সাধারণ split conformal
q1 = kth(np.abs(y[cal] - mean.predict(X.iloc[cal])))
m = mean.predict(X.iloc[te]); lo1, hi1 = m - q1, m + q1
qlo = fit(loss="quantile", quantile=0.05)               # CQR: দুই quantile model
qhi = fit(loss="quantile", quantile=0.95)
E = np.maximum(qlo.predict(X.iloc[cal]) - y[cal], y[cal] - qhi.predict(X.iloc[cal]))
q2 = kth(E)
lo2, hi2 = qlo.predict(X.iloc[te]) - q2, qhi.predict(X.iloc[te]) + q2
yt, Xt = y[te], X.iloc[te]
grp = np.where(Xt.remote == 1, "remote", np.where(Xt.km > 40, "far", "near"))
res = pd.DataFrame({"split_cov": (yt >= lo1) & (yt <= hi1), "split_w": hi1 - lo1,
                    "cqr_cov": (yt >= lo2) & (yt <= hi2), "cqr_w": hi2 - lo2})
tab = res.groupby(grp).mean()
tab.loc["all"] = res.mean()
tab["n"] = pd.Series(grp).value_counts().reindex(tab.index).fillna(len(te)).astype(int)
print(f"q_split={q1:.2f}  q_cqr={q2:.2f}"); print(tab.round(3))
fig, ax = plt.subplots(figsize=(6, 3.4)); fig.subplots_adjust(bottom=0.17)
for r, c in [(0, "tab:blue"), (1, "tab:orange")]:
    s = (Xt.remote == r).to_numpy()
    ax.scatter(Xt.km[s], (hi2 - lo2)[s], s=8, c=c, label=f"CQR, remote={r}")
ax.axhline(2 * q1, color="k", ls="--", label="split conformal (constant)")
ax.set_xlabel("distance from hub (km)"); ax.set_ylabel("interval width (hours)")
ax.legend(fontsize=8); plt.show()

q_split=7.59  q_cqr=0.88
        split_cov  split_w  cqr_cov   cqr_w    n
far         0.939   15.173    0.890  13.354  310
near        1.000   15.173    0.910   7.121  290
remote      0.727   15.173    0.880  22.790  150
all         0.920   15.173    0.896  12.831  750

৫ · Classification-এ prediction set — LAC, APS আর Mondrian

বাজারঘরের ফেরত-দল প্রতিটি ফেরতের কারণ আগেভাগে অনুমান করতে চায়: মাপ মেলেনি (size), ভাঙা বা নষ্ট (damaged), দেরি (late), নাকি ভুল পণ্য (wrong_item). এটি একটি multiclass classification (পর্ব ২৫)। একটি শ্রেণির বদলে conformal একটি সেট দেয়: “কারণটি {size, damaged}-এর মধ্যে, অন্তত ৯০% নিশ্চয়তায়।” সেট ছোট মানে model নিশ্চিত, বড় মানে দ্বিধায়; সেটের আকার নিজেই অনিশ্চয়তার মাপ।

সবচেয়ে সরল score LAC (least ambiguous set-valued classifier; Sadinle, Lei আর Wasserman, ২০১৯): s = ১ − p̂(আসল শ্রেণি)। Calibration-এ q̂ বের করার পর নতুন উদাহরণের সেট = যেসব শ্রেণির p̂ ≥ ১ − q̂. গড় সেটের আকার এতে যতটা সম্ভব ছোট হয়, কিন্তু সহজ উদাহরণে বেশি আর কঠিন (অনিশ্চিত) উদাহরণে কম coverage দেওয়ার প্রবণতা থাকে। APS (adaptive prediction sets; Romano, Sesia আর Candès, ২০২০) অন্যভাবে ভাবে: সম্ভাবনা বড় থেকে ছোট সাজিয়ে যোগ করতে থাকুন যতক্ষণ না আসল শ্রেণি আসে; সেই যোগফলই score. তাতে সেটের আকার উদাহরণের অনিশ্চয়তার সঙ্গে আরও মানানসই হয়, বিনিময়ে গড় সেট একটু বড়।

Marginal coverage এখানেও শ্রেণিগুলোর মধ্যে অসমান হতে পারে: বিরল শ্রেণির উদাহরণ সাধারণত কম সম্ভাবনা পায়, তাই তাদের score বড়, আর একটি সাধারণ q̂ তাদের প্রায়ই বাদ দেয়। প্রতিকার Mondrian conformal: প্রতিটি শ্রেণির জন্য আলাদা q̂, শুধু সেই শ্রেণির calibration উদাহরণ থেকে। তাতে প্রতিটি শ্রেণি নিজের ১ − α প্রতিশ্রুতি পায়, যদি প্রতিটি শ্রেণিতে যথেষ্ট calibration উদাহরণ থাকে। একই ভাবনা regression-এর দলেও খাটে (কাছে, দূরে, দূরবর্তী এলাকা)। নিচের ব্লক LogisticRegression দিয়ে দুটিই চালায়; ডেটায় কারণগুলোর ভাগ অসমান (class imbalance).

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from bazar_conformal import returns, REASONS

X, y = returns(n=4000)
tr, cal, te = np.split(np.random.default_rng(4).permutation(len(y)), [2000, 3000])
clf = LogisticRegression(max_iter=1000).fit(X[tr], y[tr])
P_cal, P_te = clf.predict_proba(X[cal]), clf.predict_proba(X[te])

def qhat(s, alpha=0.1):
    k = int(np.ceil((len(s) + 1) * (1 - alpha)))
    return np.inf if k > len(s) else np.sort(s)[k - 1]

s_cal = 1 - P_cal[np.arange(len(cal)), y[cal]]          # LAC score = 1 - p(আসল শ্রেণি)
q = qhat(s_cal)
sets = P_te >= 1 - q                                     # প্রতি সারিতে কোন শ্রেণি থাকবে
hit = sets[np.arange(len(te)), y[te]]
print(f"q_hat={q:.3f}  -> keep reasons with p >= {1 - q:.3f}")
print("top-1 accuracy:", round((P_te.argmax(1) == y[te]).mean(), 3),
      "| set coverage:", round(hit.mean(), 3))
print("set size counts:", pd.Series(sets.sum(1)).value_counts().sort_index().to_dict())

q_c = np.array([qhat(s_cal[y[cal] == c]) for c in range(4)])   # Mondrian
sets_m = P_te >= 1 - q_c
hit_m = sets_m[np.arange(len(te)), y[te]]
tab = pd.DataFrame({"reason": np.array(REASONS)[y[te]], "LAC": hit,
                    "Mondrian": hit_m, "size_M": sets_m.sum(1)})
out = tab.groupby("reason").agg(n=("LAC", "size"), LAC=("LAC", "mean"),
                                Mondrian=("Mondrian", "mean"), size_M=("size_M", "mean"))
print(out.round(3)); print("Mondrian thresholds:", np.round(1 - q_c, 3))

q_hat=0.876  -> keep reasons with p >= 0.124
top-1 accuracy: 0.772 | set coverage: 0.91
set size counts: {1: 669, 2: 220, 3: 90, 4: 21}
              n    LAC  Mondrian  size_M
reason                                  
damaged     327  0.957     0.917   2.260
late         82  0.622     0.915   2.220
size        532  0.979     0.906   1.957
wrong_item   59  0.424     0.864   2.288
Mondrian thresholds: [0.391 0.306 0.024 0.034]

৬ · Exchangeability ভাঙলে — সময়, drift আর rolling recalibration

পুরো প্রতিশ্রুতি একটি শর্তের ওপর দাঁড়িয়ে: calibration আর নতুন ডেটা একই রকম। সময়ের ডেটায় এটি প্রায়ই সত্য নয়। বর্ষায় দূরের রাস্তা ডুবে যায়, নতুন কুরিয়ার আসে, ঈদের ভিড় বাড়ে। Input-এর বণ্টন বদলানো (data drift) বা input আর উত্তরের সম্পর্ক বদলানো (concept drift) হলে পুরোনো মাসের score নতুন মাসের ভুলের প্রতিনিধি থাকে না। এই shift-গুলো ধরা আর সামলানোর পূর্ণ আলোচনা পরের পর্ব ৫৯ (Distribution shift — covariate, label আর concept shift ধরা ও সামলানো)-এর বিষয়; এখানে দেখব coverage-এর কী হয়।

নিচের ব্লকে model শেখে মাস ১–৩-এ। “স্থির” পদ্ধতি একবার মাস ৪ দিয়ে calibrate করে আর পরের সব মাসে সেই q̂ রাখে। “Rolling” পদ্ধতি প্রতিটি মাসের শুরুতে আগের মাসের (যার আসল সময় এখন জানা) score দিয়ে আবার calibrate করে; model একই থাকে। মাস ৬ থেকে বর্ষা: ছড়ানো ১.৬ গুণ, আর ৪০ km-এর বেশি দূরের অর্ডারে গড়ে ৩ ঘণ্টা বাড়তি।

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from bazar_conformal import monthly

df = monthly()                       # ৮ মাস × ৩০০ অর্ডার, মাস ৬ থেকে বর্ষা
cols = ["km", "rain"]
old = df[df.month <= 3]
model = LinearRegression().fit(old[cols], old.hours)
df["score"] = np.abs(df.hours - model.predict(df[cols]))
score = lambda m: df.score[df.month == m].to_numpy()

def qhat(s, alpha=0.1):
    k = int(np.ceil((len(s) + 1) * (1 - alpha)))
    return np.sort(s)[k - 1]

q_static = qhat(score(4))            # একবার, মাস ৪ দিয়ে calibrate
rows = []
for m in range(5, 9):
    q_roll = qhat(score(m - 1))      # প্রতি মাসে আগের মাস দিয়ে আবার calibrate
    rows.append([m, np.mean(score(m) <= q_static), np.mean(score(m) <= q_roll),
                 2 * q_static, 2 * q_roll])
cols_out = ["test_month", "static_cov", "rolling_cov", "static_w", "rolling_w"]
print(pd.DataFrame(rows, columns=cols_out).round(3).to_string(index=False))

far = df.km > 40
for name, mask in [("near (<=40 km)", ~far), ("far (>40 km)", far)]:
    s8 = df.score[(df.month == 8) & mask]
    print(f"month 8 {name}: rolling coverage = {np.mean(s8 <= qhat(score(7))):.3f}")

 test_month  static_cov  rolling_cov  static_w  rolling_w
          5       0.867        0.867     9.877      9.877
          6       0.693        0.743     9.877     10.959
          7       0.683        0.857     9.877     18.094
          8       0.670        0.927     9.877     21.648
month 8 near (<=40 km): rolling coverage = 1.000
month 8 far (>40 km): rolling coverage = 0.862

এই সমস্যার জন্য কয়েকটি গবেষণা-পথ আছে। Adaptive conformal inference (ACI; Gibbs আর Candès, ২০২১) প্রতিটি নতুন উত্তরের পর α-কে নিজেই সমন্বয় করে: α_t+1 = α_t + γ(α − err_t), যেখানে err_t = ১ যদি শেষ উত্তর পরিসরের বাইরে পড়ে। বেশি ভুল হলে α ছোট হয়, পরিসর চওড়া হয়; দীর্ঘমেয়াদে গড় ভুল α-র কাছে থাকে, এমনকি বণ্টন বদলালেও। Covariate shift জানা থাকলে weighted conformal (Tibshirani, Barber, Candès আর Ramdas, ২০১৯) calibration score-কে নতুন বণ্টনের অনুপাতে ওজন দেয়। আর time series-এর জন্য পর্ব ৫২-এর backtest-এর ভুল দিয়ে band বানানো আসলে এই rolling ভাবনার কাছাকাছি। কোনোটিই i.i.d.-র মতো শক্ত নিশ্চয়তা দেয় না; সব ক্ষেত্রে coverage-এর নজরদারি চালু রাখুন।

৭ · ডেটা বাঁচানোর রূপ আর MAPIE — একটি সত্যিকারের run

Split conformal-এর দাম ডেটা: calibration ভাগ model শেখায় না। ডেটা কম হলে তিনটি বিকল্প আছে। Full conformal (মূল রূপ) প্রতিটি সম্ভাব্য y-এর জন্য model আবার fit করে; প্রতিশ্রুতি একই, কিন্তু খরচ বিশাল। Jackknife+ (Barber, Candès, Ramdas আর Tibshirani, ২০২১) leave-one-out model আর তাদের residual মিলিয়ে পরিসর বানায়; তার প্রমাণিত নিশ্চয়তা ১ − ২α, যদিও বাস্তবে সাধারণত ১ − α-এর কাছে থাকে। CV+ আর cross-conformal একই কাজ K-fold দিয়ে করে, কম খরচে। এগুলো বোঝার জন্য পর্ব ২১-এর jackknife মনে রাখাই যথেষ্ট।

পাতার Python-এ MAPIE নেই, তাই একটি আলাদা ফোল্ডারে uv দিয়ে ছোট project বানিয়ে আসল run করা হয়েছে। লেখার সময়ে Python 3.13.5, MAPIE 1.6.0, scikit-learn 1.8.0. MAPIE 1.x-এর API আগের 0.x থেকে অনেক বদলেছে (পুরোনো MapieRegressor আর alpha-এর জায়গায় এখন SplitConformalRegressor আর confidence_level), তাই পুরোনো টিউটোরিয়াল বা AI-এর দেওয়া কোড মিলিয়ে নিন।

$ uv init --vcs none --author-from none bazarghar-conformal
$ cd bazarghar-conformal
$ uv add scikit-learn==1.8.0 mapie
 …
 + mapie==1.6.0
 + scikit-learn==1.8.0
 …

from mapie.regression import SplitConformalRegressor
from mapie.classification import SplitConformalClassifier
# একই ডেটা আর একই ভাগ (seed 0 আর 4), model আগে থেকেই fit করা
scr = SplitConformalRegressor(model, confidence_level=0.9, prefit=True)
scr.conformalize(X[cal], y[cal])
pred, iv = scr.predict_interval(X[te])         # iv: (n, 2, 1)
scc = SplitConformalClassifier(clf, confidence_level=0.9,
                               conformity_score="lac", prefit=True)
scc.conformalize(Xc[cal], yc[cal])
_, sets = scc.predict_set(Xc[te])              # sets: (n, 4, 1), True/False

mapie 1.6.0 | scikit-learn 1.8.0
interval shape: (500, 2, 1) | half-width: 7.09
coverage: 0.904
sets shape: (1000, 4) | coverage: 0.91 | size counts: [0, 669, 220, 90, 21]

৮ · AI-এর যুগে — conformal কোড যাচাই

“আমার model-এ ৯০% prediction interval দাও” বললে AI assistant প্রায়ই চলনসই কোড দেয়, কিন্তু এই বিষয়ের ভুলগুলো চুপচাপ থাকে। নিচের ভুলগুলো বারবার দেখা যায়।

AI-এর দেওয়া conformal কোড হাতে পেলে — যাচাইয়ের চেকলিস্ট

• Train, calibration আর test তিনটি আলাদা; calibration-এর কোনো সারি model শেখায়নি, আর সময়ের ডেটা সময় ধরে ভাগ করেছি
• q̂-এর সূচক ⌈(n + ১)(১ − α)⌉ নিজে ছোট একটি উদাহরণে হাতে মিলিয়েছি
• Coverage একটি ভাগে নয়, কয়েকটি seed বা ভাগে মেপেছি, আর গড় প্রস্থ বা সেটের আকারও রিপোর্ট করেছি
• গুরুত্বপূর্ণ দল (এলাকা, শ্রেণি, পণ্যের ধরন) ধরে আলাদা coverage দেখেছি, ছোট দলের ওঠানামা মাথায় রেখে
• Library-র API আর default নিজের version-এর docs-এ মিলিয়েছি, আর ফল হাতে-লেখা কোডের সঙ্গে মিলিয়েছি
• রিপোর্টে “marginal coverage, exchangeability ধরে” লিখেছি, “প্রতিটি ক্রেতার জন্য নিশ্চিত” লিখিনি
• আসল ক্রেতার নাম, ঠিকানা বা ফোন নম্বর AI-তে পেস্ট করিনি