Model তুলনার পরিসংখ্যান — paired test, McNemar আর একাধিক তুলনা

AI-এর ভিত্তি কোর্সের পর্ব ৩৭-এ hypothesis test, p-value, t-test আর Bonferroni-র প্রথম পরিচয় পেয়েছেন, আর পর্ব ৬৬-এ repeated আর nested CV. এই কোর্সের পর্ব ০৪-এ seed-ভেদে ছড়ানো, mean ± sd আর Nadeau–Bengio-র corrected resampled t-test দেখেছেন, সেখানে প্রতিশ্রুতি ছিল তুলনার পূর্ণ পরিসংখ্যান এখানে আসবে। এই পর্ব সেই কাজটি করে: কোন প্রশ্নে কোন test, test-টি নিজেই কতটা সৎ, আর অনেকগুলো তুলনা একসঙ্গে করলে কী হয়। চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, SciPy, scikit-learn আর statsmodels); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে, scikit-learn-এ প্রায় দশ সেকেন্ড।

১ · “০.৮৬ বনাম ০.৮৪” কেন যথেষ্ট নয় — আগে প্রশ্নটি ঠিক করা

ML-এ “model A-র accuracy ০.৮৬, B-র ০.৮৪” লাইনটি প্রতিদিন লেখা হয়। এর পেছনে অন্তত তিনটি এলোমেলোতার উৎস আছে: test set (কোন সারিগুলো মাপার জন্য পড়েছে), training set (কোন সারিগুলো দিয়ে শেখানো হয়েছে), আর algorithm-এর ভেতরের এলোমেলোতা (seed, initialization). কোন উৎস হিসাবে ধরবেন, তা নির্ভর করে প্রশ্নের ওপর। Thomas Dietterich ১৯৯৮ সালে Neural Computation-এ প্রশ্নগুলোকে মোটামুটি এভাবে আলাদা করেছিলেন:

প্রতিটি test-এর কাঠামো একই, যা পর্ব ৩৭-এ দেখেছেন: একটি null hypothesis (দুজন সমান ভালো), একটি test statistic, আর p-value = null সত্য হলে এত বড় বা আরও বড় statistic আসার সম্ভাবনা। একটি test “সৎ” যদি null সত্য হলে তার p < ০.০৫ হওয়ার হার সত্যিই প্রায় ৫% হয়; এই হারকে বলে Type I error rate. আর আসল পার্থক্য ধরার সম্ভাবনা হলো test-এর power.

২ · একই test set-এ দুই classifier — McNemar-এর 2×2 টেবিল

দুটি classifier A আর B একই n-টি test সারিতে অনুমান করেছে। প্রতিটি সারি চারটি ঘরের একটিতে পড়ে, আর এই ঘরগুলো একটি contingency table:

A-র accuracy (a + b)/n, B-র (a + c)/n, তাই পার্থক্য (b − c)/n: a আর d বাদ পড়ে যায়। b আর c-কে বলে discordant (মতভেদের) সারি। Null hypothesis হলো দুজনের ভুলের হার সমান; সেটি সত্য হলে, মতভেদের কোনো সারিতে “A ঠিক” হওয়ার সম্ভাবনা ঠিক ½. অর্থাৎ b + c-টি মতভেদ দেওয়া থাকলে b একটি binomial(b + c, ½) চলক। এখান থেকে দুটি রূপ আসে:

Exact: p = 2 · P(X ≥ max(b, c)), X ~ Binomial(b + c, ½)    χ² রূপ: χ² = (|b − c| − 1)² / (b + c), df = 1

χ² রূপটি binomial-এর normal approximation-এর বর্গ: b − c-এর গড় ০ আর variance b + c, তাই (b − c)/√(b + c) প্রায় standard normal, আর তার বর্গ ১ df-এর χ². “−১” হলো continuity correction: পূর্ণসংখ্যার গণনাকে মসৃণ বক্ররেখা দিয়ে মাপলে যে আশাবাদী ঝোঁক আসে, তা কমানো (A. L. Edwards ১৯৪৮ সালে এই সংশোধন প্রস্তাব করেন বলে প্রচলিত)। Quinn McNemar পরীক্ষাটি ১৯৪৭ সালে Psychometrika-তে প্রকাশ করেন, মনোবিজ্ঞানে একই মানুষের আগে-পরে উত্তরের তুলনার জন্য। b + c ছোট হলে (অনেক বইয়ে ২৫-এর নিচে বলা হয়, এটি একটি প্রথা) exact রূপ নিন।

৩ · McNemar হাতে-কলমে — χ², exact binomial আর statsmodels

বাজারঘরের ১,৫০০টি কাল্পনিক অর্ডার: ডেলিভারিতে কত দিন, দাম, আর এটি গ্রাহকের প্রথম অর্ডার কি না। লক্ষ্য, অর্ডারটি ফেরত আসবে কি না। প্রথম ১,০০০টিতে দুটি model train হয়, A = logistic regression, B = গভীরতা ৪-এর decision tree; শেষ ৫০০টি তালাবন্ধ test set. ব্লকটি টেবিল বানায়, তারপর তিনভাবে p-value বের করে: নিজের সূত্রে χ², SciPy-র binomtest দিয়ে exact, আর statsmodels-এর mcnemar দিয়ে যাচাই।

import numpy as np
from scipy import stats
from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from statsmodels.stats.contingency_tables import mcnemar

rng = np.random.default_rng(22)          # কাল্পনিক বাজারঘর: অর্ডার ফেরত আসবে কি না
n = 1500
days = rng.gamma(3, 1.2, n)              # ডেলিভারিতে দিন
price = rng.lognormal(6.5, 0.6, n) / 100 # দাম, শত টাকা
first = rng.integers(0, 2, n)            # প্রথম অর্ডার?
z = -3 + 0.45 * days + 0.08 * price + 0.9 * first + 0.6 * (days > 5) * first
y = (rng.random(n) < 1 / (1 + np.exp(-z))).astype(int)
X = np.column_stack([days, price, first])
tr, te = slice(0, 1000), slice(1000, None)   # প্রথম ১০০০ train, শেষ ৫০০ test

pA = LogisticRegression(C=1.0).fit(X[tr], y[tr]).predict(X[te])
pB = DecisionTreeClassifier(max_depth=4, random_state=0).fit(X[tr], y[tr]).predict(X[te])
okA, okB = pA == y[te], pB == y[te]
print(f"accuracy  A (logistic) = {okA.mean():.3f}   B (tree) = {okB.mean():.3f}")

a = int(np.sum(okA & okB)); b = int(np.sum(okA & ~okB))
c = int(np.sum(~okA & okB)); d = int(np.sum(~okA & ~okB))
print("table [[both right, only A right], [only B right, both wrong]]:")
print(np.array([[a, b], [c, d]]))
chi2 = (abs(b - c) - 1) ** 2 / (b + c)                # Edwards-এর continuity correction
print(f"chi2 (corrected) = {chi2:.3f}, p = {stats.chi2.sf(chi2, 1):.4f}")
print(f"exact binomial: p = {stats.binomtest(b, b + c, 0.5).pvalue:.4f}")
r1 = mcnemar([[a, b], [c, d]], exact=True)
r2 = mcnemar([[a, b], [c, d]], exact=False, correction=True)
print(f"statsmodels exact p = {r1.pvalue:.4f}, chi2 p = {r2.pvalue:.4f}")

accuracy  A (logistic) = 0.678   B (tree) = 0.634
table [[both right, only A right], [only B right, both wrong]]:
[[291  48]
 [ 26 135]]
chi2 (corrected) = 5.959, p = 0.0146
exact binomial: p = 0.0141
statsmodels exact p = 0.0141, chi2 p = 0.0146

p-value-এর পাশে পার্থক্যের interval

p-value বলে “শূন্য পার্থক্যের সঙ্গে ডেটা কতটা বেমানান”, কিন্তু পার্থক্য কত বড় হতে পারে বলে না। জোড়া অনুপাতের পার্থক্য d̂ = (b − c)/n-এর standard error:

SE = √[ (b + c)/n − d̂² ] / √n

এখানে (b + c)/n = ০.১৪৮ আর d̂² = ০.০০১৯৪, তাই SE = √০.১৪৬১ / √৫০০ = ০.৩৮২২/২২.৩৬ ≈ ০.০১৭১. ৯৫% interval (Wald রূপ) ০.০৪৪ ± ১.৯৬ × ০.০১৭১ = [০.০১০৫, ০.০৭৭৫]. রিপোর্টের লাইন: “একই ৫০০টি test অর্ডারে A − B = +৪.৪ শতাংশ বিন্দু (৯৫% CI ১.১ থেকে ৭.৮), McNemar exact p = ০.০১৪”.

৪ · CV fold কেন স্বাধীন নয় — corrected test আর 5×2cv

এবার দ্বিতীয় প্রশ্ন: algorithm A কি এই আকারের ডেটায় B-র চেয়ে ভালো? এখানে training set-এর এলোমেলোতাও গুনতে হয়, তাই একটি test set যথেষ্ট নয়; চাই অনেক train/test ভাগ। সবচেয়ে প্রচলিত পথ: একই k-fold ভাগে দুজনকে চালিয়ে প্রতিটি fold-এর পার্থক্য d₁, …, d_k নেওয়া, তারপর paired t-test: t = d̄ / √(s²/k), k − 1 df. সমস্যা হলো t-test ধরে নেয় d-গুলো স্বাধীন। নিচের stepper দেখায় কেন তারা নয়।

Yoshua Bengio আর Yves Grandvalet ২০০৪ সালে দেখান যে k-fold CV-র variance-এর কোনো সাধারণ unbiased অনুমান নেই; তাই প্রতিটি সংশোধনই আনুমানিক। Claude Nadeau আর Bengio-র সংশোধন (২০০৩, Machine Learning জার্নাল) এসেছিল random subsampling-এর জন্য, আর Remco Bouckaert ও Eibe Frank ২০০৪ সালে একে ১০ × ১০ repeated CV-তে প্রয়োগ করে ভালো ফল পান। পর্ব ০৪-এ এই সংশোধন দিয়ে interval চওড়া করেছিলেন; এখানে আরেকটি পথ।

Dietterich-এর 5×2cv paired t-test

ভাবনাটি হলো fold-এর overlap সীমিত রাখা। ডেটা অর্ধেক অর্ধেক ভাগ, দুজনকে প্রথম অর্ধেকে train করে দ্বিতীয়তে মাপা, তারপর উল্টো। দুটি পার্থক্য p_i⁽¹⁾ আর p_i⁽²⁾, গড় p̄_i, আর variance s_i² = (p_i⁽¹⁾ − p̄_i)² + (p_i⁽²⁾ − p̄_i)². এটি ৫ বার, প্রতিবার নতুন এলোমেলো ভাগে। Statistic:

t̃ = p₁⁽¹⁾ / √( (1/5) ∑_i=1..5 s_i² ),   ৫ df-এর t-distribution

লবে শুধু প্রথম ভাগের প্রথম পার্থক্য, কারণ Dietterich-এর বিশ্লেষণে তখন লব আর হর প্রায় স্বাধীন থাকে, যা t-distribution-এর শর্ত। দুটি অর্ধেক training set-এ কোনো মিল নেই, আর ৫টি পুনরাবৃত্তি variance-এর অনুমান স্থির করে। দাম হলো প্রতিটি model শেখে মাত্র অর্ধেক ডেটায়, আর লবে একটিমাত্র সংখ্যা থাকায় power কম।

৫ · Null সত্য হলে কে কত ভুল করে — একটি ছোট simulation

দাবি যাচাইয়ের সৎ পথ হলো উত্তর-জানা জগৎ বানানো। নিচের ব্লকে দুটি feature x₁ আর x₂ সমান শক্তিশালী, আর label নির্ভর করে x₁ + x₂-এর ওপর। A শুধু x₁ দেখে, B শুধু x₂; গঠনগতভাবে দুজন হুবহু সমান ভালো, তাই null সত্য। ২০০টি নতুন ডেটাসেট (প্রতিটিতে ২০০ সারি) বানিয়ে চারটি test চালানো হয়, আর গোনা হয় কতবার p < ০.০৫ এসেছে। সৎ test-এর হার প্রায় ০.০৫ হওয়ার কথা; ২০০ বারে দৈবের কারণে ±০.০৩-এর মতো ওঠানামা স্বাভাবিক। Classifier হিসেবে খুব সরল একটি নিয়ম (দুই class-এর গড়ের মাঝে সীমা), যাতে হাজার হাজার fit ব্রাউজারে দ্রুত চলে।

import numpy as np
from scipy import stats

def fit_predict(xtr, ytr, xte):          # ছোট classifier: দুই class-এর গড়ের মাঝখানে সীমা
    m0, m1 = xtr[ytr == 0].mean(), xtr[ytr == 1].mean()
    return ((xte > (m0 + m1) / 2) == (m1 > m0)).astype(int)

def diff(X, y, te):                      # A দেখে শুধু x1, B শুধু x2: গঠনগতভাবে সমান ভালো
    tr = np.ones(len(y), bool); tr[te] = False
    return (np.mean(fit_predict(X[tr, 0], y[tr], X[te, 0]) == y[te])
            - np.mean(fit_predict(X[tr, 1], y[tr], X[te, 1]) == y[te]))

rng = np.random.default_rng(4)
n, reps, k = 200, 200, 10
names = ["naive t, 10-fold", "naive t, 10x10 CV", "corrected t, 10x10", "5x2cv t"]
rej = np.zeros(4)
for r in range(reps):                    # প্রতিবার নতুন ডেটাসেট, null সত্য
    X = rng.normal(size=(n, 2)); y = (X[:, 0] + X[:, 1] + rng.normal(0, 1.5, n) > 0).astype(int)
    d = np.array([[diff(X, y, f) for f in np.array_split(rng.permutation(n), k)]
                  for _ in range(10)])   # ১০ বার ১০-fold = ১০০টি পার্থক্য
    D, J = d.ravel(), d.size
    t1 = d[0].mean() / np.sqrt(d[0].var(ddof=1) / k)
    t2 = D.mean() / np.sqrt(D.var(ddof=1) / J)
    t3 = D.mean() / np.sqrt(D.var(ddof=1) * (1 / J + 1 / (k - 1)))   # Nadeau–Bengio
    p1s, s2 = [], []
    for i in range(5):                   # Dietterich-এর 5×2cv
        h = rng.permutation(n)[: n // 2]
        a, b = diff(X, y, h), diff(X, y, np.setdiff1d(np.arange(n), h))
        p1s.append(a); s2.append((a - b) ** 2 / 2)
    t4 = p1s[0] / np.sqrt(np.mean(s2))
    pv = [2 * stats.t.sf(abs(t), df) for t, df in ((t1, k - 1), (t2, J - 1), (t3, J - 1), (t4, 5))]
    rej += np.array(pv) < 0.05
print(f"{reps} datasets, true difference = 0, alpha = 0.05")
for nm, v in zip(names, rej / reps):
    print(f"  {nm:20s} rejects {v:.3f}")

200 datasets, true difference = 0, alpha = 0.05
  naive t, 10-fold     rejects 0.045
  naive t, 10x10 CV    rejects 0.490
  corrected t, 10x10   rejects 0.025
  5x2cv t              rejects 0.045

কিন্তু “একক ১০-fold-এ naive t ঠিক আছে” সিদ্ধান্তটি এই একটি স্থির classifier-এর। একই কোডে শুধু fit_predict-কে 1-nearest neighbour দিয়ে বদলানো হলো (return ytr[np.abs(xte[:, None] - xtr[None, :]).argmin(1)]), যা খুব অস্থির, কারণ প্রতিটি অনুমান একটি training সারির ওপর নির্ভর করে। ডেস্কটপে সত্যিকারের run (Python 3.14.4, NumPy 2.5.4, SciPy 1.18.1, একই seed, একটি run):

$ python null_sim_1nn.py
200 datasets, true difference = 0, alpha = 0.05
  naive t, 10-fold     rejects 0.115
  naive t, 10x10 CV    rejects 0.620
  corrected t, 10x10   rejects 0.110
  5x2cv t              rejects 0.125

অস্থির learner-এ একক ১০-fold-এর naive t ১১.৫%, repeated CV-র naive t ৬২%. আর লক্ষ করুন, এখানে সংশোধিত test আর 5×2cv-ও ১১–১৩%-এ উঠেছে। কারণটি গভীর: এখানে null হলো “অনেক ডেটাসেট জুড়ে গড়ে A আর B সমান”. কিন্তু একটি নির্দিষ্ট ২০০ সারির ডেটাসেটে, 1-NN-এর কাছে দৈবক্রমে x₁ বা x₂ একটু বেশি কাজের হতে পারে। একটি ডেটাসেটের ভেতরের কোনো resampling সেই ডেটাসেট-থেকে-ডেটাসেট ওঠানামা পুরো দেখতে পায় না, আর Bengio–Grandvalet-এর ফলাফল ঠিক এটিই বলে। শিক্ষা তিনটি: repeated CV-তে naive t কখনো নয়; corrected test বা 5×2cv ব্যবহার করুন, তবে এগুলোকে আনুমানিক ধরুন; আর সম্ভব হলে একাধিক স্বাধীন ডেটাসেট বা একটি বড় তালাবন্ধ test set জোগাড় করুন। (একটি simulation-এর সংখ্যা সাধারণ সত্য নয়; অন্য learner আর ডেটায় হার আলাদা হবে।)

৬ · অনেক ডেটাসেটে অনেক model — average rank আর Friedman test

বাজারঘর আটটি জেলায় আলাদা করে ফেরত-অনুমানের model রাখতে চায়, আর চারটি algorithm-এর মধ্যে বাছাই করতে চায়। প্রতিটি জেলার ডেটা আলাদা, কঠিনতাও আলাদা। প্রথম প্রলোভন: আট জেলার accuracy-র গড় নিয়ে তুলনা। Janez Demšar ২০০৬ সালে Journal of Machine Learning Research-এ (এই বিষয়ের বহুল উদ্ধৃত প্রবন্ধ) এর বিরুদ্ধে যুক্তি দেন: সহজ ডেটাসেটের ০.০১ আর কঠিনটির ০.০১ সমান অর্থ বহন করে না, আর একটি অদ্ভুত ডেটাসেট গড় টেনে নিতে পারে। তাঁর সুপারিশ: দুই algorithm-এর জন্য Wilcoxon signed-rank test, অনেকগুলোর জন্য Friedman test.

Friedman test (অর্থনীতিবিদ Milton Friedman, ১৯৩৭) প্রতিটি ডেটাসেটে k-টি algorithm-কে rank দেয় (১ = সেরা, সমান হলে গড় rank), তারপর প্রতিটি algorithm-এর average rank R_j নেয়। Null-এ সব algorithm সমান, তাই প্রত্যাশিত rank (k + 1)/2. N-টি ডেটাসেটে:

χ²_F = [12N / (k(k + 1))] · [ ∑_j R_j² − k(k + 1)²/4 ],   প্রায় χ²_k−1

Null বাতিল হলে প্রশ্ন আসে, কোন জোড়া আলাদা? Demšar-এর প্রস্তাবিত post-hoc হলো Nemenyi test: দুটি algorithm-এর average rank-এর পার্থক্য critical difference CD = q_α√(k(k + 1)/(6N))-এর বেশি হলে তারা আলাদা। q_α আসে studentized range বণ্টন থেকে, √2 দিয়ে ভাগ করে; k = ৪-এ q_0.05 ≈ ২.৫৬৯. এই CD-ই Demšar-এর পরিচিত “critical difference diagram”-এর ভিত্তি, যেখানে algorithm-গুলো average rank-এর রেখায় বসে আর CD-র মধ্যে থাকা দলগুলো একটি দাগে যুক্ত থাকে।

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.datasets import make_classification
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB

models = {"logistic": LogisticRegression(),
          "tree": DecisionTreeClassifier(max_depth=4, random_state=0),
          "knn": KNeighborsClassifier(15, algorithm="kd_tree"),
          "nbayes": GaussianNB()}
rows = {}
for s in range(8):                       # ৮টি কাল্পনিক জেলা, প্রতিটিতে আলাদা ফেরত-ডেটা
    X, y = make_classification(n_samples=300, n_features=6, n_informative=3 + s % 3,
                               n_redundant=1, flip_y=0.05 * (1 + s % 2), random_state=s)
    cv = StratifiedKFold(5, shuffle=True, random_state=s)
    rows[f"district {s + 1}"] = {m: cross_val_score(est, X, y, cv=cv).mean()
                                 for m, est in models.items()}
acc = pd.DataFrame(rows).T.round(3)
ranks = acc.rank(axis=1, ascending=False)        # প্রতিটি সারিতে ১ = সেরা
print(acc.assign(best=acc.idxmax(axis=1)))
print("average rank:", ranks.mean().round(3).to_dict())
stat, p = stats.friedmanchisquare(*[acc[m] for m in acc.columns])
print(f"Friedman chi2 = {stat:.3f}, p = {p:.4f}")
k, N = acc.shape[1], acc.shape[0]
q = stats.studentized_range.ppf(0.95, k, 1e6) / np.sqrt(2)   # Nemenyi-র q_0.05
cd = q * np.sqrt(k * (k + 1) / (6 * N))
print(f"q_0.05 = {q:.3f}, critical difference CD = {cd:.3f}")

            logistic   tree    knn  nbayes      best
district 1     0.697  0.847  0.870   0.760       knn
district 2     0.870  0.807  0.870   0.837  logistic
district 3     0.867  0.840  0.887   0.880       knn
district 4     0.843  0.800  0.843   0.830  logistic
district 5     0.823  0.873  0.907   0.823       knn
district 6     0.837  0.827  0.910   0.810       knn
district 7     0.867  0.833  0.893   0.827       knn
district 8     0.817  0.757  0.837   0.767       knn
average rank: {'logistic': 2.438, 'tree': 3.25, 'knn': 1.125, 'nbayes': 3.188}
Friedman chi2 = 14.610, p = 0.0022
q_0.05 = 2.569, critical difference CD = 1.658

৭ · একাধিক তুলনা — FWER, FDR, Bonferroni, Holm আর BH

অংশ ৬-এর ছয়টি জোড়া, বা ২০টি জেলায় নতুন বনাম পুরোনো model: প্রতিটিতে ৫%-এ test চালালে, null সব জায়গায় সত্য হলেও অন্তত একটি মিথ্যা আবিষ্কারের সম্ভাবনা দ্রুত বাড়ে। m-টি স্বাধীন test-এ তা ১ − ০.৯৫^m; m = ২০-এ প্রায় ০.৬৪. পর্ব ৩৭-এ এটি simulation-এ দেখেছেন, আর পর্ব ০৪-এ “garden of forking paths” হিসেবে। এখানে প্রতিকারগুলো নিখুঁতভাবে। দুটি আলাদা লক্ষ্য আছে:

Bonferroni: প্রতিটি test-কে α/m-এ চালানো। যুক্তি পর্ব ০৩-এর union bound: P(কোনো একটি মিথ্যা) ≤ ∑ P(প্রতিটি মিথ্যা) ≤ m · α/m = α. কোনো স্বাধীনতার শর্ত লাগে না, তবে m বড় হলে খুব কড়া। Holm (Sture Holm, ১৯৭৯, Scandinavian Journal of Statistics): p-value-গুলো ছোট থেকে বড় সাজিয়ে p₍₁₎ ≤ … ≤ p_(m); i-তম-কে তুলনা করা α/(m − i + 1)-এর সঙ্গে, প্রথম যেখানে ব্যর্থ সেখানে থামা, আর তার আগের সবগুলো বাতিল। এটিও যেকোনো নির্ভরশীলতায় FWER ≤ α রাখে, আর সবসময় Bonferroni-র সমান বা বেশি আবিষ্কার করে; তাই Bonferroni-র বদলে Holm না নেওয়ার কারণ প্রায় নেই। Benjamini–Hochberg (Yoav Benjamini আর Yosef Hochberg, ১৯৯৫): সবচেয়ে বড় i খোঁজা যেখানে p_(i) ≤ i·q/m, আর ১ থেকে i পর্যন্ত সব বাতিল। স্বাধীন (বা এক ধরনের ধনাত্মক নির্ভরশীল) test-এ FDR ≤ (m₀/m)·q, যেখানে m₀ হলো সত্য null-এর সংখ্যা। যেকোনো নির্ভরশীলতার জন্য Benjamini আর Daniel Yekutieli ২০০১ সালে একটি রক্ষণশীল রূপ দেন।

৮ · ২০টি p-value আর জানা উত্তর — কে কত মিথ্যা আবিষ্কার করে

বাজারঘর ২০টি জেলায় নতুন সুপারিশ-model আর পুরোনোটি পাশাপাশি চালিয়েছে; প্রতিটি জেলায় ৪০টি জোড়া পার্থক্য (ধরুন একই দিনে দুই model-এর ঝুড়ির মূল্যের পার্থক্য, standardize করা)। Simulation-এ আমরা জানি: শুধু প্রথম ৫টি জেলায় আসল উন্নতি আছে (গড় ০.৬ SD), বাকি ১৫টিতে কিছুই নেই। প্রতিটি জেলায় paired t-test, তারপর চারটি নীতি। প্রথমে একটি পরীক্ষা খুঁটিয়ে, তারপর ২,০০০ বার পুনরাবৃত্তি করে দীর্ঘমেয়াদি হার। Holm আর BH নিজে vectorize করে লেখা, কারণ multipletests একবারে একটি সারি নেয় আর হাজার বার ডাকলে ব্রাউজারে ধীর; প্রথম অংশে statsmodels-এর সঙ্গে মিলিয়ে নেওয়া হয়েছে।

import numpy as np
from scipy import stats
from statsmodels.stats.multitest import multipletests

m, a = 20, 0.05                          # ২০টি জেলায় নতুন বনাম পুরোনো model; শুধু প্রথম ৫টিতে আসল উন্নতি
true = np.zeros(m, bool); true[:5] = True
def pvalues(rng, R, n=40):               # প্রতিটি জেলায় ৪০টি জোড়া পার্থক্য → paired t-test
    d = rng.normal(0.6 * true[:, None], 1.0, size=(R, m, n))
    return stats.ttest_1samp(d, 0.0, axis=2).pvalue           # আকার (R, m)

def decide(P, meth):                     # P-এর প্রতিটি সারি একটি পুরো পরীক্ষা
    if meth == "none": return P <= a
    if meth == "bonferroni": return P <= a / m
    i = np.arange(1, m + 1); ps = np.sort(P, axis=1)
    ok = ps <= (a / (m - i + 1) if meth == "holm" else a * i / m)
    k = ok.cumprod(1).sum(1) if meth == "holm" else np.where(ok, i, 0).max(1)
    return P.argsort(1).argsort(1) < k[:, None]              # সবচেয়ে ছোট k-টি বাতিল

methods = ["none", "bonferroni", "holm", "fdr_bh"]
P1 = pvalues(np.random.default_rng(22), 1)
print("smallest 7 p:", " ".join(f"{v:.4f}" for v in np.sort(P1[0])[:7]))
for meth in methods:
    rej = decide(P1, meth)[0]
    same = meth == "none" or (rej == multipletests(P1[0], a, method=meth)[0]).all()
    print(f"{meth:10s} found {rej.sum()}, false {(rej & ~true).sum()}, statsmodels agrees: {same}")

P = pvalues(np.random.default_rng(0), 2000)          # ২০০০ বার পুরো পরীক্ষা
print("2000 repeats:  P(any false)  mean FDP  power")
for meth in methods:
    rej = decide(P, meth); V = (rej & ~true).sum(1)
    fdp = V / np.maximum(rej.sum(1), 1)
    print(f"  {meth:10s} {np.mean(V > 0):9.3f} {fdp.mean():9.3f} {rej[:, true].mean():6.3f}")

smallest 7 p: 0.0000 0.0001 0.0001 0.0037 0.0045 0.0443 0.0714
none       found 6, false 1, statsmodels agrees: True
bonferroni found 3, false 0, statsmodels agrees: True
holm       found 3, false 0, statsmodels agrees: True
fdr_bh     found 5, false 0, statsmodels agrees: True
2000 repeats:  P(any false)  mean FDP  power
  none           0.542     0.120  0.959
  bonferroni     0.039     0.009  0.707
  holm           0.048     0.011  0.728
  fdr_bh         0.190     0.038  0.868

ML-এ লুকানো বহু তুলনাও আছে: ২০০টি config-এর সেরাটিকে একই validation set-এ baseline-এর সঙ্গে তুলনা, ১৫টি metric-এর যেটিতে উন্নতি সেটিই রিপোর্ট, বা একই test set বারবার দেখা। সংখ্যা জানা থাকলে সংশোধন করুন; তার চেয়েও নির্ভরযোগ্য হলো বাছাই আর মূল্যায়ন আলাদা ডেটায় রাখা, যেমন তালাবন্ধ test set বা nested CV.

৯ · Statistical বনাম practical significance, Bayesian দৃষ্টি আর কোন test কখন

p < ০.০৫ মানে শুধু “শূন্য পার্থক্যের সঙ্গে ডেটা বেমানান”. ১০ লাখ test সারিতে ০.০২ শতাংশ বিন্দুর পার্থক্যও significant হবে, অথচ তা নতুন model চালানোর খরচ, latency বা রক্ষণাবেক্ষণের দাম পোষায় না। উল্টো দিকে ২০০ সারিতে ৪ শতাংশ বিন্দুর আসল উন্নতি p > ০.০৫ দিতে পারে, কারণ power কম। American Statistical Association ২০১৬ সালে p-value নিয়ে বিবৃতিতে (Wasserstein আর Lazar) ঠিক এটিই জোর দিয়ে বলে: p-value প্রভাবের আকার বা গুরুত্ব মাপে না। তাই রিপোর্টে সবসময় effect size আর তার interval দিন, আর আগে থেকে ঠিক করুন কত বড় পার্থক্য ব্যবসার কাছে অর্থপূর্ণ।

ব্যবহারিক সীমার ধারণাটি Bayesian তুলনার কেন্দ্রে। Benavoli, Corani, Demšar আর Zaffalon ২০১৭ সালে JMLR-এ “Time for a change” শিরোনামের একটি টিউটোরিয়ালে প্রস্তাব করেন: একটি ROPE (region of practical equivalence) ঠিক করুন, যেমন accuracy-র পার্থক্য ±০.০১-এর মধ্যে হলে “কার্যত সমান”. তারপর পার্থক্যের posterior থেকে তিনটি সম্ভাবনা বের করুন: P(A ভালো), P(কার্যত সমান), P(B ভালো). উত্তর তখন “p = ০.০৩” নয়, বরং “৮৫% সম্ভাবনায় A অন্তত ১ শতাংশ বিন্দু ভালো”, যা সিদ্ধান্তের ভাষার কাছাকাছি। CV-র correlation সামলাতে Corani আর Benavoli ২০১৫ সালে একটি Bayesian correlated t-test দেন, যেখানে Nadeau–Bengio-র ρ = n_test/n একই ভূমিকা নেয়।

১০ · AI-এর যুগে — “significant উন্নতি” দাবি যাচাই

AI assistant-এর দেওয়া তুলনার কোড প্রায় সবসময় চলে। বিপদ হলো ভুল test-ও একটি p-value দেয়, যা দেখতে সঠিকটির মতোই। পরিচিত ভুলগুলো:

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

• প্রশ্নটি লিখে নিয়েছি: নির্দিষ্ট model, নাকি algorithm, নাকি অনেক ডেটাসেট
• দুই model একই সারিতে মাপা হয়েছে, আর test paired (McNemar বা fold-ভিত্তিক)
• Repeated CV-তে naive t নেই; corrected t বা 5×2cv ব্যবহার হয়েছে
• কতগুলো তুলনা হয়েছে গুনেছি, আর Holm বা BH দিয়েছি
• Test set তালাবন্ধ ছিল; tuning বা model বাছাই তাতে হয়নি
• একটি সরল baseline-এর সঙ্গে কয়েকটি seed বা fold-এ তুলনা আছে
• Effect size আর interval আছে, আর তা ব্যবহারিক সীমার সঙ্গে মিলিয়েছি
• Function-এর নাম docs-এ মিলিয়েছি, আর কোনো আসল গ্রাহকের ডেটা AI-তে পেস্ট করিনি