Linear regression — loss, gradient descent আর closed form

এই কোর্সে রেখা fit করা নতুন নয়। পর্ব ১৭-এ matrix, solve, condition number আর np.linalg.lstsq এসেছে; পর্ব ৩৪-এ বাজারঘরের ৮টি অর্ডারে এক আর দুই parameter-এর রেখা gradient descent দিয়ে fit করে lstsq-এর সঙ্গে মেলানো হয়েছে। এই পর্ব সেখান থেকে এগোয়: অনেক feature, আসল অর্থে একটি linear regression model, আর সেই model কতটা বিশ্বাস করা যায় তা যাচাইয়ের অভ্যাস। কোড পাতাতেই চলে (ব্রাউজারের ভেতরে Python, Pyodide); scikit-learn বা statsmodels প্রথমবার লোড হতে প্রায় ১০ সেকেন্ড বা তার বেশি লাগতে পারে, পরের run দ্রুত।

১ · Model-টি — ŷ = w·x + b থেকে ŷ = Xw

পর্ব ৩৪-এর model ছিল এক লাইনের: দিন ≈ w × দূরত্ব + b. এখানে w হলো রেখার ঢাল (দূরত্ব এক একক বাড়লে দিন কত বাড়ে), আর b হলো intercept, দূরত্ব শূন্য হলে রেখা যেখানে y-অক্ষ ছোঁয়। কিন্তু ডেলিভারির সময় শুধু দূরত্বে নির্ভর করে না। পার্সেল ভারী হলে বাইকে তুলতে সময় লাগে, পণ্য বেশি হলে মেলাতে সময় লাগে, বৃষ্টি হলে রাস্তা ধীর, অফিস ছুটির ভিড়ে আরও ধীর। প্রতিটি কারণ একটি feature, আর প্রতিটির জন্য একটি করে weight:

ŷ = b + w₁x₁ + w₂x₂ + … + w_px_p

ŷ (পড়া হয় “y-hat”) হলো model-এর অনুমান, x₁ … x_p হলো p-টি feature-এর মান, আর w₁ … w_p হলো তাদের coefficient. চায়ের দোকানের হিসাবের খাতা ভাবুন: “চা ১০ টাকা করে, বিস্কুট ৫ টাকা করে, আর বসার জন্য কিছু না” — মোট বিল = 10 × চা + 5 × বিস্কুট + 0. Linear regression ঠিক এই রকম একটি দামের তালিকা খোঁজে, শুধু “দাম”-গুলো আগে জানা নেই, পুরোনো বিল থেকে বের করতে হয়। “Linear” মানে w-গুলোর ওপর রৈখিক: প্রতিটি w একবার করে, একটি সংখ্যা দিয়ে গুণ হয়ে যোগ হয়।

এখন কৌশলটি, যা পর্ব ৩৪-এর শেষ কোডেও ছিল (column_stack([x, ones])): প্রতিটি অর্ডারের সারিতে সবার আগে একটি বাড়তি কলাম বসাই, যার মান সবসময় 1. তাহলে b-ও একটি সাধারণ weight হয়ে যায় (1 × b = b), আর সব অর্ডারের অনুমান একসঙ্গে একটি matrix গুণ:

ŷ = X w · X-এর shape (n, p + 1), w-এর (p + 1,), ŷ-এর (n,)

এই এক লাইনের সুবিধা বিশাল। Feature পাঁচটি হোক বা পাঁচশো, সূত্র একই; NumPy-তে একটি @, কোনো loop নেই (পর্ব ১৫-এর vectorization). পরের সব অংশে X বললে এই bias কলামসহ matrix বুঝবেন।

২ · বাজারঘরের কাজ, আর loss কেন বর্গ

বাজারঘর (এই কোর্সের কাল্পনিক অনলাইন দোকান) শহরের ভেতরের অর্ডারে ক্রেতাকে বলতে চায়, “আপনার পার্সেল আনুমানিক কত মিনিটে পৌঁছাবে।” আমরা ২০০টি অর্ডারের ডেটা নিজেরাই বানাব, একটি জানা নিয়ম দিয়ে, যাতে পরে দেখা যায় পদ্ধতিগুলো সেই নিয়ম ফিরে পায় কি না। পরীক্ষার খাতা দেখার আগে শিক্ষক যেমন উত্তরপত্র নিজে জানেন, তেমন।

সংখ্যাগুলো এই পর্বের জন্য বানানো, কোনো আসল কুরিয়ারের নয়। বাস্তবে “আসল coefficient” কেউ জানে না; এখানে জানি বলেই পদ্ধতি যাচাই করা যায়। এখন প্রতিটি অর্ডারের ভুল বা residual r_i = y_i − ŷ_i, আর loss হলো MSE:

MSE(w) = (1/n) · ∑_i (y_i − ŷ_i)² = (1/n) · ‖y − Xw‖²

‖ ‖ হলো vector-এর দৈর্ঘ্য বা L2 norm (পর্ব ১৭)। বর্গ কেন, পরম মান কেন নয়? তিনটি কারণ। এক, বর্গ মসৃণ: সব জায়গায় derivative আছে, তাই gradient descent সহজ, আর derivative শূন্য বসিয়ে সরাসরি সমাধানও পাওয়া যায় (অংশ ৩)। দুই, বর্গ বড় ভুলকে বেশি শাস্তি দেয়: ১০ মিনিটের একটি ভুল ১ মিনিটের দশটি ভুলের চেয়ে খারাপ (100 বনাম 10). ক্রেতার কাছে প্রায়ই সেটাই সত্য। তিন, গোলমাল যদি normal distribution মেনে চলে, তাহলে MSE কমানো আর maximum likelihood একই উত্তর দেয় (পর্ব ৩৮-এ MLE-এর যোগসূত্র)।

দ্বিতীয় কারণটির উল্টো পিঠও আছে। বড় ভুলকে বেশি শাস্তি দেওয়া মানে একটি অস্বাভাবিক বিন্দু, একটি outlier, রেখাকে নিজের দিকে অনেকখানি টেনে নিতে পারে। বিকল্প loss হলো MAE, ভুলের পরম মানের গড়; এটি কমালে যে রেখা আসে, তাকে least absolute deviations রেখাও বলে, আর এটি মধ্যমার (median) মতো আচরণ করে। দশটি কাল্পনিক অর্ডারে দুটিকে মেলানো হলো, তারপর একটি পার্সেলে ১২০ মিনিট যোগ করা হলো (ভুল ঠিকানায় দুই ঘণ্টা ঘুরেছে)। সত্যিকারের run, ডেস্কটপে, Python 3.14.4, scikit-learn 1.9.1 (MAE-রেখার জন্য QuantileRegressor(quantile=0.5, alpha=0)):

minutes: [21.1 16.9 27.8 30.9 36.1 41.6 43.  51.5 55.3 68.6]
clean data:
  MSE line: b =  10.71, w =  2.60
  MAE line: b =  10.36, w =  2.57
one outlier:
  MSE line: b = -13.29, w =  5.87
  MAE line: b =  10.36, w =  2.57

পরিষ্কার ডেটায় দুই রেখা প্রায় এক (আসল নিয়ম ছিল b = 12, w = 2.5). একটি বিন্দু বদলাতেই MSE-রেখার ঢাল দ্বিগুণের বেশি, 5.87, আর intercept ঋণাত্মক, −13.29: মানে ২ কিমির অর্ডারে model-এর অনুমান প্রায় −1.5 মিনিট! MAE-রেখা একটুও নড়েনি, কারণ শেষ বিন্দুটি এমনিতেই রেখার ওপরে ছিল, আর পরম মানের loss-এ “কত ওপরে” তার গুরুত্ব নেই, শুধু “ওপরে” না “নিচে”। তবে এর মূল্য আছে: MAE-এর কোনো সরল closed form নেই (linear programming লাগে), আর কোণাকার loss-এ gradient descent কম মসৃণ। মাঝামাঝি পথ হিসেবে Huber loss আছে (ছোট ভুলে বর্গ, বড় ভুলে পরম মান), scikit-learn-এ HuberRegressor.

৩ · Closed form — normal equations, হাতে তিনটি বিন্দু

MSE একটি বাটির মতো (convex) loss, পর্ব ৩৪-এর ছবির মতো; তলা একটিই, আর তলায় gradient শূন্য। পর্ব ৩৪-এর chain rule অনেক feature-এ খাটালে MSE-এর gradient হয়:

∇MSE(w) = (2/n) · Xᵀ (Xw − y)

Xᵀ হলো X-এর transpose (সারি আর কলাম অদলবদল)। এটি শূন্য বসালে Xᵀ(Xw − y) = 0, মানে:

XᵀX w = Xᵀy ⇒ w = (XᵀX)⁻¹ Xᵀy

এটিই normal equations, আর w-এর এই সরাসরি সূত্রই linear regression-এর closed form. সমীকরণটির একটি সুন্দর জ্যামিতিক অর্থ আছে: Xᵀr = 0 মানে residual vector r = y − Xw প্রতিটি কলামের সঙ্গে লম্ব (orthogonal, dot product শূন্য)। প্রথম কলাম সব 1, তাই একটি সঙ্গে সঙ্গে পাওয়া ফল: bias থাকলে residual-গুলোর যোগফল সবসময় শূন্য। ফুটবল মাঠে বলের ছায়া মাটিতে পড়ার মতো ভাবুন: Xw হলো y-এর “ছায়া” X-এর কলামগুলোর তৈরি সমতলে, আর r হলো বল থেকে ছায়া পর্যন্ত খাড়া দাগ। নাম “normal” এসেছে এই লম্ব (normal) হওয়া থেকে। এবার পরীক্ষার খাতার মতো তিনটি বিন্দুতে হাতে করি: দূরত্ব x = 2, 4, 6 কিমি, সময় y = 20, 26, 38 মিনিট।

1X = [[1,2],[1,4],[1,6]]2A = X.T @ X3c = X.T @ y4A_inv = inv(A)5b, w = A_inv @ c6r = y - X @ [b, w]

দুটি বিন্দু থাকলে কী হতো? তখন সমীকরণ দুটি, অজানা দুটি, রেখা দুই বিন্দু দিয়ে হুবহু যায়, residual শূন্য। যেমন x = 1, 3 আর y = 2, 6 হলে w = (6 − 2) ÷ (3 − 1) = 2, b = 0; normal equations-ও একই উত্তর দেবে। এই ছোট ঘটনাটি মনে রাখুন: parameter যত, বিন্দুও ততটি হলে model সব বিন্দু মুখস্থ করে ফেলে; অংশ ৯-এ polynomial-এর সঙ্গে এটি আবার আসবে।

৪ · কেন inverse নয়, lstsq — condition number আর তিনটি পথের তুলনা

কাগজে (XᵀX)⁻¹Xᵀy নিখুঁত। কম্পিউটারে এটি একটি পরিচিত ফাঁদ। পর্ব ১৭-এর এক লাইনের সারাংশ: condition number 10^k হলে float64-এর প্রায় ১৬টি দশমিক অঙ্কের মধ্যে k-টি হারানোর ঝুঁকি, আর সেই কারণে inv(A) @ b নয়, solve(A, b). এখানে একটি নতুন মোড় আছে: XᵀX বানানোর মুহূর্তেই condition number বর্গ হয়ে যায়, cond(XᵀX) ≈ cond(X)². X-এ ৭টি অঙ্ক হারানোর ঝুঁকি থাকলে XᵀX-এ ১৪টি।

কখন X খারাপ হয়? যখন দুটি কলাম প্রায় একই তথ্য দেয়। উদাহরণ: “সোজা দূরত্ব” আর “রাস্তার দূরত্ব”, যেখানে রাস্তা প্রায় হুবহু সোজা দূরত্বের 1.3 গুণ। ৫০টি কাল্পনিক অর্ডারে গোলমাল ছাড়া y বানানো হলো, যাতে সঠিক উত্তর ঠিক w = [10, 1, 1.5] জানা থাকে; তারপর তিনটি পথে সমাধান। সত্যিকারের run, ডেস্কটপে, Python 3.14.4, NumPy 2.5.4:

cond(X)     = 4.1e+07
cond(X.T X) = 1.9e+15
inv(X.T X) X.T y     w = [9.9065 0.9563 1.5352]  rel. error 1.1e-02
solve(X.T X, X.T y)  w = [10.      0.6559  1.7647]  rel. error 4.3e-02
lstsq(X, y)          w = [10.   1.   1.5]  rel. error 2.0e-11

cond(X) প্রায় 4 × 10⁷, আর XᵀX-এর প্রায় 2 × 10¹⁵, ঠিক বর্গের কাছাকাছি; float64-এর সীমার একদম কিনারায়। ফল: inverse-এর পথে আপেক্ষিক ভুল ১%, আর এখানে solve-ও বাঁচায়নি (৪%), কারণ সমস্যাটি inverse-এ নয়, XᵀX বানানোতেই; তথ্য সেখানেই হারিয়ে গেছে। lstsq XᵀX বানায়ই না; সে সরাসরি X-এর ওপর কাজ করে (NumPy-তে LAPACK-এর SVD-ভিত্তিক routine), তাই তার ভুল 2 × 10⁻¹¹, inverse-এর পথের চেয়ে প্রায় ৫০ কোটি গুণ কম। পর্ব ১৭-এর নিয়ম তাই এখানে আরও এক ধাপ এগোয়: least squares-এ inv তো নয়ই, normal equations নিজে না বানিয়ে lstsq ব্যবহার করুন। এই ভুলের অঙ্কগুলো একটি মেশিনের একটি run; অন্য BLAS বা অন্য CPU-তে শেষ অঙ্ক বদলাতে পারে, ধরন বদলাবে না।

ভালো আচরণের ডেটায় তিন পথই একই উত্তর দেয়। নিচের কোড অংশ ২-এর নিয়মে ২০০টি অর্ডার বানায়, তারপর normal equations (solve দিয়ে), lstsq আর scikit-learn-এর LinearRegression দিয়ে coefficient বের করে আসল নিয়মের পাশে রাখে। scikit-learn নিজেই intercept যোগ করে, তাই তাকে bias কলাম ছাড়া X দেওয়া হয়েছে।

import numpy as np
from sklearn.linear_model import LinearRegression

# বাজারঘরের ২০০টি কাল্পনিক অর্ডার, জানা "সত্যিকারের" নিয়ম দিয়ে বানানো
rng = np.random.default_rng(42)
n = 200
km = rng.uniform(1, 25, n)            # দূরত্ব, কিমি
kg = rng.uniform(0.2, 15, n)          # ওজন, কেজি
items = rng.integers(1, 9, n)         # পণ্যের সংখ্যা
rain = rng.binomial(1, 0.3, n)        # বৃষ্টি: 1 = হ্যাঁ
peak = rng.binomial(1, 0.35, n)       # ভিড়ের সময়: 1 = হ্যাঁ
true = np.array([12, 2.5, 0.8, 1.5, 9, 7])   # b, তারপর ৫টি w
X = np.column_stack([np.ones(n), km, kg, items, rain, peak])
y = X @ true + rng.normal(0, 4, n)    # মিনিট, সঙ্গে গোলমাল

w_ne = np.linalg.solve(X.T @ X, X.T @ y)       # normal equations
w_ls = np.linalg.lstsq(X, y)[0]                # least squares
sk = LinearRegression().fit(X[:, 1:], y)       # bias নিজে যোগ করে
w_sk = np.r_[sk.intercept_, sk.coef_]

names = ["bias", "km", "kg", "items", "rain", "peak"]
print("        আসল   normal  lstsq  sklearn")
for nm, t, a, b, c in zip(names, true, w_ne, w_ls, w_sk):
    print(f"{nm:6s} {t:6.2f} {a:7.3f} {b:6.3f} {c:7.3f}")
print("সর্বোচ্চ পার্থক্য (normal − lstsq):", f"{np.abs(w_ne - w_ls).max():.1e}")
r = y - X @ w_ls
print(f"RMSE = {np.sqrt(np.mean(r**2)):.2f} মিনিট, R² = {sk.score(X[:, 1:], y):.3f}")
print("একটি নতুন অর্ডার (10 km, 3 kg, 2 items, বৃষ্টি, ভিড় নয়):",
      round(float(np.array([1, 10, 3, 2, 1, 0]) @ w_ls), 1), "মিনিট")

        আসল   normal  lstsq  sklearn
bias    12.00  13.361 13.361  13.361
km       2.50   2.434  2.434   2.434
kg       0.80   0.775  0.775   0.775
items    1.50   1.413  1.413   1.413
rain     9.00   9.107  9.107   9.107
peak     7.00   6.870  6.870   6.870
সর্বোচ্চ পার্থক্য (normal − lstsq): 1.2e-13
RMSE = 3.78 মিনিট, R² = 0.957
একটি নতুন অর্ডার (10 km, 3 kg, 2 items, বৃষ্টি, ভিড় নয়): 52.0 মিনিট

৫ · অনেক feature-এ gradient descent — vectorised, scaling আর mini-batch

Closed form থাকলে gradient descent কেন? পর্ব ৩৪-এর উত্তর মনে করুন: বড় model-এ (logistic regression থেকে neural network) সরাসরি সূত্র থাকে না, আর লাখো feature-এ XᵀX-এর মতো (p × p) matrix-ই বিশাল। Linear regression হলো সেই জায়গা, যেখানে gradient descent-কে সঠিক উত্তরের সঙ্গে মিলিয়ে দেখা যায়। নিয়মটি এক লাইনের, অংশ ৩-এর gradient দিয়ে:

w ← w − η · (2/n) · Xᵀ(Xw − y)

η হলো learning rate. এক feature-এর জন্য পর্ব ৩৪-এ w আর b আলাদা দুই লাইনে বদলাতে হয়েছিল; এখানে একটি matrix গুণেই সব weight একসঙ্গে। অংশ ৩-এর তিনটি অর্ডারে, w = [0, 0] থেকে, η = 0.01 দিয়ে এক ধাপ:

1w = np.zeros(2)2y_hat = X @ w3r = y_hat - y4g = 2 / n * X.T @ r5w = w - eta * g6mse = np.mean((X @ w - y)**2)

শেষ ধাপের পর্যবেক্ষণটিই পর্ব ৩৪-এর “লম্বাটে বাটি”। Feature-গুলোর মাপ আলাদা হলে loss এক দিকে খুব খাড়া, অন্য দিকে প্রায় সমতল; খাড়া দিকের জন্য η ছোট রাখতে হয়, আর তাতে সমতল দিকে চলা কচ্ছপের গতিতে। প্রতিকার feature scaling, প্রচলিত রূপ standardization: প্রতিটি কলাম থেকে গড় বাদ দিয়ে standard deviation দিয়ে ভাগ। নিচের কোড অংশ ৪-এর একই ২০০টি অর্ডারে (একই seed, একই ক্রমে random সংখ্যা) দুইভাবে চালায়। দুই ক্ষেত্রেই η = 1 ÷ (loss-এর সবচেয়ে বড় বাঁক), যা নিরাপদ সীমার (2 ÷ বাঁক) ভেতরে, যাতে তুলনাটি ন্যায্য হয়। থামে gradient-এর দৈর্ঘ্য 10⁻⁶-এর নিচে নামলে।

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)                 # আগের মতোই একই ডেটা
n = 200
F = np.column_stack([rng.uniform(1, 25, n), rng.uniform(0.2, 15, n),
                     rng.integers(1, 9, n), rng.binomial(1, 0.3, n),
                     rng.binomial(1, 0.35, n)]).astype(float)
y = np.column_stack([np.ones(n), F]) @ [12, 2.5, 0.8, 1.5, 9, 7] + rng.normal(0, 4, n)

def gd(F, y, max_iter=100000, tol=1e-6):
    X = np.column_stack([np.ones(len(y)), F])
    H = 2 / len(y) * X.T @ X                     # loss-এর বাঁকের matrix
    eta = 1 / np.linalg.eigvalsh(H).max()       # নিরাপদ পায়ের মাপ
    w, hist = np.zeros(X.shape[1]), []
    for t in range(max_iter):
        r = X @ w - y                            # residual, shape (n,)
        g = 2 / len(y) * X.T @ r                 # vectorised gradient
        hist.append(np.mean(r ** 2))
        if np.linalg.norm(g) < tol:
            break
        w -= eta * g
    return w, t, hist, np.linalg.cond(H)

Z = (F - F.mean(axis=0)) / F.std(axis=0)         # standardize
for name, A in [("কাঁচা", F), ("standardized", Z)]:
    w, t, hist, c = gd(A, y)
    print(f"{name:12s}: cond = {c:7.1f}, ধাপ = {t:5d}, MSE = {hist[-1]:.3f}")
    plt.loglog(np.arange(1, len(hist) + 1), hist, label=f"{'raw' if name == 'কাঁচা' else 'standardized'} ({t} steps)")
print("standardized-এ w (প্রতি ১ SD-তে মিনিট):", w[1:].round(2))
print("মূল এককে ফেরত (প্রতি km, kg ...):", (w[1:] / F.std(axis=0)).round(3))
plt.xlabel("iteration"); plt.ylabel("training MSE")
plt.title("Same data, same rule: scaling changes the speed"); plt.legend()
plt.show()

কাঁচা       : cond =  3640.6, ধাপ = 51809, MSE = 14.314
standardized: cond =     1.2, ধাপ =    10, MSE = 14.314
standardized-এ w (প্রতি ১ SD-তে মিনিট): [16.45  3.3   3.28  4.33  3.23]
মূল এককে ফেরত (প্রতি km, kg ...): [2.434 0.775 1.413 9.107 6.87 ]

Mini-batch-এর জায়গা। ওপরের প্রতিটি ধাপে ২০০টি সারিই পড়া হয়েছে (full batch). কোটি সারিতে প্রতি পায়ে কয়েকটি সারির এক মুঠো নিয়ে gradient আন্দাজ করাই mini-batch, আর একটি করে সারি নিলে SGD; পর্ব ৩৪-এর অংশ ৯-এ ৪০০টি অর্ডারে তিন পদ্ধতির আসল তুলনা আছে। scikit-learn-এ এর তৈরি রূপ SGDRegressor, যেটিও scaling ছাড়া ভালো কাজ করে না। ছোট বা মাঝারি ডেটায় linear regression-এর জন্য lstsq-ই দ্রুত আর নিখুঁত।

৬ · Coefficient পড়া — “বাকি সব স্থির রেখে”, আর multicollinearity

অংশ ৪-এর km coefficient 2.434. এর সঠিক অর্থ: “বাকি সব feature (ওজন, পণ্য, বৃষ্টি, ভিড়) একই থাকলে, এক কিমি বেশি দূরত্বে model গড়ে প্রায় ২.৪ মিনিট বেশি অনুমান করে।” বৃষ্টির 9.107: “একই দূরত্ব, ওজন, পণ্য আর ভিড়ের অবস্থায়, বৃষ্টির দিনে প্রায় ৯ মিনিট বেশি।” তিনটি সূক্ষ্ম কথা এতে লুকানো আছে। এক, “স্থির রেখে” অংশটি বাদ দিলে অর্থ বদলে যায়। দুই, coefficient-এর একক আছে (মিনিট প্রতি কিমি, মিনিট প্রতি কেজি), তাই km-এর 2.4 আর kg-এর 0.78 সরাসরি তুলনা করে “দূরত্ব বেশি গুরুত্বপূর্ণ” বলা যায় না; তুলনা করতে হলে standardize করা coefficient দেখুন (অংশ ৫-এর 16.45, 3.3 …)। তিন, এটি পর্যবেক্ষণের ডেটায় একটি সম্পর্ক, কারণ নয়; বৃষ্টির দিনে হয়তো রাইডারও কম থাকে, আর model সেই প্রভাবও বৃষ্টির ঘাড়ে চাপায় (পর্ব ৩৭-এর confounder).

“স্থির রেখে” কথাটি তখনই অর্থহীন হয়, যখন দুটি feature একসঙ্গে ওঠানামা করে। অংশ ৪-এর “সোজা দূরত্ব” আর “রাস্তার দূরত্ব” ফিরিয়ে আনি: রাস্তার দূরত্ব বাস্তবে না বাড়িয়ে সোজা দূরত্ব বাড়ানো যায় না। এই অবস্থার নাম multicollinearity. ডেটা যেভাবে বানানো: আসল নিয়মে সময় নির্ভর করে শুধু রাস্তার দূরত্বে (প্রতি কিমিতে ২ মিনিট) আর বৃষ্টিতে (৮ মিনিট)। কিন্তু model-কে দুটি দূরত্বই দেওয়া হলো। নিচের কোড ১০টি ভিন্ন seed-এ ৬০টি করে অর্ডার বানিয়ে fit করে, তারপর একটি seed-এ statsmodels-এর OLS দিয়ে প্রতিটি coefficient-এর standard error আর ৯৫% confidence interval ছাপে (CI-এর অর্থ পর্ব ৩৭-এ)। statsmodels প্রথমবার লোড হতে কয়েক সেকেন্ড লাগে।

import numpy as np
import statsmodels.api as sm

def orders(seed, n=60):
    rng = np.random.default_rng(seed)
    km = rng.uniform(1, 25, n)                     # সোজা দূরত্ব
    road = 1.3 * km + rng.normal(0, 0.3, n)        # রাস্তার দূরত্ব, প্রায় একই তথ্য
    rain = rng.binomial(1, 0.3, n)
    y = 10 + 2.0 * road + 8 * rain + rng.normal(0, 4, n)
    return np.column_stack([km, road, rain]), y

print("seed  w_km w_road    মোট  road একা")     # মোট = w_km + 1.3·w_road
for seed in range(10):
    F, y = orders(seed)
    w = np.linalg.lstsq(sm.add_constant(F), y)[0]
    w1 = np.linalg.lstsq(sm.add_constant(F[:, 1:]), y)[0]
    print(f"{seed:4d} {w[1]:5.2f} {w[2]:6.2f} {w[1] + 1.3 * w[2]:6.2f} {w1[1]:6.2f}")

F, y = orders(0)
print("corr(km, road) =", round(np.corrcoef(F[:, 0], F[:, 1])[0, 1], 4))
res = sm.OLS(y, sm.add_constant(F)).fit()
print(res.summary(xname=["const", "km", "road", "rain"]).tables[1])
res1 = sm.OLS(y, sm.add_constant(F[:, 1:])).fit()
print("শুধু road রেখে 95% CI:\n", res1.conf_int().round(2))

seed  w_km w_road    মোট  road একা
   0  1.13   1.11   2.57   1.97
   1  2.92  -0.15   2.72   2.09
   2  1.41   1.09   2.83   2.18
   3 -1.00   2.77   2.60   2.01
   4 -2.07   3.56   2.56   1.95
   5  0.68   1.45   2.56   1.97
   6 -1.72   3.36   2.64   2.04
   7  1.34   1.02   2.67   2.05
   8 -3.33   4.51   2.54   1.94
   9 -1.30   3.07   2.68   2.07
corr(km, road) = 0.9995
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          9.9829      1.266      7.884      0.000       7.446      12.519
km             1.1286      2.525      0.447      0.657      -3.930       6.187
road           1.1068      1.943      0.570      0.571      -2.786       5.000
rain           7.4621      1.227      6.083      0.000       5.005       9.920
==============================================================================
শুধু road রেখে 95% CI:
 [[ 7.46 12.49]
 [ 1.85  2.1 ]
 [ 5.    9.87]]

সারকথা: multicollinearity অনুমানের (prediction) ক্ষতি সাধারণত করে না, কিন্তু coefficient-এর ব্যাখ্যা নষ্ট করে। ধরার উপায়: feature-এর correlation matrix, অংশ ৪-এর condition number, আর VIF (variance inflation factor), যা বলে অন্য feature-এর কারণে একটি coefficient-এর variance কত গুণ ফুলেছে। ওপরের seed 0-এ statsmodels-এর variance_inflation_factor ডেস্কটপে দিয়েছে km 937.1, road 937.6, rain 1.0 (সত্যিকারের run, statsmodels 0.15.0); মোটা নিয়মে ৫ বা ১০-এর বেশি হলে সাবধান হওয়ার সংকেত, যদিও এই সীমা সর্বজনীন নয়। প্রতিকার: একটি কলাম বাদ দেওয়া, দুটিকে মিলিয়ে একটি বানানো, বা regularization (পর্ব ৬১-এ ridge).

৭ · Category feature — one-hot encoding আর dummy-variable trap

বৃষ্টি আর ভিড় দুটি হ্যাঁ/না feature, তাই 0 আর 1 বসালেই চলেছে। কিন্তু “এলাকা” যদি ঢাকা, চট্টগ্রাম আর সিলেট হয়? ঢাকা = 0, চট্টগ্রাম = 1, সিলেট = 2 বসালে model ধরে নেবে চট্টগ্রাম ঠিক ঢাকা আর সিলেটের মাঝামাঝি, আর সিলেটের প্রভাব চট্টগ্রামের ঠিক দ্বিগুণ; এমন কোনো ক্রম বাস্তবে নেই। সমাধান one-hot encoding: প্রতিটি এলাকার জন্য একটি আলাদা 0/1 কলাম। Pandas-এ pd.get_dummies.

এখানেই একটি পরিচিত ফাঁদ। তিনটি এলাকার তিনটি কলামের যোগফল প্রতিটি সারিতে ঠিক 1, কারণ প্রতিটি অর্ডার কোনো একটি এলাকারই। আর bias কলামও সব 1. অর্থাৎ একটি কলাম অন্যগুলো থেকে হুবহু বের করা যায়: নিখুঁত multicollinearity. এর নাম dummy-variable trap. ৯০টি কাল্পনিক অর্ডারে (আসল নিয়ম: ঢাকায় বাড়তি 0, চট্টগ্রামে 6, সিলেটে 10 মিনিট) সত্যিকারের run, ডেস্কটপে, pandas 3.1.1, NumPy 2.5.4, scikit-learn 1.9.1:

['km', 'zone_Ctg', 'zone_Dhaka', 'zone_Sylhet']
columns: 5  rank: 4  cond: 1.0e+17
trap, sklearn: [20.02  2.52  0.81 -5.56  4.75]
trap, lstsq  : [15.01  2.52  5.81 -0.55  9.75]
['km', 'zone_Dhaka', 'zone_Sylhet']  rank: 4
b, km, Dhaka, Sylhet = [20.83  2.52 -6.36  3.94]

প্রথম তিন লাইন ফাঁদে পড়া অবস্থা: পাঁচটি কলাম, কিন্তু rank মাত্র 4 (একটি কলাম বাড়তি), আর cond প্রায় 10¹⁷, মানে XᵀX-এর inverse কার্যত নেই। মজার ব্যাপার, কোনো error আসেনি! scikit-learn আর lstsq দুজনেই উত্তর দিয়েছে, কিন্তু আলাদা উত্তর: একজন বলছে intercept 20.02, অন্যজন 15.01; একজনের ঢাকা −5.56, অন্যজনের −0.55. দুটি উত্তরই একই অনুমান দেয় (যেমন চট্টগ্রাম − ঢাকা: 0.81 − (−5.56) = 6.37, আর 5.81 − (−0.55) = 6.36), কারণ অসীম সংখ্যক উত্তর সমান ভালো, আর প্রতিটি library নিজের নিয়মে একটি বেছে নেয়। অনুমান ঠিক, ব্যাখ্যা অর্থহীন। শেষ দুই লাইনে প্রতিকার: drop_first=True একটি এলাকা (বর্ণক্রমে প্রথম, এখানে Ctg) বাদ দেয়; সে হয় “ভিত্তি” (baseline). এখন উত্তর অনন্য: intercept 20.83 (আসলে 15 + 6 = 21, চট্টগ্রামের ভিত্তি), ঢাকা −6.36 (আসলে চট্টগ্রামের চেয়ে 6 কম), সিলেট +3.94 (আসলে 4 বেশি)। প্রতিটি এলাকার coefficient পড়তে হবে ভিত্তি-এলাকার সঙ্গে তুলনা হিসেবে।

৮ · Residual analysis — ভুলের ছবি পড়া, আর একটি নতুন feature

ডাক্তার শুধু জ্বরের মাপ দেখে রোগ ধরেন না; কোন সময়ে জ্বর বাড়ে, কোন উপসর্গ সঙ্গে আসে, সেই নকশা দেখেন। RMSE বা R²-এর একটি সংখ্যাও তেমনি সব বলে না। Model কী ভুল করছে, তা দেখার সবচেয়ে কাজের ছবি হলো residual plot: x-অক্ষে model-এর অনুমান (fitted value), y-অক্ষে residual (আসল − অনুমান)। ভালো model-এ এই ছবি হওয়া উচিত শূন্য রেখার চারপাশে একটি নকশাহীন মেঘ। নকশা দেখা গেলেই, model ডেটার এমন কিছু ধরেনি যা ধরা সম্ভব। চারটি প্রচলিত নকশা ধাপে ধাপে দেখুন (ছকে আঁকা নমুনা, seed-করা বানানো সংখ্যা দিয়ে, কোনো আসল ডেটা নয়)।

এবার আসল run. নতুন কাল্পনিক ডেটা: শহরের জ্যামে দূরত্ব বাড়লে সময় বাড়ে দ্রুততর হারে, আসল নিয়ম y = 15 + 0.6·km + 0.09·km² + 8·rain + গোলমাল। প্রথমে সরলরেখা (km আর rain), তারপর km²-কে তৃতীয় feature হিসেবে যোগ করে। km² নতুন feature হলেও model এখনো linear, কারণ w-গুলোর ওপর রৈখিক; শুধু x-এর জায়গায় একটি নতুন কলাম। কোডে R², adjusted R² আর RMSE-ও হিসাব হয় (অর্থ পরের অংশে)।

import numpy as np
import matplotlib.pyplot as plt

# নতুন কাল্পনিক ডেটা: শহরের জ্যামে দূরত্ব বাড়লে সময় বাড়ে দ্রুততর হারে
rng = np.random.default_rng(7)
n = 120
km = rng.uniform(1, 25, n)
rain = rng.binomial(1, 0.3, n)
y = 15 + 0.6 * km + 0.09 * km**2 + 8 * rain + rng.normal(0, 3, n)

def fit(cols):
    X = np.column_stack([np.ones(n)] + cols)
    w = np.linalg.lstsq(X, y)[0]
    fitted = X @ w
    r = y - fitted
    p = X.shape[1] - 1                             # feature-এর সংখ্যা
    r2 = 1 - np.sum(r**2) / np.sum((y - y.mean())**2)
    adj = 1 - (1 - r2) * (n - 1) / (n - p - 1)
    print(f"{p} feature: w = {w.round(3)}")
    print(f"   R² = {r2:.3f}, adjusted R² = {adj:.3f}, RMSE = {np.sqrt(np.mean(r**2)):.2f} মিনিট")
    return fitted, r

f1, r1 = fit([km, rain])                           # সরলরেখা
f2, r2 = fit([km, rain, km**2])                    # km² যোগ করে
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4), sharey=True)
for ax, f, r, t in [(axes[0], f1, r1, "y ~ km + rain"),
                    (axes[1], f2, r2, "y ~ km + rain + km^2")]:
    ax.scatter(f, r, s=14, alpha=0.7)
    ax.axhline(0, color="k", lw=1)
    ax.set_title(t); ax.set_xlabel("fitted minutes")
axes[0].set_ylabel("residual (actual - fitted)")
plt.show()

2 feature: w = [4.073 2.948 7.603]
   R² = 0.950, adjusted R² = 0.949, RMSE = 4.80 মিনিট
3 feature: w = [14.521  0.671  7.667  0.088]
   R² = 0.983, adjusted R² = 0.982, RMSE = 2.83 মিনিট

বাঁক যদি উল্টো দিকের হতো (দূরের অর্ডারে হাইওয়ে, তাই প্রতি কিমিতে সময় কমে), U-এর বদলে উল্টো U দেখা যেত, আর log(km) বা √km বেশি মানানসই হতো। কোন রূপান্তর, তা ঠিক করে ডোমেইনের জ্ঞান আর residual-এর ছবি মিলিয়ে; অন্ধভাবে km³, km⁴ যোগ করে নয়, সেই বিপদ পরের অংশে।

৯ · R², RMSE, polynomial feature আর অনুমানগুলো

R², adjusted R² আর RMSE — সাদা বাংলায়

সবচেয়ে সরল baseline হলো সব অর্ডারে একই অনুমান, y-এর গড় (পর্ব ৩৪-এর অনুশীলনে দেখেছেন, MSE-কে গড়ই সবচেয়ে কম করে)। R² বলে, model সেই baseline-এর ভুলের কত অংশ দূর করেছে:

R² = 1 − SSE ÷ SST, SSE = ∑(y_i − ŷ_i)², SST = ∑(y_i − ȳ)²

ȳ হলো y-এর গড়। অংশ ৩-এর তিনটি অর্ডারে: ȳ = 28, SST = 64 + 4 + 100 = 168, SSE = 6, তাই R² = 1 − 6 ÷ 168 ≈ 0.964: গড়ের তুলনায় ভুলের প্রায় ৯৬% দূর হয়েছে। R² 1 মানে নিখুঁত, 0 মানে গড়ের চেয়ে ভালো নয়, আর নতুন ডেটায় মাপলে ঋণাত্মকও হতে পারে (গড়ের চেয়েও খারাপ)। একটি ফাঁদ: train ডেটায় নতুন feature যোগ করলে R² কখনো কমে না, এমনকি feature-টি সম্পূর্ণ এলোমেলো হলেও। Adjusted R² এর জন্য শাস্তি দেয়: 1 − (1 − R²)(n − 1) ÷ (n − p − 1), যেখানে p feature-এর সংখ্যা; অকেজো feature যোগ করলে এটি কমতে পারে। আর RMSE = √MSE সবচেয়ে সহজে বোঝানো যায়, কারণ এর একক মিনিট: “model সাধারণত প্রায় ৩–৪ মিনিট এদিক-ওদিক করে।” কোন metric কখন, train নয় test ডেটায় মাপা, আর MAE বনাম RMSE বাছাই পর্ব ৬০-এ (মূল্যায়ন); এখানে শুধু পরিচয়।

Polynomial feature — overfitting-এর সেতু

km² কাজ করেছে; তাহলে km³, km⁴ … km¹² যোগ করলে আরও ভালো? এগুলোকে বলে polynomial features. মাত্র ১৫টি train অর্ডারে (অংশ ৮-এর নিয়মে, বৃষ্টি বাদ) বিভিন্ন degree-র polynomial fit করে ২০০টি নতুন অর্ডারে যাচাই। সত্যিকারের run, ডেস্কটপে, NumPy 2.5.4 (numpy.polynomial.Polynomial.fit):

degree  1: train RMSE =  3.68, test RMSE =     5.47
degree  2: train RMSE =  2.61, test RMSE =     3.70
degree  4: train RMSE =  1.66, test RMSE =     4.60
degree  8: train RMSE =  1.56, test RMSE =     8.38
degree 12: train RMSE =  0.99, test RMSE =  2210.82

Train RMSE প্রতিবার কমেছে, ৩.৬৮ থেকে ০.৯৯; test RMSE সবচেয়ে কম degree 2-এ (৩.৭০, যা আসল নিয়মের degree), তারপর বেড়েছে, আর degree 12-এ ২,২১০ মিনিট! ১৩টি parameter আর ১৫টি বিন্দু: model প্রায় প্রতিটি train বিন্দু মুখস্থ করেছে, আর বিন্দুগুলোর ফাঁকে ও প্রান্তের বাইরে বুনো দোল খেয়েছে (কিছু test অর্ডারের দূরত্ব train-এর সীমার বাইরে, সেখানে দোল সবচেয়ে বেশি)। এটিই overfitting, আর ছোট ডেটায় এক run-এর সংখ্যা বলে seed বদলালে এগুলো বেশ বদলাবে। কেন এমন হয় (bias আর variance), learning curve আর regularization দিয়ে কীভাবে সামলানো যায়, তা পর্ব ৬১-এ (Overfitting, bias আর variance).

অনুমানগুলো — নিয়ম নয়, নির্দেশিকা

Linear regression-এর coefficient পেতে কোনো অনুমান লাগে না; lstsq যেকোনো ডেটায় সবচেয়ে কম বর্গের রেখা দেয়। অনুমান লাগে যখন সেই রেখা দিয়ে কিছু দাবি করবেন: standard error, CI, p-value, বা “এই coefficient-ই আসল প্রভাব”। পাঠ্যবইয়ে প্রচলিত তালিকা, মোটামুটি গুরুত্বের ক্রমে:

• Linearity: y-এর গড় feature-গুলোর রৈখিক যোগফল (দরকার হলে km²-এর মতো রূপান্তরিত feature-সহ)। না মানলে অংশ ৮-এর U.
• Independent errors: এক অর্ডারের ভুল আরেকটির ভুলের ওপর নির্ভর করে না। একই দিনের বা একই রাইডারের অর্ডার, বা সময়ের ধারায় পরপর দিন (পর্ব ২৬), প্রায়ই এটি ভাঙে; তখন standard error অতিরিক্ত আশাবাদী হয়।
• Constant variance: ভুলের ছড়ানো সব জায়গায় মোটামুটি এক; না মানলে ফানেল, আর CI ভুল।
• Normal errors (সবচেয়ে কম জরুরি): ছোট নমুনায় CI আর p-value-এর জন্য লাগে; বড় নমুনায় central limit theorem-এর কারণে সাধারণত কম গুরুত্বপূর্ণ।

এগুলো কঠোর পাস-ফেল পরীক্ষা নয়; বাস্তবের ডেটা কোনোটিই নিখুঁতভাবে মানে না। প্রশ্নটি হলো, ভাঙন কতটা, আর আপনার সিদ্ধান্তে তার প্রভাব কতটা। শুধু অনুমান (prediction) চাইলে অনেক ভাঙন সহনীয়, কারণ শেষ বিচার হয় test ডেটায়; coefficient-এর অর্থ বা CI দাবি করলে এগুলো গুরুত্বপূর্ণ। কোন পরীক্ষা কতটা জরুরি, তা নিয়ে পরিসংখ্যানবিদদের মধ্যেও মতভেদ আছে।

১০ · AI-এর যুগে — AI সহকারীর regression কোড আর ব্যাখ্যা যাচাই

“এই CSV-তে linear regression করে বলো কোন feature সবচেয়ে গুরুত্বপূর্ণ” বললে AI সহকারী কয়েক সেকেন্ডে কোড, টেবিল আর আত্মবিশ্বাসী ব্যাখ্যা দেয়। কোড প্রায়ই চলে; সমস্যা থাকে চোখের আড়ালে। যা প্রায়ই দেখা যায়:

• সূত্র হুবহু কোডে: np.linalg.inv(X.T @ X) @ X.T @ y, অংশ ৪-এর ফাঁদ; বা statsmodels-এ add_constant ভুলে যাওয়া, ফলে intercept ছাড়া model.
• পুরোনো বা অস্তিত্বহীন API: scikit-learn-এর LinearRegression(normalize=True) (১.২ version-এ সরানো হয়েছে), বা অস্তিত্বহীন কোনো parameter; ভুল error চাপা দিতে warnings.filterwarnings("ignore").
• Coefficient-কে কারণ বলা, আর একক না দেখে “বড় coefficient মানে বেশি গুরুত্বপূর্ণ”; multicollinear feature-এর ঋণাত্মক coefficient-এর গল্প বানানো।
• সব category-র dummy রেখে intercept-ও রাখা; scaling-এর গড় পুরো ডেটা থেকে নেওয়া (test-এর তথ্য train-এ ঢোকা); আর train ডেটাতেই R² মেপে “৯৮% নির্ভুল” দাবি।
• Residual plot না দেখা, আর না-চালানো আউটপুট বা p-value বানিয়ে লেখা।

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

• সমাধান lstsq বা scikit-learn দিয়ে, inv দিয়ে নয়; statsmodels-এ add_constant আছে
• জানা coefficient-এর বানানো ডেটায় চালিয়ে দেখেছি model সেগুলো মোটামুটি ফিরে পায়
• X-এর shape আর প্রথম কয়েকটি সারি দেখেছি; category-তে drop_first বা encoder-এর drop আছে, আর condition number বা VIF মেপেছি
• ডেটা আগে ভাগ করেছি; scaling শুধু train থেকে; RMSE আর R² test-এ, baseline-এর পাশে
• Residual plot দেখেছি; নকশা থাকলে feature বা রূপান্তর ভেবেছি, শুধু সংখ্যায় ভরসা করিনি
• Coefficient-এর ব্যাখ্যায় একক আর “বাকি সব স্থির রেখে” লিখেছি, কারণ দাবি করিনি; অচেনা API অফিসিয়াল docs-এ মিলিয়েছি, আর গ্রাহকের আসল ডেটা AI-তে দিইনি