Model বোঝা — permutation importance, PDP আর ICE

এতদিন প্রশ্ন ছিল “model কতটা ভালো অনুমান করে?” এই মডিউলের প্রশ্ন আলাদা: model কীসের ওপর ভর করে, তার অনুমান কতটা নিশ্চিত, আর কার জন্য ন্যায্য। প্রথম পর্বে আমরা পুরো model-কে একসঙ্গে বুঝব। ভিত্তি তিনটি: AI-এর ভিত্তি কোর্সের পর্ব ৬৩-এ permutation importance-এর প্রথম পরিচয়, পর্ব ৩৩-এ impurity বনাম permutation importance আর high-cardinality পক্ষপাত (সেই demo এখানে আর দেখানো হবে না), আর পর্ব ৩৬-এ partial_dependence-এর প্রথম ব্যবহার। চারটি চালানো যায় এমন ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, pandas, scikit-learn আর matplotlib); প্রথমবার scikit-learn নামতে প্রায় দশ সেকেন্ড লাগে।

১ · Model বোঝার প্রশ্ন — global বনাম local, শুরু থেকে বনাম পরে

Interpretability মানে মানুষ model-এর আচরণ কতটা বুঝতে পারে। দুটি অক্ষে ভাবা কাজের। প্রথম অক্ষ: global ব্যাখ্যা (পুরো ডেটায় গড় আচরণ) বনাম local ব্যাখ্যা (একটি সারি)। দ্বিতীয় অক্ষ: model নিজেই পড়া যায় এমন, অর্থাৎ শুরু থেকেই ব্যাখ্যাযোগ্য (ছোট linear model, ছোট গাছ, GAM), নাকি একটি black box-এর ওপর পরে জোড়া post-hoc ব্যাখ্যা. AI-এর ভিত্তি কোর্সের পর্ব ৭২ explainability-কে দায়িত্বশীল AI-এর অংশ হিসেবে পরিচয় করিয়েছে; এখানে যন্ত্রগুলো গভীরে।

সব post-hoc পদ্ধতির একটি সাধারণ সীমা আছে: এরা model কী করে তা বর্ণনা করে, দুনিয়া কীভাবে চলে তা নয়। Model যদি ভুল কারণে ঠিক উত্তর দেয় (যেমন leakage-এর কলাম ধরে), importance সেই ভুল কারণকেই “গুরুত্বপূর্ণ” দেখাবে। Cynthia Rudin ২০১৯ সালে Nature Machine Intelligence-এ যুক্তি দেন, উচ্চ ঝুঁকির সিদ্ধান্তে black box ব্যাখ্যা করার চেয়ে শুরু থেকেই ব্যাখ্যাযোগ্য model বানানো ভালো, কারণ tabular ডেটায় দুই ধরনের model-এর accuracy প্রায়ই কাছাকাছি। সবাই এই অবস্থানের সঙ্গে পুরো একমত নন (ছবি, লেখা বা খুব জটিল সম্পর্কে ফাঁক বড় হতে পারে), তবে অংশ ৮-এ দেখবেন, এই পর্বের ডেটায় যুক্তিটি সংখ্যায় খাটে।

২ · বাজারঘরের ডেলিভারি-সময়ের ডেটা আর একটি black box

বাজারঘর জানতে চায় একটি অর্ডার কত মিনিটে পৌঁছাবে। নিচের module ৩,০০০টি কাল্পনিক ডেলিভারি বানায়। ছয়টি feature: দোকান থেকে বাড়ির সোজা দূরত্ব (distance_km), map app-এর আন্দাজি রাস্তার দূরত্ব (road_km), বৃষ্টি (rain), পণ্যের সংখ্যা (items), অর্ডারের ঘণ্টা (hour) আর কুপন (coupon). আসল নিয়ম simulation-এ জানা, তাই প্রতিটি ব্যাখ্যা যাচাই করা যাবে: সময় নির্ভর করে আসল পথের দৈর্ঘ্যের ওপর, বৃষ্টিতে প্রতি km ৭০% বেশি সময় (একটি interaction), প্রতিটি পণ্যে ১ মিনিট, বিকেল ৫টা থেকে রাত ৯টার আগ পর্যন্ত ভিড়ে ৭ মিনিট। কুপনের কোনো প্রভাব নেই। আসল পথ ডেটায় নেই; সোজা দূরত্ব আর app-এর আন্দাজ দুটিই তার অসম্পূর্ণ ছায়া, তাই দুটি খুব ঘনিষ্ঠ (r = ০.৯৮১)। রাস্তার দূরত্ব কখনো সোজা দূরত্বের চেয়ে কম হয় না; আসল ডেটায় অনুপাত ১.০৩ থেকে ১.৬-এর মধ্যে। এই তথ্য অংশ ৭-এ কাজে লাগবে।

# বাজারঘর (কাল্পনিক দোকান): ডেলিভারি কত মিনিটে, seeded ডেটা
import numpy as np
import pandas as pd

COLS = ["distance_km", "road_km", "rain", "items", "hour", "coupon"]

def deliveries(n=3000, seed=56):
    rng = np.random.default_rng(seed)
    distance = rng.uniform(0.5, 10.0, n)          # সোজা দূরত্ব
    route = distance * rng.uniform(1.15, 1.45, n) # আসল পথ (ডেটায় নেই)
    road = route * rng.uniform(0.9, 1.1, n)       # app-এর আন্দাজি রাস্তার দূরত্ব
    rain = rng.binomial(1, 0.3, n)
    items = rng.integers(1, 9, n)
    hour = rng.integers(9, 23, n)
    coupon = rng.binomial(1, 0.4, n)              # সময়ে কোনো প্রভাব নেই
    rush = ((hour >= 17) & (hour <= 20)).astype(int)
    minutes = (12 + 2.4 * route * (1 + 0.7 * rain) + 1.0 * items
               + 7 * rush + rng.normal(0, 3, n))  # আসল নিয়ম: পথ × বৃষ্টি
    X = pd.DataFrame({"distance_km": distance.round(2), "road_km": road.round(2),
                      "rain": rain, "items": items, "hour": hour, "coupon": coupon})
    return X, minutes.round(1)

def split():
    X, y = deliveries()
    return X[:2000], X[2000:], y[:2000], y[2000:]

Black box হিসেবে নেওয়া হলো পর্ব ৩৬-এর HistGradientBoostingRegressor, ২০০টি গাছ। প্রথম ২,০০০ সারিতে শেখানো, শেষ ১,০০০ সারি held-out test. সব ব্লকে একই ভাগ, একই seed.

৩ · Permutation importance — held-out ডেটায়, বারবার, ছড়ানোসহ

ধারণাটি Leo Breiman ২০০১ সালে random forest-এর প্রবন্ধে দেন, আর পরে এটি যেকোনো model-এ খাটানো হয়। শেখানো model স্থির রাখুন। Held-out ডেটায় একটি কলামের মানগুলো সারিগুলোর মধ্যে এলোমেলো করুন: কলামের বণ্টন একই থাকে, কিন্তু label-এর সঙ্গে সম্পর্ক ভেঙে যায়। এবার score কতটা খারাপ হলো, সেটিই permutation importance:

PI_j = (১/R) ∑_r=১^R [ MAE(এলোমেলো j, বার r) − MAE(আসল) ]

এখানে j হলো কলাম, R হলো কতবার এলোমেলো করা হলো (n_repeats). প্রতিটি এলোমেলো করা একটি দৈব ঘটনা, তাই R বার করে গড় আর standard deviation দুটিই দেখতে হয়। MAE ব্যবহার করায় importance-এর একক মিনিট: “এই কলামের তথ্য হারালে গড় ভুল কত মিনিট বাড়ে”। নিচের stepper পাঁচটি সারিতে পুরো হিসাব করে।

তিনটি সিদ্ধান্ত ফলকে বদলায়। (১) কোন ডেটায়: train ডেটায় মাপলে model যা মুখস্থ করেছে তাও “গুরুত্ব” পায়; held-out ডেটায় মাপলে শুধু যা নতুন ডেটায় কাজে লাগে। (২) কোন score: R², MAE বা log loss আলাদা ক্রম দিতে পারে; যে মাপে model-কে বিচার করেন, সেটিই নিন। (৩) কতবার: ছোট test সেটে spread বড়। প্রথম ব্লক নিজে লেখা কোড sklearn-এর সঙ্গে মেলায়, train বনাম test তুলনা করে, আর পরের অংশের জন্য correlated জোড়ার grouped আর drop-column সংখ্যাও ছাপে।

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor as HGB
from sklearn.inspection import permutation_importance
from sklearn.metrics import mean_absolute_error as mae
from bazar_eta import split, COLS

Xtr, Xte, ytr, yte = split()
m = HGB(max_iter=200, random_state=0).fit(Xtr, ytr)
base = mae(yte, m.predict(Xte))
print(f"test MAE {base:.2f} min (train {mae(ytr, m.predict(Xtr)):.2f})")

def perm_imp(cols, n_repeats=10, seed=0):    # নিজে লেখা, sklearn-এর বীজের নিয়মে
    s = np.random.RandomState(seed).randint(2**31)
    rs, idx, Xp = np.random.RandomState(s), np.arange(len(Xte)), Xte.copy()
    d = []
    for _ in range(n_repeats):
        rs.shuffle(idx)
        Xp[cols] = Xp[cols].to_numpy()[idx]  # এক বা একাধিক কলাম একসঙ্গে এলোমেলো
        d.append(mae(yte, m.predict(Xp)) - base)
    return np.array(d)

mine = np.array([perm_imp([c]) for c in COLS])
kw = dict(scoring="neg_mean_absolute_error", n_repeats=10, random_state=0)
sk = permutation_importance(m, Xte, yte, **kw)
print("max |mine - sklearn| =", np.abs(mine - sk.importances).max())
tr = permutation_importance(m, Xtr, ytr, **kw).importances_mean
print(pd.DataFrame({"test mean": mine.mean(1), "test sd": mine.std(1),
                    "train mean": tr}, index=COLS).round(3))
g = perm_imp(["distance_km", "road_km"])
print(f"grouped distance+road: {g.mean():.3f} +/- {g.std():.3f}")
for cols in [["distance_km"], ["road_km"], ["distance_km", "road_km"]]:
    r = HGB(max_iter=200, random_state=0).fit(Xtr.drop(columns=cols), ytr)
    print("drop-column", cols, f"{mae(yte, r.predict(Xte.drop(columns=cols))) - base:.3f}")

test MAE 2.74 min (train 1.64)
max |mine - sklearn| = 0.0
             test mean  test sd  train mean
distance_km      3.330    0.109       3.869
road_km          4.606    0.156       5.243
rain             3.997    0.227       4.162
items            0.955    0.049       1.586
hour             1.679    0.073       2.321
coupon          -0.007    0.010       0.102
grouped distance+road: 9.616 +/- 0.248
drop-column ['distance_km'] 0.109
drop-column ['road_km'] 0.209
drop-column ['distance_km', 'road_km'] 6.900

৪ · Correlated feature — importance ভাগ হয়, grouped permutation আর drop-column

প্রথম ব্লকের শেষ পাঁচ লাইন একটি ধাঁধা দেখায়। দুটি দূরত্ব-কলাম আলাদাভাবে ৩.৩৩ আর ৪.৬১, যোগ ৭.৯৪. কিন্তু দুটি একই এলোমেলো ক্রমে একসঙ্গে সরালে ৯.৬২, যোগের চেয়েও বেশি। আর একটিকে পুরো বাদ দিয়ে model আবার শেখালে (drop-column importance) ক্ষতি মাত্র ০.১০৯ আর ০.২০৯ মিনিট; দুটি বাদ দিলে ৬.৯০০.

তিনটি পদ্ধতি তিনটি আলাদা প্রশ্নের উত্তর দেয়। একক permutation: “এই শেখানো model এই কলামের ওপর কতটা ভর করে?” Model দুটি কলামেই কিছু ভর দিয়েছে, তাই দুটিই মাঝারি, কিন্তু এলোমেলো করার সময় তৈরি হয় এমন সারি যা বাস্তবে নেই (১ km সোজা দূরত্ব, ১২ km রাস্তা)। Giles Hooker, Lucas Mentch ও Siyu Zhou ২০২১ সালে দেখান, এমন extrapolation importance-কে ফুলিয়ে বা বিকৃত করতে পারে। Grouped permutation: “এই তথ্যের দল হারালে কত ক্ষতি?” দুটি কলাম একই ক্রমে সরালে তাদের পারস্পরিক সম্পর্ক অটুট থাকে, অসম্ভব জোড়া তৈরি হয় না। ঘনিষ্ঠ কলামগুলোকে দলে বাঁধার একটি সরল উপায় correlation-এর ভিত্তিতে hierarchical clustering. Drop-column: “এই কলাম ছাড়া নতুন model কতটা খারাপ?” এটি ডেটায় কলামটির অনন্য তথ্যের মাপ, শেখানো model-এর নয়। দাম চড়া: প্রতিটি কলামে (বা দলে) একবার পুরো training, আর ১০০টি feature-এ ১০০ বার; সৎ হতে হলে প্রতিবার একই hyperparameter আর cross-validation.

৫ · Partial dependence আর ICE — feature বদলালে অনুমান কোন দিকে

Importance বলে “কতটা”, কিন্তু দিক বা আকার বলে না। Jerome Friedman ২০০১ সালে gradient boosting-এর প্রবন্ধে partial dependence দেন। একটি feature S-কে একটি মান v-তে বসান, বাকি সব কলাম (C) প্রতিটি সারিতে যেমন আছে তেমন রেখে, তারপর অনুমানের গড়:

PD_S(v) = (১/n) ∑_i=১ⁿ f̂(v, x_i,C)

অর্থাৎ “সব অর্ডারের রাস্তা যদি v km হতো, গড়ে কত মিনিট লাগত?” v বদলে আঁকা রেখাই PDP. গড় নেওয়ার আগের প্রতিটি সারির নিজের রেখাকে Alex Goldstein ও সহকর্মীরা ২০১৫ সালে নাম দেন ICE (individual conditional expectation): সারি i-এর জন্য v ↦ f̂(v, x_i,C). PDP হলো ICE রেখাগুলোর গড়।

ICE-এর রেখাগুলো প্রায়ই আলাদা উচ্চতায় শুরু হয় (একটি অর্ডারে ৮টি পণ্য, আরেকটিতে ১টি), তাই আকারের পার্থক্য চোখে পড়ে না। সমাধান centred ICE: প্রতিটি রেখা থেকে তার প্রথম বিন্দুর মান বাদ দেওয়া, যাতে সব রেখা শূন্য থেকে শুরু হয়। তখন রেখাগুলো ছড়িয়ে গেলে বোঝা যায়, feature-টির প্রভাব সারিভেদে আলাদা, অর্থাৎ অন্য কোনো feature-এর সঙ্গে interaction আছে। দ্বিতীয় ব্লক ৩০০টি test সারিতে ICE আর PDP নিজে বানায়, partial_dependence(method="brute", kind="both")-এর সঙ্গে মেলায়, আর বৃষ্টি অনুযায়ী রং করে দুটি ছবি আঁকে।

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import HistGradientBoostingRegressor as HGB
from sklearn.inspection import partial_dependence
from bazar_eta import split

Xtr, Xte, ytr, yte = split()
m = HGB(max_iter=200, random_state=0).fit(Xtr, ytr)
S = Xte[:300]                                  # ৩০০টি test সারি
grid = np.arange(1.0, 15.0)                    # road_km = 1, 2, ..., 14
ice = np.empty((len(S), len(grid)))            # প্রতিটি সারির নিজের রেখা
for j, v in enumerate(grid):
    Z = S.copy()
    Z["road_km"] = v                           # সব সারিতে একই মান বসানো
    ice[:, j] = m.predict(Z)
pdp = ice.mean(axis=0)                         # PDP = ICE রেখাগুলোর গড়
sk = partial_dependence(m, S, ["road_km"], custom_values={"road_km": grid},
                        method="brute", kind="both")
print("max |PDP - sklearn| =", np.abs(pdp - sk["average"][0]).max())
print("max |ICE - sklearn| =", np.abs(ice - sk["individual"][0]).max())
print("PDP at 2, 6, 10, 14 km:", pdp[[1, 5, 9, 13]].round(1))
wet = S.rain.to_numpy() == 1
slope = (ice[:, 11] - ice[:, 1]) / 10          # 2 থেকে 12 km, মিনিট প্রতি km
print(f"ICE slope min/km  dry {slope[~wet].mean():.2f}  rain {slope[wet].mean():.2f}"
      f"  PDP {(pdp[11] - pdp[1]) / 10:.2f}")
cice = ice - ice[:, [0]]                       # centred ICE: 1 km-এ সবাই শূন্য
fig, ax = plt.subplots(1, 2, figsize=(8, 3.4), sharex=True)
for a, Y, t in [(ax[0], ice, "ICE and PDP"), (ax[1], cice, "centred ICE")]:
    a.plot(grid, Y[~wet].T, color="tab:orange", alpha=0.15, lw=0.8)
    a.plot(grid, Y[wet].T, color="tab:blue", alpha=0.15, lw=0.8)
    a.set(title=t, xlabel="road_km (orange = dry, blue = rain)")
ax[0].plot(grid, pdp, color="black", lw=3, label="PDP")
ax[0].set_ylabel("predicted minutes"); ax[0].legend()
fig.subplots_adjust(wspace=0.25, bottom=0.16)

max |PDP - sklearn| = 4.263256414560601e-14
max |ICE - sklearn| = 0.0
PDP at 2, 6, 10, 14 km: [30.4 38.1 43.3 47.1]
ICE slope min/km  dry 1.33  rain 1.98  PDP 1.54

৬ · Interaction মাপা — 2-D PDP আর Friedman-এর H-statistic

ICE-এর ছড়ানো বলে “interaction আছে”, কিন্তু কার সঙ্গে, তা বলে না। একটি পথ দুটি feature একসঙ্গে বদলানো: 2-D partial dependence PD_jk(v_j, v_k), যা PartialDependenceDisplay.from_estimator(m, X, [("road_km", "rain")]) একটি heatmap বা দুই রেখা হিসেবে আঁকে। Interaction না থাকলে 2-D তলটি দুটি 1-D রেখার যোগফল: PD_jk = PD_j + PD_k (তিনটিই গড় শূন্যে কেন্দ্র করে)। স্টেপারের খেলনা model-এ রাস্তা ২ থেকে ৬ km হলে শুকনোয় ৮ মিনিট বাড়ে, বৃষ্টিতে ১৪; যোগফলের নিয়ম মানলে দুটিই সমান হতো।

Jerome Friedman আর Bogdan Popescu ২০০৮ সালে এই ফাঁককে একটি সংখ্যায় আনেন, H-statistic:

H^২_jk = ∑_i [PD_jk(x_ij, x_ik) − PD_j(x_ij) − PD_k(x_ik)]^২ / ∑_i PD_jk(x_ij, x_ik)^২

অর্থাৎ যৌথ প্রভাবের variance-এর কত ভাগ যোগফল দিয়ে ব্যাখ্যা হয় না। H = ০ মানে কোনো interaction নেই, ১-এর কাছে মানে প্রভাবের প্রায় সবটাই যৌথ। সীমা তিনটি। হিসাব দামি: প্রতিটি সারির বিন্দুতে PD_jk লাগে, তাই n সারিতে প্রায় n^২টি অনুমান; চর্চায় কয়েকশো সারির নমুনা নেওয়া হয়। সংখ্যাটি নমুনাভেদে ওঠানামা করে, আর কোনো সহজ “তাৎপর্যের” সীমা নেই। আর এটিও PD-র ওপর দাঁড়ানো, তাই correlated feature-এ PDP-র সব দুর্বলতা পায়। লেখার সময়ে scikit-learn-এ H-এর তৈরি ফাংশন নেই; R-এর কিছু package-এ আর কয়েকটি Python library-তে আছে, তবে নিজে লেখাও কঠিন নয়। পর্ব ৩৬-এর interaction_cst উল্টো দিক থেকে একই প্রশ্নে হাত দেয়: কোন feature-গুলো একসঙ্গে কাজ করতে পারবে, তা শেখার আগেই বেঁধে দেওয়া।

৭ · PDP যখন ভুল বলে — অসম্ভব সারি আর ALE

PDP-র সূত্রে একটি লুকানো অনুমান আছে: S-কে v-তে বসিয়ে বাকি কলাম অপরিবর্তিত রাখা যায়, অর্থাৎ S আর C স্বাধীন। road_km আর distance_km-এর মতো ঘনিষ্ঠ জোড়ায় এটি ভাঙে। নিচের ছবিতে ধূসর কীলকটি আসল ডেটার জায়গা: রাস্তা সোজা দূরত্বের ১.০৩ থেকে ১.৬ গুণ।

Daniel Apley আর Jingyu Zhu-র প্রবন্ধে (২০২০ সালে Journal of the Royal Statistical Society, Series B-তে প্রকাশিত; আগে preprint হিসেবে ছড়িয়েছিল) accumulated local effects (ALE) এর সমাধান। ধাপগুলো: feature-এর পরিসরকে quantile দিয়ে K খোপে ভাগ; প্রতিটি খোপে শুধু সেই খোপের সারি নিয়ে, তাদের feature-কে খোপের নিচের কিনারা z_k−1 আর উপরের কিনারা z_k-এ বসিয়ে অনুমানের পার্থক্যের গড়; তারপর খোপ থেকে খোপে যোগ করে চলা, আর শেষে গড় শূন্যে কেন্দ্র করা:

ALE(z_k) = ∑_m=১^k (১/n_m) ∑_{i ∈ খোপ m} [ f̂(z_m, x_i,C) − f̂(z_m−১, x_i,C) ] − ধ্রুবক

এখানে n_m হলো খোপ m-এর সারির সংখ্যা। PDP প্রশ্ন করে “সবাইকে v-তে নিলে?”; ALE প্রশ্ন করে “যারা আসলেই v-এর কাছে আছে, তাদের একটু বদলালে?” নিচের stepper একটি খেলনা model-এ পার্থক্যটি সংখ্যায় দেখায়।

খেলনা model-টি ইচ্ছা করে চরম। বাস্তবে কি এমন হয়? তৃতীয় ব্লক একই বাজারঘর ডেটায় দ্বিতীয় একটি black box শেখায়: standardise করে degree-4 polynomial feature, তারপর পর্ব ১৭-এর ridge. Test MAE প্রায় HGB-এর সমান। তারপর road_km-এর ALE নিজে লেখে (১০টি quantile খোপ), একই বিন্দুতে PDP-র সঙ্গে তুলনা করে, আর গুনে দেখে PDP-র প্রান্তের বিন্দুতে কত ভাগ সারি অসম্ভব।

import numpy as np
import matplotlib.pyplot as plt
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_absolute_error as mae
from bazar_eta import split
Xtr, Xte, ytr, yte = split()
m = make_pipeline(StandardScaler(), PolynomialFeatures(4), Ridge(1.0)).fit(Xtr, ytr)
print(f"degree-4 polynomial ridge: test MAE {mae(yte, m.predict(Xte)):.2f} min")
f, K = "road_km", 10
z = np.quantile(Xte[f], np.linspace(0, 1, K + 1))     # quantile সীমানা
b = np.clip(np.searchsorted(z, Xte[f], side="left") - 1, 0, K - 1)
lo, hi = Xte.copy(), Xte.copy()
lo[f], hi[f] = z[b], z[b + 1]                         # শুধু নিজের খোপের দুই কিনারা
d = m.predict(hi) - m.predict(lo)                     # স্থানীয় পরিবর্তন
n = np.bincount(b, minlength=K)
acc = np.r_[0, np.cumsum(np.bincount(b, weights=d, minlength=K) / n)]
ale = acc - np.sum((acc[:-1] + acc[1:]) / 2 * n) / n.sum()   # গড় শূন্যে কেন্দ্র
pdp = np.array([m.predict(Xte.assign(**{f: v})).mean() for v in z])
pdp -= pdp.mean()
print("grid :", z.round(1))
print("ALE  :", ale.round(1))
print("PDP  :", pdp.round(1))
print("corr(distance, road) =", round(Xte.distance_km.corr(Xte.road_km), 3))
for v in [z[0], z[-1]]:                    # আসল ডেটায় road/distance 1.03 থেকে 1.6
    r = v / Xte.distance_km
    print(f"road = {v:.1f} km: {((r < 1.03) | (r > 1.6)).mean():.0%} rows impossible")
print("impossible row (distance 1 km, road 15 km) ->",
      m.predict(Xte[:1].assign(distance_km=1.0, road_km=15.0)).round(1), "min")
fig, ax = plt.subplots(figsize=(6, 3.4))
ax.plot(z, ale, "o-", label="ALE (stays on the data)")
ax.plot(z, pdp, "s--", label="PDP (averages impossible rows)")
ax.set(xlabel="road_km", ylabel="centred effect (minutes)")
ax.legend(); fig.subplots_adjust(left=0.13, bottom=0.16)

degree-4 polynomial ridge: test MAE 2.88 min
grid : [ 0.6  1.9  3.2  4.6  5.8  7.1  8.3  9.5 10.7 11.8 15.2]
ALE  : [-13.5 -10.1  -6.7  -3.6  -0.8   1.4   3.4   4.9   6.5   8.1   7.3]
PDP  : [-19.8  -8.5  -1.8   1.9   4.2   6.    7.8   9.2   9.8   8.8 -17.6]
corr(distance, road) = 0.981
road = 0.6 km: 99% rows impossible
road = 15.2 km: 95% rows impossible
impossible row (distance 1 km, road 15 km) -> [10.3] min

৮ · শুরু থেকেই ব্যাখ্যাযোগ্য model আর global surrogate

দুটি পথ বাকি। প্রথমটি Rudin-এর পথ: black box-ই না বানানো। পর্ব ২৩-এ দেখেছেন linear model-এ interaction আর বানানো কলাম যোগ করা যায়। সব কলাম standardise করলে (standardised coefficient) সংখ্যাগুলো তুলনীয় হয়: “এক standard deviation বদলালে কত মিনিট”। কিন্তু দুটি সতর্কতা। Coefficient অন্য সব কলাম স্থির রেখে পড়া হয়, যা correlated কলামে আবার সেই অসম্ভব পরিস্থিতি। আর multicollinearity-তে coefficient অস্থির: r = ০.৯৮১ হলে variance inflation factor প্রায় ১/(১ − ০.৯৮১²) ≈ ২৭, অর্থাৎ ঐ দুটি coefficient-এর variance স্বাধীন কলামের তুলনায় প্রায় ২৭ গুণ।

দ্বিতীয় পথ: black box রেখে তার একটি পড়ার মতো নকল বানানো। Global surrogate হলো একটি ছোট, ব্যাখ্যাযোগ্য model (প্রায়ই পর্ব ৩২-এর ছোট decision tree), যা আসল label নয়, black box-এর অনুমান শেখে। এর মান মাপা হয় fidelity দিয়ে: held-out ডেটায় surrogate black box-কে কতটা নকল করে, যেমন black box-এর অনুমানের বিপরীতে R². চতুর্থ ব্লক দুটি পথ পাশাপাশি চালায়।

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor as HGB
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.tree import DecisionTreeRegressor, export_text
from sklearn.metrics import r2_score, mean_absolute_error as mae
from bazar_eta import split, COLS

Xtr, Xte, ytr, yte = split()
# ক) শুরু থেকেই ব্যাখ্যাযোগ্য: দুটি বানানো কলামসহ standardised linear model
def aug(X):
    return X.assign(road_x_rain=X.road_km * X.rain, rush=X.hour.between(17, 20) * 1)
for half in [slice(0, 1000), slice(1000, 2000)]:          # train-এর দুই অর্ধেক
    lin = make_pipeline(StandardScaler(), LinearRegression())
    lin.fit(aug(Xtr[half]), ytr[half])
    coef = pd.Series(lin[-1].coef_, index=aug(Xtr).columns).round(2)
    print(coef[["distance_km", "road_km", "road_x_rain", "rush"]].to_dict(),
          f"test MAE {mae(yte, lin.predict(aug(Xte))):.2f}")

# খ) post-hoc: black box-এর অনুমানের ওপর একটি ছোট গাছ (global surrogate)
m = HGB(max_iter=200, random_state=0).fit(Xtr, ytr)
print(f"black box test MAE {mae(yte, m.predict(Xte)):.2f}, R2 vs truth "
      f"{r2_score(yte, m.predict(Xte)):.3f}")
bb_tr, bb_te = m.predict(Xtr), m.predict(Xte)
for depth in [1, 2, 3, 4, 6, 8]:
    t = DecisionTreeRegressor(max_depth=depth, random_state=0).fit(Xtr, bb_tr)
    print(f"depth {depth}: leaves {t.get_n_leaves():3d}  fidelity R2 "
          f"{r2_score(bb_te, t.predict(Xte)):.3f}  R2 vs truth "
          f"{r2_score(yte, t.predict(Xte)):.3f}")
t3 = DecisionTreeRegressor(max_depth=3, random_state=0).fit(Xtr, bb_tr)
print(export_text(t3, feature_names=COLS, decimals=1))

{'distance_km': 4.69, 'road_km': 3.49, 'road_x_rain': 6.26, 'rush': 3.08} test MAE 2.52
{'distance_km': 4.63, 'road_km': 3.83, 'road_x_rain': 6.07, 'rush': 3.22} test MAE 2.49
black box test MAE 2.74, R2 vs truth 0.930
depth 1: leaves   2  fidelity R2 0.503  R2 vs truth 0.477
depth 2: leaves   4  fidelity R2 0.760  R2 vs truth 0.721
depth 3: leaves   8  fidelity R2 0.855  R2 vs truth 0.806
depth 4: leaves  16  fidelity R2 0.893  R2 vs truth 0.838
depth 6: leaves  64  fidelity R2 0.956  R2 vs truth 0.900
depth 8: leaves 250  fidelity R2 0.974  R2 vs truth 0.917
|--- distance_km <= 5.3
|   |--- road_km <= 3.9
|   |   |--- road_km <= 2.6
|   |   |   |--- value: [23.6]
|   |   |--- road_km >  2.6
|   |   |   |--- value: [28.2]
|   |--- road_km >  3.9
|   |   |--- rain <= 0.5
|   |   |   |--- value: [31.4]
|   |   |--- rain >  0.5
|   |   |   |--- value: [40.4]
|--- distance_km >  5.3
|   |--- rain <= 0.5
|   |   |--- distance_km <= 8.2
|   |   |   |--- value: [39.5]
|   |   |--- distance_km >  8.2
|   |   |   |--- value: [47.2]
|   |--- rain >  0.5
|   |   |--- road_km <= 9.1
|   |   |   |--- value: [51.4]
|   |   |--- road_km >  9.1
|   |   |   |--- value: [63.9]

৯ · AI-এর যুগে — ব্যাখ্যার কোড আর ব্যাখ্যার দাবি যাচাই

“আমার model-এর feature importance দেখাও” বললে AI assistant সাধারণত মিনিটে একটি সুন্দর bar chart দেয়। সমস্যা chart-এ নয়, তার নিচের সিদ্ধান্তগুলোতে: কোন ডেটায়, কোন score-এ, কতবার, আর correlated কলামের কী হলো। নিচের ভুলগুলো প্রায়ই চোখে পড়ে।

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

• Importance held-out ডেটায়, কয়েকবার এলোমেলো করে, sd-সহ মেপেছি
• একটি এলোমেলো কলাম যোগ করে দেখেছি, সেটি প্রায় শূন্য পায় কি না
• Correlated কলাম দলে বেঁধে grouped importance দেখেছি
• PDP-র পাশে ICE আর ALE দেখেছি, আর ডেটা যেখানে পাতলা সেখানে সতর্ক থেকেছি
• ছোট খেলনা ডেটায় একটি importance আর একটি PDP বিন্দু হাতে করে কোডের সঙ্গে মিলিয়েছি
• রিপোর্টে কোথাও “X-এর কারণে” লিখিনি, লিখেছি “model X-এর ওপর নির্ভর করে”