Linear model গভীরে — basis function, interaction আর regularization path

AI-এর ভিত্তি কোর্সের পর্ব ৫৮-এ normal equation, residual plot আর statsmodels-এর OLS দেখেছেন; পর্ব ৬১-এ PolynomialFeatures-এর ডিগ্রি বাড়িয়ে overfitting দেখেছেন; আর পর্ব ৬৩-এ log, ratio আর interaction feature বানিয়েছেন। এই কোর্সের পর্ব ১৭-এ ridge আর lasso-র তত্ত্ব আর path-এর প্রথম ছবি এসেছে। এই পর্ব সেগুলো জুড়ে একটি পূর্ণ model বানানোর পদ্ধতি দেখায়: বাঁকা সম্পর্ক কোন basis দিয়ে ধরবেন, coefficient কীভাবে পড়বেন, অনেক candidate থেকে কীভাবে বাছবেন, আর অনুমানের পাশে কোন পরিসর বলবেন। চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python); প্রথমবার scikit-learn আর statsmodels নামতে দশ সেকেন্ডের মতো লাগতে পারে।

১ · “Linear” মানে parameter-এ linear — basis expansion-এর মূল কথা

একটি ইনপুট x (ধরুন ডেলিভারির দূরত্ব) থেকে K-টি নতুন কলাম বানানো যাক: φ₁(x), φ₂(x), …, φ_K(x). Model:

ŷ = β₀ + β₁φ₁(x) + β₂φ₂(x) + … + β_Kφ_K(x)

এখানে φ_j হলো আগে থেকে ঠিক করা নিয়ম (যেমন x², বা একটি ঢিবি), আর β_j হলো শেখার মতো সংখ্যা। ŷ x-এর বাঁকা function হতে পারে, কিন্তু β-র দিক থেকে এটি একটি রৈখিক যোগফল। এর ফল অনেক বড়: design matrix Φ (n সারি, K কলাম) বানানোর পর বাকি সব পুরোনো হাতিয়ার হুবহু কাজ করে। ŵ = (ΦᵀΦ)⁻¹Φᵀy, পর্ব ০৭-এর projection, পর্ব ১৭-এর ridge আর lasso, আর statsmodels-এর standard error, সবই। এই কৌশলের নাম basis expansion. পর্ব ০৬-এর ভাষায়, φ-গুলো একটি function space-এর basis, আর model সেই space-এর মধ্যে সবচেয়ে কাছের function খোঁজে।

কোনটি linear model নয়? যেখানে parameter নিজেই বাঁকা জায়গায় বসে, যেমন ŷ = β₀·e^β₁x, বা knot-এর অবস্থানও শিখতে হয়। তখন loss আর quadratic থাকে না, আর closed form হারায়। Neural network-কে এভাবে দেখা যায়: শেষ স্তরটি basis-এর ওপর linear model, কিন্তু basis-ও শেখা হয়। এখানে basis আগে থেকে ঠিক, তাই সমস্যাটি convex থাকে (পর্ব ১২)।

২ · তিন রকম basis — polynomial, spline আর radial basis function

Polynomial. φ_j(x) = x^j: সবচেয়ে পরিচিত, polynomial features. দুটি সমস্যা। প্রথমত প্রতিটি কলাম পুরো x-রেখায় ছড়ানো, তাই ১০ km-এর একটি ডেটা বিন্দু ১ km-এর অনুমানও নাড়ায়। দ্বিতীয়ত উঁচু ডিগ্রিতে প্রান্তে বড় দোলা আসে। গণিতবিদ Carl Runge ১৯০১ সালে সমান দূরত্বের বিন্দুতে polynomial interpolation-এ এটি দেখান, তাই একে Runge-এর ঘটনা বলে। x, x², x³ … কলামগুলো পরস্পর খুব সম্পর্কিত, তাই ΦᵀΦ-এর condition number বিশাল; standardize না করলে x¹⁰-এর মান ১৫^১০ ছাড়ায়।

Spline. x-রেখাকে কয়েকটি knot-এ ভাগ করে প্রতিটি টুকরোয় ছোট ডিগ্রির polynomial বসানো, আর জোড়াগুলো মসৃণ রাখা। নামটি এসেছে নকশাকারদের নমনীয় কাঠের পাত থেকে, যা পেরেকের সঙ্গে বেঁকে মসৃণ বক্ররেখা আঁকত; Isaac Schoenberg ১৯৪৬ সালে গণিতে শব্দটি চালু করেন। একই spline অনেক basis-এ লেখা যায়। সহজটি truncated power basis: 1, x, আর প্রতিটি knot k-এর জন্য (x − k)₊, যা x < k হলে শূন্য। এর coefficient সরাসরি পড়া যায় (“৪ km-এর পর ঢাল কত বদলায়”), কিন্তু কলামগুলো খুব সম্পর্কিত। কম্পিউটারে পছন্দের রূপ B-spline: প্রতিটি কলাম মাত্র degree + 1টি টুকরোয় অশূন্য, আর প্রতিটি x-এ সবগুলোর যোগফল ১. এর হিসাবের স্থিতিশীল recursion Cox আর de Boor ১৯৭২ সালের কাছাকাছি সময়ে আলাদাভাবে প্রকাশ করেন। scikit-learn ১.০ থেকে SplineTransformer এই basis বানায়; default-এ ৫টি knot আর degree ৩, তাই প্রতি feature-এ ৫ + ৩ − ১ = ৭টি কলাম।

Radial basis function. φ_j(x) = exp(−((x − c_j)/w)²): c_j কেন্দ্রে একটি ঢিবি, w তার চওড়া। স্থানীয়, মসৃণ, আর অনেক মাত্রায় সহজে যায় (দূরত্ব ‖x − c‖ যেকোনো মাত্রায় মাপা যায়)। Broomhead আর Lowe ১৯৮৮ সালে RBF network নামে এটি ML-এ জনপ্রিয় করেন। প্রতিটি ডেটা বিন্দুকেই কেন্দ্র বানালে আর ridge যোগ করলে আসে kernel ridge regression; সেটি পর্ব ২৯ (Kernel trick — RBF, polynomial আর Mercer-এর অন্তর্দৃষ্টি) আর পর্ব ৩০ (SVR আর kernel ridge regression)-এর বিষয়।

নিচের stepper degree-১ B-spline দিয়ে হাতে সারি বানায়, তারপর একই function truncated power basis-এ লেখে, শেষে scikit-learn-এর default cubic basis-এর একটি সারি দেখায়।

৩ · কতটা নমনীয়? Cross-validation দিয়ে basis বাছাই

কলাম বাড়ালে training ভুল কমতেই থাকে, কিন্তু নতুন ডেটার ভুল একসময় বাড়ে; এটিই bias–variance trade-off, যা AI-এর ভিত্তি কোর্সের পর্ব ৬১-এ U-আকৃতির curve হিসেবে দেখেছেন। এখানে প্রশ্নটি একটু বড়: একই সংখ্যক কলামে কোন basis ভালো? নিচের ব্লক বাজারঘরের ৮০টি কাল্পনিক ডেলিভারি বানায়। শহরের ভেতরে (৫ km পর্যন্ত) প্রতি km-এ ৬ মিনিট, তারপর ১.৫ মিনিট, আর ১০ km-এর কাছে একটি ব্রিজের জট (প্রায় ১২ মিনিটের ঢিবি)। তিন রকম basis-এ ৪ থেকে ৩০টি কলাম, প্রতিবার একই pipeline-এ খুব ছোট ridge (α = 10⁻⁶, শুধু সংখ্যার স্থিতির জন্য), আর ৫-fold CV-র RMSE. RBF-এর জন্য scikit-learn-এ সরাসরি transformer নেই (RBFSampler অন্য জিনিস, random Fourier feature), তাই FunctionTransformer দিয়ে নিজেই বানানো হয়েছে।

import numpy as np, pandas as pd
import matplotlib.pyplot as plt
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import (PolynomialFeatures, SplineTransformer,
                                   StandardScaler, FunctionTransformer)
from sklearn.linear_model import Ridge
from sklearn.model_selection import KFold, cross_val_score

rng = np.random.default_rng(23)              # কাল্পনিক বাজারঘর: ৮০টি ডেলিভারি
d = rng.uniform(0.5, 15, 80)                 # দূরত্ব (km)
truth = lambda x: np.where(x < 5, 10 + 6 * x, 40 + 1.5 * (x - 5)) + 12 * np.exp(-((x - 10) / 0.8) ** 2)
y = truth(d) + rng.normal(0, 3, 80)          # ডেলিভারির মিনিট
X = d.reshape(-1, 1)
def rbf(k):                                  # k-টি Gaussian ঢিবি, সমান দূরে কেন্দ্র
    c = np.linspace(0.5, 15, k); w = c[1] - c[0]
    return FunctionTransformer(lambda Z: np.exp(-((Z - c) / w) ** 2))
make = {"poly":   lambda k: make_pipeline(StandardScaler(), PolynomialFeatures(k, include_bias=False), StandardScaler(), Ridge(1e-6)),
        "spline": lambda k: make_pipeline(SplineTransformer(n_knots=k - 2, degree=3), Ridge(1e-6)),
        "rbf":    lambda k: make_pipeline(rbf(k), Ridge(1e-6))}
cv = KFold(5, shuffle=True, random_state=42)
ks = list(range(4, 31, 2))                   # কলামের সংখ্যা: spline-এ n_knots + 2
res = pd.DataFrame({name: [-cross_val_score(m(k), X, y, cv=cv, scoring="neg_root_mean_squared_error").mean()
                           for k in ks] for name, m in make.items()}, index=ks)
res.index.name = "columns"
print(res.loc[[4, 8, 12, 16, 20, 24, 30]].round(2))
print("best k:", res.idxmin().to_dict(), " best RMSE:", res.min().round(2).to_dict())
fig, ax = plt.subplots(1, 2, figsize=(9, 3.6), layout="constrained")
res.plot(ax=ax[0], marker="o"); ax[0].set(ylim=(3, 12), xlabel="number of basis columns", ylabel="5-fold CV RMSE (min)")
g = np.linspace(0.5, 15, 300).reshape(-1, 1)
ax[1].scatter(d, y, s=8, c="0.6")
for name, k in [("poly", 16), ("spline", res["spline"].idxmin()), ("rbf", res["rbf"].idxmin())]:
    ax[1].plot(g, make[name](k).fit(X, y).predict(g), lw=2, label=f"{name}, k={k}")
ax[1].plot(g, truth(g), "k--", lw=1, label="truth"); ax[1].legend(fontsize=8)
ax[1].set(xlabel="distance (km)", ylabel="delivery minutes"); plt.show()

         poly  spline    rbf
columns
4        3.92    3.92   4.16
8        3.57    3.81   3.65
12       3.64    3.42   3.33
16       7.26    3.47   3.46
20       9.87    3.59   3.54
24      11.28    5.12   4.09
30      15.64    4.72  33.42
best k: {'poly': 8, 'spline': 14, 'rbf': 12}  best RMSE: {'poly': 3.57, 'spline': 3.35, 'rbf': 3.33}

ব্যবহারিক দুটি পরামর্শ। প্রথম, knot গুনে বাছাই করার চেয়ে প্রায়ই ভালো হলো অনেক knot রেখে ridge-ধাঁচের শাস্তি দেওয়া, আর শাস্তির মাপ λ CV-তে বাছা। Paul Eilers আর Brian Marx ১৯৯৬ সালে এই P-spline প্রস্তাব করেন: পাশাপাশি B-spline coefficient-এর পার্থক্যের ওপর শাস্তি, যাতে বক্ররেখা মসৃণ থাকে। তখন নমনীয়তার মাপ হয় effective degrees of freedom, ঠিক পর্ব ১৭-এর মতো। দ্বিতীয়, knot-এর বাইরে কী হবে আগে ভাবুন। SplineTransformer-এর extrapolation default-এ "constant", তাই ১৫ km-এর বেশি দূরত্বে অনুমান প্রান্তের মানে আটকে যায়; "linear" দিলে শেষ ঢাল ধরে চলে। ঘণ্টার মতো চক্রাকার feature-এ extrapolation="periodic" ২৩টা আর ০টাকে পাশাপাশি রাখে।

৪ · Transform-এর পরে coefficient পড়া, আর GAM-এর ধারণা

Basis বদলালে coefficient বদলায় কিন্তু অনুমান বদলায় না, এটি stepper-এ দেখেছেন। তাই রূপান্তরের পরে “এই coefficient মানে কী” প্রশ্নের উত্তর সবসময় রূপান্তর জেনে দিতে হয়। পাঁচটি সাধারণ ক্ষেত্র:

• Spline বা RBF: একক coefficient প্রায় অর্থহীন; cubic B-spline-এর β_j কোনো সহজ ঢাল বা মান নয়। পড়ুন বক্ররেখা: একটি grid-এ অনুমান আঁকুন, বা দুটি বিন্দুর পার্থক্য বলুন (“৩ থেকে ৪ km-এ কত মিনিট বাড়ে”)। বিশ্বাসযোগ্যতার জন্য বক্ররেখার পাশে CI বা bootstrap ব্যান্ড দিন (পর্ব ২১)।
• Polynomial: ŷ = β₀ + β₁x + β₂x² হলে x-এর প্রভাব β₁ + 2β₂x, অর্থাৎ x-এর মানের ওপর নির্ভর করে। β₁ একা মানে “x = ০-তে ঢাল”, যা হয়তো ডেটার বাইরে।
• Log target: log(অর্ডারের মূল্য)-এ promo-র β = ০.১২ মানে মূল্য গুণ হয় e^০.১২ = ১.১২৭৫ দিয়ে, অর্থাৎ প্রায় ১২.৭৫% বেশি, টাকার কোনো স্থির অঙ্ক নয়। সাবধান: exp(log-এর অনুমান) দেয় median-ধাঁচের মান, গড় নয়; গড় চাইলে Duan (১৯৮৩)-এর “smearing” সংশোধন লাগে।
• Log feature: y-তে log x-এর coefficient β মানে x ১% বাড়লে y প্রায় β/১০০ একক বাড়ে।
• Standardized feature: standardization-এর পরে β মানে “এক standard deviation বাড়লে”; মূল এককে ফিরতে β-কে ওই কলামের sd দিয়ে ভাগ করুন। Pipeline-এর শেষ ধাপের coef_ সরাসরি km বা টাকায় নয়।

অনেক feature-এ বাঁকা সম্পর্ক: generalized additive model

দূরত্ব ছাড়াও ডেলিভারি সময় দিনের ঘণ্টা (সন্ধ্যার ভিড়) আর বৃষ্টির ওপর নির্ভর করে। প্রতিটির জন্য আলাদা একটি মসৃণ function শিখে যোগ করলে পাওয়া যায় generalized additive model (GAM): ŷ = β₀ + f₁(দূরত্ব) + f₂(ঘণ্টা) + β·বৃষ্টি. প্রতিটি f_j একটি penalized spline, তাই পুরো model এখনো parameter-এ linear, শুধু শাস্তিসহ। Trevor Hastie আর Robert Tibshirani ১৯৮৬ সালে ধারণাটি প্রকাশ করেন, ১৯৯০ সালে বই; R-এর mgcv package (Simon Wood) আজ এর সবচেয়ে প্রচলিত বাস্তবায়ন। GAM-এর বড় সুবিধা পড়তে পারা: প্রতিটি f_j আলাদা করে আঁকা যায়। দাম হলো, feature-দের মধ্যে interaction নিজে থেকে আসে না; চাইলে আলাদা করে যোগ করতে হয় (যেমন দূরত্ব × বৃষ্টির tensor spline).

Python-এ pygam package এটি দেয়, কিন্তু ব্রাউজারের Python-এ নেই। তাই ডেস্কটপে সত্যিকারের run: লেখার সময়ে Python 3.14.4, pygam 0.12.0, NumPy 2.5.4 আর SciPy 1.16.3, একটি venv-এ। ডেটা উপরের মতোই, শুধু ৬০০টি অর্ডার আর দুটি নতুন feature:

import numpy as np
from pygam import LinearGAM, s, f

rng = np.random.default_rng(23)                 # কাল্পনিক বাজারঘর: ৬০০টি ডেলিভারি
n = 600
km = rng.uniform(0.5, 15, n)
hour = rng.uniform(8, 23, n)
rain = rng.integers(0, 2, n)
f_km = np.where(km < 5, 10 + 6 * km, 40 + 1.5 * (km - 5)) + 12 * np.exp(-((km - 10) / 0.8) ** 2)
f_hour = 8 * np.exp(-((hour - 19) / 1.5) ** 2)  # সন্ধ্যার ভিড়
y = f_km + f_hour + 6 * rain + rng.normal(0, 3, n)
X = np.column_stack([km, hour, rain])

gam = LinearGAM(s(0, n_splines=20) + s(1, n_splines=20) + f(2)).gridsearch(X, y, progress=False)
gam.summary()
for x in (3.0, 9.0, 10.0, 13.0):                 # দূরত্বের partial effect: ১ km বাড়লে কত মিনিট?
    a, b = gam.partial_dependence(term=0, X=np.array([[x, 0, 0], [x + 1, 0, 0]]))
    print(f"km {x:>4} -> {x + 1:>4}: {b - a:+.1f} min")

~/bazarghar-ml$ python gam_delivery.py
…/gam_delivery.py:15: UserWarning: KNOWN BUG: p-values computed in this summary are likely much smaller than they should be.
…
LinearGAM
=============================================== ==========================================================
Distribution:                        NormalDist Effective DoF:                                     29.8614
Link Function:                     IdentityLink Log Likelihood:                                 -1484.4614
Number of Samples:                          600 AIC:                                             3030.6455
                                                AICc:                                             3034.107
                                                GCV:                                                9.5209
                                                Scale:                                              2.9448
                                                Pseudo R-Squared:                                   0.9518
==========================================================================================================
Feature Function                  Lambda               Rank         EDoF         P > x        Sig. Code
================================= ==================== ============ ============ ============ ============
s(0)                              [0.2512]             20           15.1         1.11e-16     ***
s(1)                              [0.2512]             20           13.8         1.11e-16     ***
f(2)                              [0.2512]             2            1.0          1.11e-16     ***
intercept                                              1            0.0          1.11e-16     ***
==========================================================================================================
… (Significance codes-এর লাইন আর p-value নিয়ে দুটি WARNING অনুচ্ছেদ)
km  3.0 ->  4.0: +6.9 min
km  9.0 -> 10.0: +8.8 min
km 10.0 -> 11.0: -5.2 min
km 13.0 -> 14.0: +1.6 min

প্রতিটি অংশের rank (কলাম) ২০, কিন্তু শাস্তির পরে effective df দূরত্বে ১৫.১ আর ঘণ্টায় ১৩.৮; পুরো model প্রায় ৩০টি parameter-এর সমান নমনীয়। gridsearch GCV (generalized cross-validation) কমিয়ে সব অংশে একই λ = ০.২৫১২ বেছেছে। শেষ চার লাইন coefficient না পড়ে বক্ররেখা পড়ার উদাহরণ: ৩ থেকে ৪ km-এ +৬.৯ মিনিট (আসল ৬), ১৩ থেকে ১৪-এ +১.৬ (আসল ১.৫), আর ব্রিজের দুই পাশে +৮.৮ আর −৫.২ (আসল প্রায় +১১ আর −৮): শাস্তি ঢিবিটিকে একটু চ্যাপ্টা করেছে, যা শাস্তির স্বাভাবিক bias. pygam নিজেই সতর্ক করছে যে এই summary-র p-value (সব ১.১১e-16) বিশ্বাসযোগ্য নয়, কারণ λ ডেটা থেকে বাছা হয়েছে; এই ধরনের সতর্কবার্তা চাপা না দিয়ে পড়ুন। scikit-learn-এ কাছাকাছি রূপ: ColumnTransformer-এ প্রতিটি কলামে SplineTransformer, তারপর RidgeCV.

৫ · Interaction আর hierarchy — একটি coefficient-এর মানে কীভাবে বদলায়

Basis expansion একটি feature-কে বাঁকায়; interaction দুটি feature-কে জোড়ে। বাজারঘরের প্রশ্ন: promo code কি প্রথমবারের ক্রেতার অর্ডার বেশি বাড়ায়? Model ŷ = b₀ + b₁·promo + b₂·first + b₃·(promo × first), যেখানে দুটিই ০/১. promo × first কলামটি কেবল “প্রথমবারের ক্রেতা আর promo পেয়েছেন” ঘরে ১. ফলে b₁-এর মানে আর “promo-র প্রভাব” নয়, হয়ে যায় “first = ০ হলে promo-র প্রভাব”. নিচের stepper চারটি ঘর একে একে ভরে।

Hierarchy নিয়ম: model-এ promo × first রাখলে promo আর first-এর main effect-ও রাখুন। John Nelder ১৯৭০-এর দশকে একে “marginality” নামে জোর দিয়ে বলেছেন। কারণ: main effect বাদ দিলে আপনি গোপনে একটি শর্ত চাপান (যেমন “promo ছাড়া প্রথমবার আর ফেরত আসা ক্রেতার গড় সমান”), আর interaction-এর coefficient সেই ভুল শর্তের ক্ষতিপূরণ করতে গিয়ে অন্য মানে নেয়। তা ছাড়া কোডিং বদলালে (০/১-এর বদলে ১/০, বা center করলে) hierarchy-ছাড়া model-এর অনুমানই বদলে যায়, যা hierarchy-সহ model-এ হয় না। নিচের ব্লক statsmodels-এর OLS দিয়ে চারটি model তুলনা করে।

import numpy as np
import pandas as pd
import statsmodels.api as sm

rng = np.random.default_rng(7)               # কাল্পনিক বাজারঘর: ২০০০টি অর্ডার
n = 2000
promo = rng.integers(0, 2, n)                # promo code দেওয়া হয়েছিল কি না
first = (rng.random(n) < 0.3).astype(int)    # প্রথমবারের ক্রেতা কি না (~30%)
value = 850 + 120 * promo - 200 * first + 180 * promo * first + rng.normal(0, 150, n)
df = pd.DataFrame({"promo": promo, "first": first, "value": value})
print(df.pivot_table("value", "first", "promo").round(0))   # চারটি ঘরের গড়

def fit(cols, label):
    X = sm.add_constant(df[cols].astype(float))
    r = sm.OLS(df["value"], X).fit()
    print(f"\n{label}\n" + r.params.round(1).to_string() + f"\n(SE of promo: {r.bse['promo']:.1f})")

df["promo_x_first"] = df["promo"] * df["first"]
fit(["promo", "first", "promo_x_first"], "A. hierarchy-সহ পূর্ণ model")
df["first_c"] = df["first"] - df["first"].mean()      # first-কে center করা
df["promo_x_first_c"] = df["promo"] * df["first_c"]
fit(["promo", "first_c", "promo_x_first_c"], "B. first center করে")
fit(["promo", "first"], "C. interaction ছাড়া")
fit(["promo", "promo_x_first"], "D. first-এর main effect বাদ (hierarchy ভাঙা)")

promo      0      1
first
0      853.0  970.0
1      650.0  939.0

A. hierarchy-সহ পূর্ণ model
const            852.8
promo            117.0
first           -203.1
promo_x_first    172.3
(SE of promo: 8.0)

B. first center করে
const              792.7
promo              168.0
first_c           -203.1
promo_x_first_c    172.3
(SE of promo: 6.7)

C. interaction ছাড়া
const    827.3
promo    168.0
first   -115.1
(SE of promo: 6.9)

D. first-এর main effect বাদ (hierarchy ভাঙা)
const            793.9
promo            175.9
promo_x_first    -30.8
(SE of promo: 8.0)

সংখ্যার feature-এও একই কথা: দূরত্ব × বৃষ্টি রাখলে বৃষ্টির main effect মানে “০ km-এ বৃষ্টির প্রভাব”; দূরত্ব center করলে “গড় দূরত্বে”. আর বেশি interaction-এ সাবধান: ১০টি feature থেকে জোড়ায় জোড়ায় ৪৫টি interaction, তিনটি করে ১২০টি। সব যোগ করলে model দ্রুত noise ধরতে শুরু করে; পরের অংশ এই সমস্যাটিই সমাধান করে।

৬ · অনেক candidate feature — regularization path আর one-standard-error নিয়ম

বাস্তবে প্রশ্নটি প্রায়ই এমন: ১০টি মূল feature আছে, তাদের জোড়ায় জোড়ায় interaction মিলে ৫৫টি candidate, আর ৩০০টি training সারি। কোনগুলো রাখবেন? পর্ব ১৭-এ দেখেছেন, lasso আর elastic net-এর α বড় থেকে ছোট করলে একটি একটি করে feature ঢোকে; এই পুরো যাত্রাই regularization path. ElasticNetCV প্রতিটি fold-এ পুরো path চালায় (আগের α-র সমাধান থেকে শুরু করে, যাকে warm start বলে), আর প্রতিটি α-র validation MSE রাখে mse_path_-এ। তারপর প্রশ্ন: কোন α?

সবচেয়ে কম গড় CV ভুলের α স্বাভাবিক উত্তর, কিন্তু CV ভুল নিজেই অনিশ্চিত; কাছাকাছি কয়েকটি α-র পার্থক্য প্রায়ই noise. One-standard-error নিয়ম বলে: সর্বনিম্ন গড় ভুলের সঙ্গে তার standard error যোগ করে একটি সীমা ধরুন, আর সেই সীমার নিচে থাকা সবচেয়ে সরল model (সবচেয়ে বড় α) নিন। Breiman, Friedman, Olshen আর Stone-এর CART বই (১৯৮৪) tree ছাঁটাইয়ের জন্য এই নিয়ম দেয়; পরে Hastie, Tibshirani আর Friedman-এর বই আর R-এর glmnet (lambda.1se) একে lasso-তে জনপ্রিয় করে। এটি একটি প্রথা, উপপাদ্য নয়।

এবার আসল কাজ। নিচের ব্লকে ৪০০টি কাল্পনিক ডেলিভারি, ১০টি মূল feature, যার মাত্র চারটি (দূরত্ব, পণ্যের সংখ্যা, বৃষ্টি, ভিড়ের সময়) আর একটি interaction (বৃষ্টি × দূরত্ব) সময়কে প্রভাবিত করে। বাকিগুলো (promo, প্রথমবারের ক্রেতা, সপ্তাহান্ত, রাইডারের অভিজ্ঞতা, নতুন app, আর একটি নিছক noise) কিছুই করে না। ২৫% সারি শুরুতেই test হিসেবে তালাবন্ধ। Pipeline: standardize (center), তারপর জোড়ায় জোড়ায় গুণফল, তারপর আবার scale, তারপর ElasticNetCV (l1_ratio ০.৫ আর ১.০, ৫০টি α)। alphas=50 লেখা হয়েছে, কারণ n_alphas scikit-learn ১.৭-এ deprecated হয়েছে আর ১.৯-এ বাদ যাবে; max_iter=50000 যাতে কোনো ConvergenceWarning না আসে।

import numpy as np
import matplotlib.pyplot as plt
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import ElasticNetCV, ElasticNet
from sklearn.model_selection import KFold, train_test_split

rng = np.random.default_rng(5)               # কাল্পনিক বাজারঘর: ৪০০টি ডেলিভারি, ১০টি মূল feature
names = ["dist", "items", "rain", "peak", "promo", "first", "weekend", "rider_exp", "app_new", "noise"]
n = 400
B = np.column_stack([rng.uniform(1, 12, n), rng.integers(1, 9, n)] + [rng.integers(0, 2, n) for _ in range(5)]
                    + [rng.uniform(0, 36, n), rng.integers(0, 2, n), rng.normal(size=n)]).astype(float)
y = 12 + 3 * B[:, 0] + 1.5 * B[:, 1] + 5 * B[:, 2] + 6 * B[:, 3] + 0.8 * B[:, 2] * B[:, 0] + rng.normal(0, 4, n)
B_tr, B_te, y_tr, y_te = train_test_split(B, y, test_size=0.25, random_state=42)   # test আগে আলাদা
poly = PolynomialFeatures(2, interaction_only=True, include_bias=False)
feat = poly.fit(B_tr).get_feature_names_out(names)           # 10 + 45 = 55টি candidate
pipe = lambda est: make_pipeline(StandardScaler(), poly, StandardScaler(), est)   # center → গুণ → scale
cv = KFold(5, shuffle=True, random_state=42)
m = pipe(ElasticNetCV(l1_ratio=[0.5, 1.0], alphas=50, cv=cv, max_iter=50000)).fit(B_tr, y_tr)[-1]
i = list(m.l1_ratio).index(m.l1_ratio_)
mse, alphas = m.mse_path_[i], m.alphas_[i]                    # (50 alpha × 5 fold), alpha বড় → ছোট
mean, se = mse.mean(1), mse.std(1, ddof=1) / np.sqrt(mse.shape[1])
j = mean.argmin(); j1 = np.flatnonzero(mean <= mean[j] + se[j]).min()   # 1-SE: সবচেয়ে বড় alpha
print(f"l1_ratio={m.l1_ratio_}  alpha_min={alphas[j]:.3f}  alpha_1se={alphas[j1]:.3f}")
print(f"CV MSE: min {mean[j]:.2f} ± SE {se[j]:.2f}, threshold {mean[j] + se[j]:.2f}, at 1se {mean[j1]:.2f}")
for tag, a in [("min", alphas[j]), ("1se", alphas[j1])]:
    p = pipe(ElasticNet(alpha=a, l1_ratio=m.l1_ratio_, max_iter=50000)).fit(B_tr, y_tr)
    nz = np.flatnonzero(p[-1].coef_)
    rmse = np.sqrt(np.mean((p.predict(B_te) - y_te) ** 2))
    print(f"\n{tag}: {len(nz)} of {len(feat)} features, test RMSE {rmse:.2f}\n" + ", ".join(feat[nz]))
plt.figure(figsize=(6, 3.4), layout="constrained")
plt.errorbar(np.log10(alphas), mean, yerr=se, fmt="o-", ms=3, capsize=2)
plt.axvline(np.log10(alphas[j]), ls="--", c="k", label="min CV")
plt.axvline(np.log10(alphas[j1]), c="r", label="1-SE rule")
plt.yscale("log"); plt.xlabel("log10(alpha)"); plt.ylabel("5-fold CV MSE (log scale)"); plt.legend(); plt.show()

l1_ratio=1.0  alpha_min=0.395  alpha_1se=0.695
CV MSE: min 17.55 ± SE 1.64, threshold 19.20, at 1se 18.66

min: 8 of 55 features, test RMSE 3.69
dist, items, rain, peak, weekend, rider_exp, dist rain, rain rider_exp

1se: 5 of 55 features, test RMSE 3.92
dist, items, rain, peak, dist rain

একটি ছোট পরিবর্তন অনেক কিছু বদলায়। ডেস্কটপে একটি স্থানীয় run-এ (scikit-learn 1.8.0) একই ডেটায় শুরুর standardize ধাপটি বাদ দিয়ে, অর্থাৎ center না করেই গুণফল বানিয়ে, one-SE model ১১টি feature রেখেছে, যার ছয়টি ভুয়া interaction (যেমন items × rain, peak × weekend), আর test RMSE ৪.০৫. কারণ: center না করা গুণফল rain × dist নিজের উপাদান rain-এর সঙ্গে খুব সম্পর্কিত, আর lasso সম্পর্কিত কলামের মধ্যে খেয়ালখুশিমতো ওজন ভাগ করে (পর্ব ১৭-এর bootstrap পরীক্ষা মনে করুন)। আরেকটি সমস্যা: lasso hierarchy মানে না, main effect বাদ দিয়ে interaction রাখতে পারে। তখন হয় বাদ পড়া main effect-গুলো ফিরিয়ে আনুন, নয়তো hierarchy মানে এমন পদ্ধতি নিন (Bien, Taylor আর Tibshirani-র ২০১৩ সালের hierNet, বা group lasso).

৭ · Confidence interval বনাম prediction interval — গড়ের পরিসর আর একটি অর্ডারের পরিসর

একটি নতুন অর্ডারের ইনপুট x₀ (সামনে ১ বসানো সারি)। অনুমান ŷ₀ = x₀ᵀβ̂. দুটি প্রশ্ন আলাদা: (ক) এই দূরত্বের সব অর্ডারের গড় সময় কোথায়? (খ) এই অর্ডারটির সময় কোথায় পড়বে? Residual-এর standard deviation s, আর h₀ = x₀ᵀ(XᵀX)⁻¹x₀ (x₀ ডেটার কেন্দ্র থেকে কত দূরে, তার মাপ; পর্ব ১৭-এর hat matrix-এর একটি ঘর)। তাহলে:

CI: ŷ₀ ± t_n−p · s√h₀     PI: ŷ₀ ± t_n−p · s√(1 + h₀)

এখানে t_n−p হলো n − p degrees of freedom-এর t বণ্টনের ৯৭.৫ শতাংশ বিন্দু (৯৫% পরিসরের জন্য)। এক feature-এর regression-এ h₀ = 1/n + (x₀ − x̄)²/∑(x_i − x̄)². Confidence interval-এ কেবল β̂-র অনিশ্চয়তা, যা n বাড়লে ১/√n হারে শূন্যের দিকে যায়। Prediction interval-এ বাড়তি “১”, অর্থাৎ নতুন অর্ডারের নিজের noise, যা n যত বড়ই হোক থাকে। তাই n → ∞-এ CI-এর চওড়া → ০, কিন্তু PI-এর চওড়া → 2 × 1.96 × σ. দুটি সূত্রই ধরে নেয় model ঠিক, noise স্বাধীন আর সব x-এ একই variance (homoscedastic), আর PI-এর জন্য noise প্রায় normal. নিচের ব্লক statsmodels-এর get_prediction দিয়ে দুটি পরিসর বানায়, হাতে সূত্র মেলায়, আর ৫০০০টি ভবিষ্যৎ অর্ডারে যাচাই করে।

import numpy as np
import pandas as pd
import statsmodels.api as sm

def orders(n, seed):                          # কাল্পনিক বাজারঘর: শহরের ভেতরের ডেলিভারি
    r = np.random.default_rng(seed)
    km = r.uniform(1, 8, n)
    return km, 10 + 5 * km + r.normal(0, 6, n)    # আসল গড় 10 + 5·km, noise-এর sd 6
fit = lambda k, yy: sm.OLS(yy, sm.add_constant(k)).fit()

km, y = orders(100, 1)                        # ১০০টি পুরোনো অর্ডারে fit
res = fit(km, y)
x0 = np.array([[1.0, 5.0]])                   # নতুন অর্ডার: ৫ km
print(res.get_prediction(x0).summary_frame(alpha=0.05).round(2).T.to_string(header=False))
h = (x0 @ np.linalg.inv(res.model.exog.T @ res.model.exog) @ x0.T).item()
s = np.sqrt(res.scale)                        # residual-এর standard deviation
print(f"s = {s:.2f}  h0 = {h:.4f}  SE(mean) = {s*np.sqrt(h):.2f}  SE(new) = {s*np.sqrt(1 + h):.2f}\n")

km_new, y_new = orders(5000, 99)              # ভবিষ্যতের ৫০০০টি অর্ডার (held-out)
rows = []
for n in (20, 100, 1000):
    f = fit(*orders(n, 1)).get_prediction(sm.add_constant(km_new)).summary_frame(alpha=0.05)
    rows.append({"n": n, "CI width": (f.mean_ci_upper - f.mean_ci_lower).mean(),
                 "PI width": (f.obs_ci_upper - f.obs_ci_lower).mean(),
                 "new y in CI": ((f.mean_ci_lower <= y_new) & (y_new <= f.mean_ci_upper)).mean(),
                 "new y in PI": ((f.obs_ci_lower <= y_new) & (y_new <= f.obs_ci_upper)).mean()})
print(pd.DataFrame(rows).round(3).to_string(index=False))

hit_ci = hit_pi = 0                           # ৫০০ বার নতুন training set: ৫ km-এ কতবার ঢাকে?
for rep in range(500):
    fr = fit(*orders(100, 1000 + rep)).get_prediction(x0).summary_frame(alpha=0.05).iloc[0]
    y5 = 35 + np.random.default_rng(rep).normal(0, 6)       # ৫ km-এর একটি নতুন অর্ডার
    hit_ci += fr.mean_ci_lower <= 35 <= fr.mean_ci_upper
    hit_pi += fr.obs_ci_lower <= y5 <= fr.obs_ci_upper
print(f"\n500 repeats: CI covers true mean 35 in {hit_ci/500:.3f}, PI covers a new order in {hit_pi/500:.3f}")

mean           34.65
mean_se         0.61
mean_ci_lower  33.44
mean_ci_upper  35.86
obs_ci_lower   22.71
obs_ci_upper   46.59
s = 5.99  h0 = 0.0104  SE(mean) = 0.61  SE(new) = 6.02

   n  CI width  PI width  new y in CI  new y in PI
  20     9.766    31.878        0.564        0.988
 100     3.286    24.000        0.220        0.957
1000     1.055    24.268        0.071        0.959

500 repeats: CI covers true mean 35 in 0.960, PI covers a new order in 0.966

এই সুন্দর ফলাফলের পেছনে ডেটা নিজেই সূত্রের অনুমান মেনে বানানো। বাস্তবে ডেলিভারির ভুল দূরত্বের সঙ্গে বাড়ে (heteroscedasticity), আর কিছু অর্ডার হঠাৎ অনেক দেরি হয়। তখন একই চওড়া PI কাছের অর্ডারে বেশি চওড়া আর দূরের অর্ডারে কম চওড়া হয়, আর সামগ্রিক ৯৫% ঠিক থাকলেও প্রতিটি দলে ভুল হয়। এর দুটি প্রতিকার আসছে: পর্ব ২৪ (Robust আর quantile regression — outlier আর অনুমানের পরিসর) সরাসরি ১০ম আর ৯০তম percentile শেখে, আর পর্ব ৫৮ (Conformal prediction) প্রায় কোনো অনুমান ছাড়াই coverage-এর নিশ্চয়তা দেয়। আর ridge বা lasso দিয়ে fit করা model-এ উপরের সূত্র খাটে না, কারণ shrinkage-এর bias আছে; সেখানে bootstrap (পর্ব ২১) বা conformal পদ্ধতি লাগে।

৮ · AI-এর যুগে — basis, interaction আর interval-এর দাবি যাচাই

“একটু polynomial যোগ করো”, “সব interaction দিয়ে lasso চালাও”, “৯৫% interval দাও”: AI assistant এই তিনটি কাজ দ্রুত করে দেয়, আর কোড প্রায় সবসময় চলে। ভুল থাকে অর্থে: কোন পরিসর, কোন coefficient, কোন ডেটায় বাছাই। নিচের টেবিলে সাধারণ ভুলগুলো।

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

• Test set শুরুতেই আলাদা ছিল; basis, scaler আর বাছাই সব pipeline-এর ভেতরে, CV-র প্রতিটি fold-এ আলাদা fit
• ডিগ্রি, knot বা λ CV-তে বেছেছি, আর একটি সরল baseline (শুধু সরলরেখা) কয়েকটি seed বা fold-এ তুলনা করেছি
• প্রতিটি interaction-এর main effect আছে; coefficient-এর অর্থ রূপান্তর আর centering জেনে লিখেছি
• One-SE নিয়মের হিসাব একটি ছোট টেবিলে হাতে মিলিয়েছি, আর alphas_-এর ক্রম দেখেছি
• একক অনুমানের পাশে PI দিয়েছি, CI নয়, আর held-out ডেটায় coverage মেপেছি
• Library-র parameter নিজের version-এর ডকুমেন্টেশনে আছে কি না দেখেছি; কোনো warning চাপা দিইনি
• কোনো আসল ক্রেতার ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি