Probability calibration — Platt, isotonic আর reliability diagram

AI-এর ভিত্তি কোর্সের পর্ব ৫৯-এ calibration-এর প্রথম দেখা পেয়েছেন: calibration_curve দিয়ে পাঁচটি ঝুড়ি, আর কড়া penalty দিলে logistic regression-এর সম্ভাবনা গড়ের দিকে চেপে যায়। সেখানে কথা দেওয়া ছিল, সারানোর পদ্ধতি পরের কোর্সে আসবে। এই কোর্সের পর্ব ০২-এ দেখেছেন log loss আর Brier score proper scoring rule, আর hinge loss-এর সেরা উত্তর সম্ভাবনা নয়। এই পর্ব সেই সুতোগুলো এক জায়গায় বাঁধে: সম্ভাবনা কতটা সৎ তা কীভাবে মাপবেন, কেন model ভুল বলে, আর কীভাবে সারাবেন। চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, scikit-learn 1.8.0); প্রথমবার scikit-learn নামতে প্রায় দশ সেকেন্ড লাগে।

১ · ০.৭ বললে কি ৭০%? — calibration আর discrimination

সূত্রে: model প্রতিটি অর্ডারের জন্য একটি সংখ্যা p̂(x) দেয়। Model পুরোপুরি calibrated যদি প্রতিটি p-এর জন্য

P(Y = 1 | p̂(X) = p) = p

এখানে Y আসল ফল (ফেরত এলে ১); বাঁ দিক মানে “যেসব অর্ডারে model p বলেছে, তাদের ফেরতের আসল হার”. শর্তটি একটি অর্ডারের নয়, একই সংখ্যা পাওয়া দলের।

এর পাশে আরেকটি গুণ আছে: discrimination, অর্থাৎ ফেরত-অর্ডারকে বেশি আর না-ফেরত অর্ডারকে কম সংখ্যা দেওয়ার ক্ষমতা। AUC এটিই মাপে (AI-এর ভিত্তি কোর্সের পর্ব ৬০)। দুটি গুণ স্বাধীন, দুটি চরম উদাহরণে তা পরিষ্কার:

• Model ক: প্রতিটি ফেরত-অর্ডারকে ০.৫১, প্রতিটি না-ফেরতকে ০.৪৯ বলে। ক্রম নিখুঁত, AUC = ১. কিন্তু “০.৫১” বলা অর্ডারগুলোর ১০০% ফেরত এসেছে: ভয়ানক ভুল calibrated.
• Model খ: সবাইকে দোকানের গড় ফেরতের হার ০.৩২ বলে। “০.৩২” বলা অর্ডারগুলোর ঠিক ৩২% ফেরত: সংজ্ঞা অনুযায়ী calibrated. কিন্তু AUC = ০.৫, কোনো কাজের নয়।

Score-এর ওপর যেকোনো একঘেয়ে বাড়ন্ত রূপান্তর ক্রম বদলায় না, তাই AUC-ও না। Calibration ঠিক এমন রূপান্তর: ক্রম ভালো করে না, শুধু সংখ্যাগুলোকে সৎ করে।

২ · Reliability diagram — ঝুড়ি ধরে মিলিয়ে দেখা

সংজ্ঞার “p̂ = p” শর্তটি বাস্তবে সরাসরি মাপা যায় না, কারণ ঠিক ০.৬৫০০০ পাওয়া অর্ডার হয়তো একটিও নেই। তাই কাছাকাছি সংখ্যাগুলোকে একটি ঝুড়িতে (bin) রাখা হয়। প্রতিটি ঝুড়ির জন্য দুটি সংখ্যা: গড় অনুমান (x অক্ষ) আর বাস্তবে ফেরতের হার (y অক্ষ)। বিন্দুগুলো কর্ণরেখা y = x-এর ওপর থাকলে model calibrated. নিচে একটি ঝুড়ি ভরার পুরো হিসাব:

ঝুড়ি কীভাবে কাটবেন — সমান প্রস্থ না সমান সংখ্যা

scikit-learn-এর calibration_curve(y_true, y_prob, n_bins=5, strategy="uniform") দুই রকম ঝুড়ি দেয়। "uniform" মানে [০, ০.১), [০.১, ০.২) … সমান প্রস্থ; পড়া সহজ, কিন্তু অনুমান কম পড়া অঞ্চলের ঝুড়িতে দুই-চারটি অর্ডার, তাই বিন্দু লাফায়। "quantile" মানে প্রতিটি ঝুড়িতে প্রায় সমান অর্ডার; বিন্দুগুলো সমান নির্ভরযোগ্য, কিন্তু অনুমান ঘন হলে গাদাগাদি। খালি ঝুড়ি বাদ পড়ে। দুটিই দেখুন, সঙ্গে অনুমানের একটি histogram: কোন অঞ্চলে কত অর্ডার, তা না জানলে বাঁকা বিন্দুর গুরুত্ব বোঝা যায় না।

৩ · ECE, MCE আর Brier score — calibration এক সংখ্যায়

ছবি পড়তে ভালো, কিন্তু অনেক model তুলনা বা রিপোর্টের জন্য একটি সংখ্যা দরকার। সবচেয়ে প্রচলিত হলো expected calibration error (ECE): প্রতিটি ঝুড়ির ফাঁক, ঝুড়ির অর্ডারের অনুপাত দিয়ে ওজন করে যোগ।

ECE = ∑_b (n_b / N) · | ō_b − p̄_b |

এখানে b একটি ঝুড়ি, n_b তার অর্ডারসংখ্যা, N মোট অর্ডার, p̄_b ঝুড়ির গড় অনুমান আর ō_b বাস্তব হার। সবচেয়ে বড় ফাঁকটিকে বলে MCE (maximum calibration error). এই দুটি নাম জনপ্রিয় হয় Mahdi Pakdaman Naeini, Gregory Cooper আর Milos Hauskrecht-এর ২০১৫ সালের AAAI প্রবন্ধ থেকে।

Brier score — ক্রম আর সততা একসঙ্গে

Glenn Brier ১৯৫০ সালে Monthly Weather Review-এ বৃষ্টির পূর্বাভাস যাচাইয়ের জন্য যে মাপ প্রস্তাব করেন, binary ক্ষেত্রে তা হলো:

BS = (1/N) ∑_i (p_i − y_i)²

যেখানে p_i অনুমান আর y_i ০ বা ১. ০.৮ বলে ফেরত এলে (০.৮ − ১)² = ০.০৪, না এলে (০.৮ − ০)² = ০.৬৪. এটি log loss-এর মতোই proper scoring rule (পর্ব ০২): আসল সম্ভাবনা সৎভাবে বললেই প্রত্যাশিত Brier সবচেয়ে কম। তবে log loss আত্মবিশ্বাসী ভুলকে অসীম পর্যন্ত শাস্তি দিতে পারে, Brier সর্বোচ্চ ১. ঝুড়ি লাগে না বলে binning-এর পছন্দও ঢোকে না। তাই model তুলনায় Brier বা log loss, রোগনির্ণয়ে reliability diagram আর ECE.

৪ · চার model, একই ফেরত-ডেটা — reliability diagram হাতে-কলমে

এই পর্বের ডেটা: বাজারঘরের ৫,০০০টি কাল্পনিক অর্ডার, seed দিয়ে বানানো, কোনো আসল দোকানের নয়। ছয়টি feature (ডেলিভারির দিন, দূরত্ব, courier-এর দেরির সূচক, দাম, পোশাক কি না, প্রথম অর্ডার কি না), লক্ষ্য ফেরত আসবে কি না। দূরত্ব আর দেরির সূচক ইচ্ছা করে ডেলিভারির দিনের সঙ্গে জড়ানো, আর ফেরতের আসল নিয়মে তারা নেই। ভাগ: ২,৫০০ train, ১,০০০ calibration, ১,৫০০ তালাবন্ধ test. চারটি ব্লকই এই module থেকে ডেটা নেয়।

# bazar_returns.py — বাজারঘরের কাল্পনিক অর্ডার আর ফেরত (seed দিয়ে বানানো)
import numpy as np

def make_returns(n=5000, seed=26, rate_shift=0.0):
    rng = np.random.default_rng(seed)
    days = rng.gamma(3, 1.2, n)                       # ডেলিভারিতে দিন
    dist = 0.8 * days + rng.normal(0, 0.6, n)         # hub থেকে দূরত্ব (১০০ km), days-এর সঙ্গে জড়ানো
    late = days + rng.normal(0, 0.8, n)               # courier-এর দেরির সূচক, এটিও জড়ানো
    price = rng.lognormal(0.0, 0.6, n)                # দাম, হাজার টাকা
    fashion = rng.integers(0, 2, n)                   # পোশাক বা জুতা?
    first = rng.integers(0, 2, n)                     # গ্রাহকের প্রথম অর্ডার?
    z = (-4.0 + rate_shift + 0.45 * days + 0.3 * price + 1.1 * fashion
         + 0.7 * first + 0.6 * fashion * first)
    y = (rng.random(n) < 1 / (1 + np.exp(-z))).astype(int)
    X = np.column_stack([days, dist, late, price, fashion, first])
    return X, y

def split(X, y):                                      # train ২৫০০ · calibration ১০০০ · test ১৫০০
    return (X[:2500], y[:2500]), (X[2500:3500], y[2500:3500]), (X[3500:], y[3500:])

def ece(p, y, bins=10):                               # সমান-প্রস্থে ECE (ব্লক ১-এর একই function)
    idx = np.minimum((p * bins).astype(int), bins - 1)
    return sum(np.mean(idx == b) * abs(p[idx == b].mean() - y[idx == b].mean())
               for b in range(bins) if np.any(idx == b))

প্রথম ব্লক চারটি model শেখায়: logistic regression, Gaussian Naive Bayes, ২০০ গাছের random forest, আর RBF SVM. SVM-এর decision function-কে সরাসরি sigmoid দিয়ে ০–১-এ আনা হয়েছে, যা অনেকে না ভেবেই করেন। তারপর দুই রকম ঝুড়িতে reliability diagram, আর নিজে লেখা ECE ও Brier.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.calibration import calibration_curve
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.naive_bayes import GaussianNB
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC
from bazar_returns import make_returns, split

X, y = make_returns()
(Xtr, ytr), _, (Xte, yte) = split(X, y)
def ece(p, y, bins=10):                      # নিজে লেখা ECE, সমান-প্রস্থের ঝুড়ি
    idx = np.minimum((p * bins).astype(int), bins - 1)
    return sum(np.mean(idx == b) * abs(p[idx == b].mean() - y[idx == b].mean())
               for b in range(bins) if np.any(idx == b))

models = {"logistic": LogisticRegression(max_iter=1000), "gaussian NB": GaussianNB(),
          "forest": RandomForestClassifier(200, random_state=0)}
probs = {k: m.fit(Xtr, ytr).predict_proba(Xte)[:, 1] for k, m in models.items()}
svm = make_pipeline(StandardScaler(), SVC()).fit(Xtr, ytr)
probs["SVM sigmoid"] = 1 / (1 + np.exp(-svm.decision_function(Xte)))  # score-এ সরাসরি sigmoid

fig, axes = plt.subplots(1, 2, figsize=(9, 4.2), sharey=True)
for ax, strat in zip(axes, ["uniform", "quantile"]):
    ax.plot([0, 1], [0, 1], "k--", lw=1, label="perfect")
    for name, p in probs.items():
        frac, mean = calibration_curve(yte, p, n_bins=10, strategy=strat)
        ax.plot(mean, frac, "o-", ms=3, label=name)
    ax.set(title=f"strategy='{strat}'", xlabel="mean predicted probability")
axes[0].set_ylabel("observed return rate"); axes[0].legend(fontsize=8)
print(f"test orders = {len(yte)}, returned = {yte.mean():.3f}\nmodel         Brier   ECE")
for name, p in probs.items():
    print(f"{name:12s}  {np.mean((p - yte) ** 2):.4f}  {ece(p, yte):.4f}")

test orders = 1500, returned = 0.321
model         Brier   ECE
logistic      0.1668  0.0248
gaussian NB   0.2117  0.1624
forest        0.1768  0.0306
SVM sigmoid   0.1806  0.0872

৫ · Model কেন ভুল calibrated হয়

Calibration নির্ভর করে model কোন loss কমিয়েছে, কী ধরে নিয়েছে, আর ডেটা কতটা মিলেছে তার ওপর। Alexandru Niculescu-Mizil আর Rich Caruana ২০০৫ সালের ICML প্রবন্ধে (“Predicting good probabilities with supervised learning”) দশটির মতো algorithm অনেক ডেটাসেটে তুলনা করে চেনা আকারগুলো দেখান। মূল কারণগুলো:

Naive Bayes-এর হিসাবটি দেখুন। AI-এর ভিত্তি কোর্সের পর্ব ৫১-এ দেখেছেন, Naive Bayes প্রতিটি feature-এর likelihood ratio গুণ করে। ধরুন “ডেলিভারিতে ৮ দিন” একাই ফেরতের odds ৩ গুণ করে। আমাদের ডেটায় দূরত্ব আর দেরির সূচক প্রায় একই তথ্য বহন করে, তাই তারাও প্রত্যেকে প্রায় ৩ গুণ বলে। Naive Bayes তিনটিকে আলাদা প্রমাণ ধরে odds ২৭ গুণ করে, অথচ আসল প্রমাণ ৩ গুণের। ফল: odds-এর log মাপে প্রায় তিন গুণ চরম অনুমান। ক্রম তুলনামূলক কম ক্ষতিগ্রস্ত হয়।

SVM-এর গাদাগাদি। SVM-এর decision function f(x) margin-এর ভাষায় মাপা: যে বিন্দু সঠিক দিকে margin-এর বাইরে, তার |f| ≥ ১, আর soft-margin SVM অনেক বিন্দুকে ঠিক f ≈ ±১-এ বসায় (পর্ব ২৮, Soft margin — hinge loss, C আর dual problem-এ এর কারণ আসবে)। আমাদের run-এ test অর্ডারের score-এর মধ্যমা প্রায় −১.০০, আর σ(−১) = ১/(১ + e) ≈ ০.২৬৯. তাই প্রায় অর্ধেক অর্ডার একই “০.২৭” পায়, যদিও তাদের আসল হার এর চেয়ে কম। Score-এর মাপ C, kernel আর scaling-এর ওপর নির্ভর করে; σ(f)-এর কোনো সম্ভাবনা-অর্থ নেই।

৬ · Brier score ভাঙা — reliability, resolution আর uncertainty

Allan Murphy ১৯৭৩ সালে (Journal of Applied Meteorology) দেখান, ঝুড়ি ধরে হিসাব করলে Brier score তিনটি অর্থপূর্ণ অংশে ভাঙে। এই Brier decomposition বলে model-এর ক্ষতি কোথা থেকে আসছে:

BS ≈ REL − RES + UNC

• REL (reliability) = ∑_b (n_b/N)(p̄_b − ō_b)²: ঝুড়ির গড় অনুমান আর বাস্তব হারের ফাঁকের বর্গ। এটিই calibration-এর ভুল; কম হলে ভালো। ECE-র সঙ্গে পার্থক্য শুধু পরম মানের বদলে বর্গ।
• RES (resolution) = ∑_b (n_b/N)(ō_b − ō)²: ঝুড়িগুলোর বাস্তব হার সামগ্রিক হার ō থেকে কত দূরে। Model অর্ডারগুলোকে সত্যিই আলাদা দলে ভাগ করতে পারলে বড়; বেশি হলে ভালো, তাই বিয়োগ।
• UNC (uncertainty) = ō(১ − ō): সমস্যাটি নিজে কতটা অনিশ্চিত। Model এটি বদলাতে পারে না; ō = ০.৫ হলে সর্বোচ্চ ০.২৫.

Model খ সবাইকে ō বলে: একটিই ঝুড়ি, REL = RES = ০, তাই BS = UNC. এজন্য Brier skill score = ১ − BS/UNC মাপে “সবাইকে গড় হার বলা”-র চেয়ে কত ভালো। মূল শিক্ষা: calibration শুধু REL কমায়; RES বাড়াতে লাগে ভালো feature আর model.

“≈” কেন? ঝুড়ির ভেতরে অনুমানগুলো হুবহু এক না হলে একটি ছোট অতিরিক্ত পদ থাকে। এক ঝুড়ির ভেতরে বর্গ ভেঙে লিখলে হুবহু অভেদটি দাঁড়ায় BS = REL − RES + UNC + (WBV − ২·WBC), যেখানে WBV = (1/N)∑(p_i − p̄_b)² ঝুড়ির ভেতরের অনুমানের ছড়ানো, আর WBC = (1/N)∑(p_i − p̄_b)(y_i − ō_b) ঝুড়ির ভেতরে অনুমান আর ফলের সহ-ভেদ। ব্লকটি তিনটি অংশ আর এই অতিরিক্ত পদ হিসাব করে মিলিয়ে দেখায়।

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.naive_bayes import GaussianNB
from bazar_returns import make_returns, split

X, y = make_returns()
(Xtr, ytr), _, (Xte, o) = split(X, y)

def murphy(f, o, bins=10):                 # Brier-কে তিন ভাগে ভাঙা, সমান-প্রস্থের ঝুড়ি
    k = np.minimum((f * bins).astype(int), bins - 1)
    N, obar = len(f), o.mean()
    rel = res = within = 0.0
    for b in np.unique(k):
        m = k == b
        fk, ok = f[m].mean(), o[m].mean()
        rel += m.sum() * (fk - ok) ** 2 / N          # reliability: অনুমান বনাম বাস্তব
        res += m.sum() * (ok - obar) ** 2 / N        # resolution: ঝুড়ির হার গড় থেকে কত দূরে
        within += np.sum((f[m] - fk) ** 2 - 2 * (f[m] - fk) * (o[m] - ok)) / N
    return rel, res, obar * (1 - obar), within

models = {"logistic": LogisticRegression(max_iter=1000), "gaussian NB": GaussianNB()}
for name, mdl in models.items():
    f = mdl.fit(Xtr, ytr).predict_proba(Xte)[:, 1]
    rel, res, unc, within = murphy(f, o)
    bs = np.mean((f - o) ** 2)
    print(f"{name}: Brier = {bs:.4f}, skill vs base rate = {1 - bs / unc:.3f}")
    print(f"  REL {rel:.4f} - RES {res:.4f} + UNC {unc:.4f} = {rel - res + unc:.4f}")
    print(f"  + within-bin {within:+.4f} = {rel - res + unc + within:.4f}")
base = np.full(len(o), o.mean())                      # সবাইকে একই গড় হার বলা
rel, res, unc, _ = murphy(base, o)
print(f"always {o.mean():.3f}: REL {rel:.4f}, RES {res:.4f}, Brier {np.mean((base - o) ** 2):.4f}")

logistic: Brier = 0.1668, skill vs base rate = 0.235
  REL 0.0013 - RES 0.0518 + UNC 0.2181 = 0.1676
  + within-bin -0.0007 = 0.1668
gaussian NB: Brier = 0.2117, skill vs base rate = 0.029
  REL 0.0347 - RES 0.0394 + UNC 0.2181 = 0.2134
  + within-bin -0.0017 = 0.2117
always 0.321: REL 0.0000, RES 0.0000, Brier 0.2181

৭ · Platt scaling আর isotonic regression — score থেকে সৎ সম্ভাবনা

সারানোর মূল ভাবনা একটাই: model-এর score s রেখে দিন, শুধু তার ওপর একটি একঘেয়ে function g শিখুন যাতে g(s) সৎ সম্ভাবনা হয়। g শেখানো হয় এমন ডেটায় যা model দেখেনি, কারণ train ডেটায় model-এর score নিজেই আশাবাদী। একে বলে calibration set.

Platt scaling — দুটি সংখ্যার sigmoid

John Platt ১৯৯৯ সালে SVM-এর score-কে সম্ভাবনায় আনতে প্রস্তাব করেন (Advances in Large Margin Classifiers বইয়ের অধ্যায়, ২০০০ সালে প্রকাশিত):

p = σ(a · s + b) = 1 / (1 + e^{−(a·s + b)})

a আর b শেখা হয় calibration set-এ log loss কমিয়ে; অর্থাৎ score-কে একমাত্র feature ধরে একটি এক-মাত্রার logistic regression. a score-কে প্রসারিত বা সংকুচিত করে (a < ১ মানে model অতি-আত্মবিশ্বাসী ছিল), b সরায়। ছোট ডেটায় sigmoid খুব খাড়া না হতে Platt লক্ষ্য ১ আর ০-এর বদলে (N₊ + ১)/(N₊ + ২) আর ১/(N₋ + ২) নেন; scikit-learn-ও তাই করে।

দুটি সূক্ষ্ম বিষয়। এক, s কোনটি? SVM-এ decision function; যে model নিজেই সম্ভাবনা দেয়, সেখানে log-odds (logit) log(p/(১ − p)) স্বাভাবিক, কারণ sigmoid তার উল্টো। কিন্তু CalibratedClassifierCV decision function না পেলে কাঁচা predict_proba-কেই s নেয় (1.8.0-এর docs আর source). দুই, Platt ধরে নেয় বিকৃতিটি sigmoid-আকারের; SVM বা boosting-এর S-আকারে তা মেলে, অন্য আকারে পুরোটা সারে না।

Isotonic regression — শুধু “একঘেয়ে” শর্ত

Bianca Zadrozny আর Charles Elkan ২০০২ সালে (KDD) প্রস্তাব করেন: g-এর কোনো আকার ধরে নেবেন না, শুধু চান g একঘেয়ে বাড়ন্ত হোক। Calibration set-এ (s_i, y_i) জোড়াগুলোয় এমন একঘেয়ে সিঁড়ি খুঁজুন যার বর্গ-ভুল সবচেয়ে কম। উত্তর দেয় pool adjacent violators (PAV) algorithm (Ayer ও সহকর্মীরা, ১৯৫৫): score অনুযায়ী সাজিয়ে বাঁ থেকে ডানে যান, কোথাও ক্রম ভাঙলে পাশের মানগুলো মিলিয়ে গড় করুন, যতক্ষণ না সব একঘেয়ে হয়। নিচে ৮টি বিন্দুতে পুরো হিসাব, আর একই ৮টিতে Platt:

অনেক শ্রেণি হলে — temperature scaling

Platt আর isotonic দুই শ্রেণির জন্য; অনেক শ্রেণিতে scikit-learn প্রতিটি শ্রেণির one-vs-rest calibrator শেখায়, তারপর যোগফল ১ করে (পর্ব ২৫, Multiclass আর multilabel — OvR, OvO, softmax আর classifier chain). Chuan Guo ও সহকর্মীরা ২০১৭ সালের ICML-এ সরলতর পথ জনপ্রিয় করেন: temperature scaling, অর্থাৎ p = softmax(z/T), সব logit একটি সংখ্যা T দিয়ে ভাগ। T > ১ সম্ভাবনা নরম করে, কিন্তু সবচেয়ে বড় logit বড়ই থাকে, তাই accuracy বদলায় না। scikit-learn 1.8-এ method="temperature" যোগ হয়েছে। Scikit-learn-এর সঙ্গে আসা digits ডেটায় (১০ শ্রেণি) প্রায় শাস্তিহীন logistic regression, আলাদা ভাগে T; ডেস্কটপে সত্যিকারের run (Python 3.14.4, scikit-learn 1.8.0, একটি run):

~/bazarghar-ml$ python temp_scaling.py
scikit-learn 1.8.0
raw          accuracy 0.963  log loss 0.428  mean max-prob 0.993
temperature  accuracy 0.963  log loss 0.165  mean max-prob 0.940
inverse temperature beta = 0.206, so T = 4.84

Accuracy হুবহু এক (০.৯৬৩), কিন্তু log loss ০.৪২৮ থেকে ০.১৬৫. কাঁচা model গড়ে ৯৯.৩% নিশ্চয়তা দেখাচ্ছিল অথচ ঠিক ছিল ৯৬.৩%; T ≈ ৪.৮-এ ভাগের পর ৯৪.০%. ভেতরে রাখা beta_ = ১/T একটি private attribute, যার নাম বদলাতে পারে।

৮ · হাতে-কলমে — আলাদা calibration ডেটা, CalibratedClassifierCV আর ফাঁদ

CalibratedClassifierCV(estimator, method=..., cv=5)-এ ডিফল্ট ensemble="auto" এখানে True: প্রতিটি fold-এ চার ভাগে model-এর একটি কপি, বাকি ভাগে তার calibrator; শেষে পাঁচ জোড়ার সম্ভাবনার গড়। ফলে আলাদা calibration ভাগ লাগে না (AI-এর ভিত্তি কোর্সের পর্ব ৬৬-এর cross-validation-এর মতো)। আগে শেখানো model হলে FrozenEstimator-এ মুড়ে দিন; তখন দেওয়া সব ডেটা শুধু calibration-এ যায়, আর ভাগ আলাদা রাখার দায় আপনার। পুরোনো লেখায় (আর AI-এর কোডে) যে cv="prefit" দেখবেন, তা 1.8.0-এ আর নেই; ডেস্কটপে সত্যিকারের run:

~/bazarghar-ml$ python prefit.py
…
sklearn.utils._param_validation.InvalidParameterError: The 'cv' parameter of CalibratedClassifierCV must be an int in the range [2, inf), an object implementing 'split' and 'get_n_splits', an iterable or None. Got 'prefit' instead.

নিচের ব্লক Naive Bayes-কে চারভাবে সারায়: নিজে লেখা Platt (logit score-এর ওপর এক-মাত্রার logistic), নিজে লেখা isotonic (দুটিই ১,০০০টি calibration অর্ডারে), আর CalibratedClassifierCV-এর sigmoid ও isotonic (train আর calibration মিলিয়ে ৩,৫০০টিতে ৫-fold). তারপর ফাঁদটি: forest-এর calibrator তার নিজের train অর্ডারেই শেখানো।

import numpy as np
from sklearn.calibration import CalibratedClassifierCV
from sklearn.ensemble import RandomForestClassifier
from sklearn.isotonic import IsotonicRegression
from sklearn.linear_model import LogisticRegression
from sklearn.naive_bayes import GaussianNB
from bazar_returns import make_returns, split, ece

X, y = make_returns()
(Xtr, ytr), (Xc, yc), (Xte, yte) = split(X, y)
nb = GaussianNB().fit(Xtr, ytr)                    # শুধু train অংশে শেখা
logit = lambda p: np.log(np.clip(p, 1e-12, 1) / np.clip(1 - p, 1e-12, 1))
s_c, s_te = (logit(nb.predict_proba(A)[:, 1]) for A in (Xc, Xte))

platt = LogisticRegression(C=1e10).fit(s_c[:, None], yc)   # ১-মাত্রার logistic = Platt
a, b = platt.coef_[0, 0], platt.intercept_[0]
iso = IsotonicRegression(out_of_bounds="clip").fit(s_c, yc)
res = {"raw NB": 1 / (1 + np.exp(-s_te)),
       "Platt (own)": platt.predict_proba(s_te[:, None])[:, 1],
       "isotonic (own)": iso.predict(s_te)}
for m in ("sigmoid", "isotonic"):                  # train + calibration একসঙ্গে, ৫-fold
    cc = CalibratedClassifierCV(GaussianNB(), method=m, cv=5)
    res[f"CCV {m}"] = cc.fit(np.r_[Xtr, Xc], np.r_[ytr, yc]).predict_proba(Xte)[:, 1]
print(f"Platt: p = sigmoid({a:.3f} * score {b:+.3f}); isotonic knots = {len(iso.X_thresholds_)}")
print("on test             Brier   ECE")
for k, p in res.items():
    print(f"{k:16s}  {np.mean((p - yte) ** 2):.4f}  {ece(p, yte):.4f}")

rf = RandomForestClassifier(200, random_state=0).fit(Xtr, ytr)  # ফাঁদ: train-এর ওপরেই calibration
for name, (A, t) in {"on train (leak)": (Xtr, ytr), "on held-out": (Xc, yc)}.items():
    cal = IsotonicRegression(out_of_bounds="clip").fit(rf.predict_proba(A)[:, 1], t)
    p = cal.predict(rf.predict_proba(Xte)[:, 1])
    print(f"forest isotonic {name:16s} Brier {np.mean((p - yte) ** 2):.4f}  ECE {ece(p, yte):.4f}")
p = rf.predict_proba(Xte)[:, 1]
print(f"forest raw                       Brier {np.mean((p - yte) ** 2):.4f}  ECE {ece(p, yte):.4f}")

Platt: p = sigmoid(0.356 * score -0.395); isotonic knots = 30
on test             Brier   ECE
raw NB            0.2117  0.1624
Platt (own)       0.1880  0.0870
isotonic (own)    0.1787  0.0453
CCV sigmoid       0.1872  0.1005
CCV isotonic      0.1785  0.0429
forest isotonic on train (leak)  Brier 0.2332  ECE 0.2221
forest isotonic on held-out      Brier 0.1768  ECE 0.0325
forest raw                       Brier 0.1768  ECE 0.0306

৯ · Calibration কেন দরকার — খরচের সিদ্ধান্ত আর prior বদল

পর্ব ০২-এর অংশ ৭ আর AI-এর ভিত্তি কোর্সের পর্ব ৫৪-এ দেখেছেন, সিদ্ধান্ত নিতে হয় expected cost দেখে। বাজারঘর ঝুঁকির অর্ডারে ডেলিভারির আগে ফোন করে অগ্রিম চার্জ চায়। কাল্পনিক সংখ্যায়: প্রতিটি ফোনে খরচ ৬০ টাকা, আর ফোন করলে ফেরত আটকায়; ফোন না করলে ফেরত এলে ৩০০ টাকা ক্ষতি (দুই দিকের courier আর আবার প্যাক)। ফেরতের সম্ভাবনা p হলে:

ফোন করলে খরচ ৬০  ·  না করলে প্রত্যাশিত খরচ ৩০০ · p  ⇒  ফোন করুন যদি p > ৬০/৩০০ = ০.২

সাধারণ রূপে, চারটি ঘরের cost matrix-এ threshold t* = (c_FP − c_TN) / [(c_FP − c_TN) + (c_FN − c_TP)]; এখানে c_FP = c_TP = ৬০, c_FN = ৩০০, c_TN = ০, তাই ৬০/(৬০ + ২৪০) = ০.২. Charles Elkan ২০০১ সালের IJCAI প্রবন্ধে (“The foundations of cost-sensitive learning”) এটি গুছিয়ে লেখেন। লুকানো শর্ত: p সৎ হতে হবে।

Prior বদলালে — odds দিয়ে সংশোধন

ঈদের মৌসুমে পোশাকের অর্ডার বাড়ে, আর ফেরতের সামগ্রিক হার π = ০.৩১ থেকে ধরা যাক π′ = ০.৫০ হয়। যদি ফেরত-অর্ডার আর না-ফেরত অর্ডারগুলো দেখতে আগের মতোই থাকে, অর্থাৎ P(x | y) না বদলায় (একে বলে label shift বা prior shift), তবে Bayes-এর সূত্র থেকে নতুন সম্ভাবনা আসে model আবার না শিখিয়েই:

p′ = p·(π′/π) / [ p·(π′/π) + (1 − p)·(1 − π′)/(1 − π) ]

অর্থাৎ odds′ = odds × [π′/(১ − π′)] / [π/(১ − π)]. Marco Saerens ও সহকর্মীরা ২০০২ সালে (Neural Computation) এটি আর π′ অজানা হলে EM দিয়ে অনুমানের পদ্ধতি দেন। সূত্রটি ধরে নেয় p আগে থেকেই calibrated. নিচের ব্লক দুটি কাজই দেখায়।

import numpy as np
from sklearn.calibration import CalibratedClassifierCV
from sklearn.linear_model import LogisticRegression
from sklearn.naive_bayes import GaussianNB
from bazar_returns import make_returns, split, ece

X, y = make_returns()
(Xtr, ytr), (Xc, yc), (Xte, yte) = split(X, y)
CALL, RETURN = 60, 300                     # টাকা: ফোন করে নিশ্চিত করা, আর ফেরতের ক্ষতি
t = CALL / RETURN                          # p > 0.2 হলে ফোন করাই সস্তা
raw = GaussianNB().fit(Xtr, ytr).predict_proba(Xte)[:, 1]
cal = CalibratedClassifierCV(GaussianNB(), method="isotonic", cv=5).fit(
    np.r_[Xtr, Xc], np.r_[ytr, yc]).predict_proba(Xte)[:, 1]
for name, p in {"raw NB": raw, "calibrated NB": cal}.items():
    act = p > t
    real = np.where(act, CALL, RETURN * yte).sum()          # আসলে কত খরচ হলো
    plan = np.where(act, CALL, RETURN * p).sum()            # model নিজে যা ভেবেছিল
    print(f"{name:14s} calls {act.sum():4d}  planned {plan:7.0f}  real {real:7.0f} taka")
print(f"call nobody: {RETURN * yte.sum()} taka, call everyone: {CALL * len(yte)} taka")

# ঈদের মৌসুম: ফেরতের হার বাড়ে, কিন্তু ফেরত-অর্ডারগুলোর চেহারা একই (label shift)
Xp, yp = make_returns(n=20000, seed=99)
rng = np.random.default_rng(1)
pos, neg = rng.permutation(np.flatnonzero(yp == 1)), rng.permutation(np.flatnonzero(yp == 0))
idx = np.r_[pos[:750], neg[:750]]                          # নতুন হার 0.50
lr = LogisticRegression(max_iter=1000).fit(Xtr, ytr)
p_old = lr.predict_proba(Xp[idx])[:, 1]
pi_old, pi_new = ytr.mean(), 0.50
r1, r0 = pi_new / pi_old, (1 - pi_new) / (1 - pi_old)
p_new = p_old * r1 / (p_old * r1 + (1 - p_old) * r0)       # odds-কে prior-এর অনুপাতে গুণ
print(f"train rate {pi_old:.3f} -> Eid rate {yp[idx].mean():.3f}")
for name, p in {"as trained": p_old, "prior-adjusted": p_new}.items():
    print(f"{name:15s} mean p {p.mean():.3f}  Brier {np.mean((p - yp[idx]) ** 2):.4f}"
          f"  ECE {ece(p, yp[idx]):.4f}")

raw NB         calls  481  planned   49979  real   91260 taka
calibrated NB  calls  864  planned   71350  real   71940 taka
call nobody: 144600 taka, call everyone: 90000 taka
train rate 0.309 -> Eid rate 0.500
as trained      mean p 0.356  Brier 0.2023  ECE 0.1444
prior-adjusted  mean p 0.499  Brier 0.1777  ECE 0.0242

১০ · AI-এর যুগে — সম্ভাবনার কোড আর দাবি যাচাই

AI assistant-এর calibration কোড চলে আর সুন্দর সংখ্যা ছাপে। বিপদ: ভুল জায়গায় শেখানো calibrator-ও চমৎকার ECE দেখায়, যদি একই ডেটায় মাপা হয়। পরিচিত ভুলগুলো:

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

• Calibrator এমন ডেটায় শেখানো, যা model দেখেনি; আর মাপা হয়েছে তৃতীয় একটি তালাবন্ধ test-এ
• Reliability diagram দেখেছি, ঝুড়ির অর্ডারসংখ্যাসহ; শুধু একটি ECE সংখ্যা নয়
• Brier বা log loss আগে-পরে তুলনা করেছি, আর AUC বদলায়নি তা মিলিয়েছি
• Calibration set-এর আকার দেখে Platt বা isotonic বেছেছি, কয়েকটি seed-এ
• Resampling বা class weight থাকলে আসল হারের ডেটায় calibration করেছি
• Threshold-এর সূত্র একটি ছোট সংখ্যার উদাহরণে নিজে হিসাব করে মিলিয়েছি
• Function আর argument-এর নাম নিজের version-এর docs-এ মিলিয়েছি, আর কোনো আসল গ্রাহকের ডেটা AI-তে পেস্ট করিনি