Robust আর quantile regression — outlier আর অনুমানের পরিসর

AI-এর ভিত্তি কোর্সের পর্ব ৫৮-এ residual plot-এ “একা দূরের বিন্দু” দেখেছেন, আর পর্ব ৩৭-এ z-score ও IQR দিয়ে একটি column-এর outlier ধরেছেন। এই কোর্সের পর্ব ০২-এ Huber আর pinball loss-এর আকার দেখেছেন, আর দেখেছেন যে absolute loss median আর pinball loss একটি quantile শেখায়। পর্ব ০৭-এ least squares-কে projection আর hat matrix হিসেবে পড়েছেন। এই পর্ব সেগুলো regression-এ নিয়ে আসে: কোন বিন্দু রেখাকে সত্যিই টেনে নেয়, কীভাবে টান সীমিত রাখা যায়, আর গড়ের বদলে পুরো পরিসর কীভাবে অনুমান করা যায়। চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, pandas, statsmodels, scikit-learn আর matplotlib); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে, scikit-learn আর statsmodels-এ প্রায় দশ সেকেন্ড।

১ · একটি ভুল সারির জোর — outlier, leverage আর influence

বাজারঘরের ডেলিভারি দল প্রতিটি অর্ডারের দুটি সংখ্যা রাখে: hub থেকে দূরত্ব (km) আর ডেলিভারিতে কত ঘণ্টা লেগেছে। সাধারণ least squares, অর্থাৎ OLS, এর মধ্যে একটি সরলরেখা টানে। তিন ধরনের অস্বাভাবিক বিন্দু সেই রেখাকে তিনভাবে প্রভাবিত করে:

• Outlier: y-এর মান রেখা থেকে অনেক দূরে, যেমন মাঝারি দূরত্বের একটি অর্ডারে ৫০ ঘণ্টা, কারণ ক্রেতা দুদিন বাড়িতে ছিলেন না। x মাঝামাঝি হলে রেখা একটু ওপরে ওঠে, কিন্তু ঢাল তেমন বদলায় না।
• উঁচু leverage: x-এর মান বাকিদের থেকে অনেক দূরে, যেমন ১৫০ km দূরের একমাত্র অর্ডার। সেটি রেখার ধারা মেনে চললে ক্ষতি নেই, বরং ঢালের অনুমান আরও নিশ্চিত হয় (“ভালো leverage”).
• প্রভাবশালী বিন্দু (influence): উঁচু leverage আর রেখার ধারা ভাঙা, দুটি একসঙ্গে। এটিই রেখাকে নিজের দিকে টেনে নেয়, আর টেনে নেয় বলে তার নিজের residual ছোট দেখায়।

দুটি শব্দের পার্থক্য মনে রাখুন। Outlier একটি অবস্থা (মান অস্বাভাবিক), influence একটি প্রভাব (বাদ দিলে ফল বদলায়)। Outlier মানেই ভুল নয়: ঈদের আগের সপ্তাহে সত্যিই ডেলিভারি দেরি হয়। আবার ভুল মানেই outlier নয়: ১৫-কে ১৫০ লিখলে বিন্দুটি y-এ একদম স্বাভাবিক দেখায়। তাই প্রথম কাজ মাপা, তারপর কারণ খোঁজা, আর শেষে সিদ্ধান্ত।

২ · Hat value আর Cook's distance — হাতে হিসাব থেকে statsmodels

OLS-এর অনুমান ŷ = Hy, যেখানে H = X(XᵀX)⁻¹Xᵀ হলো hat matrix (y-এর মাথায় “টুপি” পরায়)। H-এর কর্ণের i-তম মান h_ii-কে বলে i-তম সারির hat value বা leverage: y_i এক একক বদলালে তার নিজের অনুমান ŷ_i কতটা বদলায়। একটি feature আর intercept থাকলে সূত্রটি সরল:

h_ii = ১/n + (x_i − x̄)² / Σ_j(x_j − x̄)²

প্রতিটি h_ii ০ থেকে ১-এর মধ্যে, আর সবগুলোর যোগফল p, অর্থাৎ parameter-এর সংখ্যা (H একটি projection, তার trace তার rank). তাই গড় leverage p/n, আর প্রচলিত একটি নিয়ম: h_ii > ২p/n হলে সারিটি দেখার মতো। h_ii শুধু x-এর ওপর নির্ভর করে, y-এর ওপর একটুও নয়। নিচের stepper পাঁচটি অর্ডারে হিসাবটি হাতে করে, তারপর Cook's distance পর্যন্ত যায়।

Cook's distance (R. Dennis Cook, ১৯৭৭) প্রশ্ন করে: i-তম সারি বাদ দিয়ে আবার fit করলে সব অনুমান মিলিয়ে কতটা সরে? আবার fit না চালিয়েই:

D_i = e_i² / (p · MSE) × h_ii / (1 − h_ii)²

প্রথম গুণনীয়ক residual কত বড়, দ্বিতীয়টি leverage কত বড়। দুটির যেকোনো একটি শূন্যের কাছে হলে D ছোট থাকে: এ কারণেই A বা B ধরনের বিন্দু D-তে উঁচু আসে না, C আসে। প্রচলিত সীমা D > ৪/n বা D > ১; এগুলো দেখার সংকেত, মোছার নিয়ম নয়। আরেকটি সহায়ক মাপ studentized residual: residual-কে তার নিজের standard error σ̂_(i)√(1 − h_ii) দিয়ে ভাগ করা, যেখানে σ̂_(i) সারিটি বাদ দিয়ে মাপা। উঁচু leverage-এর বিন্দুর সাধারণ residual ছোট দেখায়, কিন্তু √(1 − h) দিয়ে ভাগ করলে সে বড় হয়ে ধরা পড়ে।

এবার আসল আকারে: ৩০টি সাধারণ অর্ডার, সঙ্গে একটি অর্ডার যার দূরত্ব ভুল করে ১৫-এর জায়গায় ১৫০ লেখা হয়েছে; তার ঘণ্টা ১৫ km-এর অর্ডারের মতোই, ১০.২. সব ডেটা seed দিয়ে বানানো, কোনো আসল দোকানের নয়; আসল ঢাল ০.২৫ ঘণ্টা প্রতি km, যা model জানে না। চারটি ব্লকই নিচের module থেকে ডেটা নেয়।

# বাজারঘর (কাল্পনিক দোকান): ডেলিভারির সময় বনাম hub থেকে দূরত্ব, seed দিয়ে বানানো
import numpy as np

TRUE_SLOPE = 0.25                    # আসল ঢাল: প্রতি km-এ 0.25 ঘণ্টা (model জানে না)

def clean(n=30, seed=24):
    rng = np.random.default_rng(seed)
    km = rng.uniform(5, 60, n).round(1)
    hours = (6 + TRUE_SLOPE * km + rng.normal(0, 2, n)).round(1)
    return km, hours

def gross_errors(n=200, frac=0.10, seed=24):
    km, hours = clean(n, seed)       # তারপর 10% সারিতে বড় ভুল বসানো
    rng = np.random.default_rng(seed + 1)
    far = np.flatnonzero(km > 35)    # দূরের এলাকায় স্ট্যাটাস দেরিতে আপডেট
    bad = rng.choice(far, int(frac * n), replace=False)
    hours[bad] += rng.uniform(25, 50, bad.size).round(1)
    return km, hours, bad

def fan(n=400, seed=24):
    rng = np.random.default_rng(seed)
    km = rng.uniform(2, 80, n).round(1)
    noise = (rng.gamma(2.0, 1.0, n) - 2.0) / np.sqrt(2)   # গড় 0, SD 1, ডানে লেজ
    hours = 6 + TRUE_SLOPE * km + (0.5 + 0.08 * km) * noise
    return km, hours.round(2)

তিনটি function তিন ধরনের ডেটা দেয়: clean পরিষ্কার ডেটা, gross_errors দূরের এলাকার ১০% সারিতে ২৫ থেকে ৫০ ঘণ্টার বাড়তি ভুল (কুরিয়ার নেটওয়ার্ক না পেয়ে স্ট্যাটাস দেরিতে দিয়েছেন), আর fan এমন ডেটা যেখানে দূরত্ব বাড়লে সময়ের ছড়ানোও বাড়ে, ডানে লম্বা লেজসহ। প্রথম ব্লক statsmodels-এর OLSInfluence দিয়ে প্রতিটি সারির hat value, studentized residual আর Cook's distance বের করে। ছবির লেখা ইংরেজিতে, কারণ পাতার matplotlib বাংলা অক্ষর ঠিকমতো আঁকে না।

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import OLSInfluence
from bazar_delivery import clean

km, hours = clean(n=30)
km = np.append(km, 150.0)            # ভুল লেখা: 15 km-এর জায়গায় 150
hours = np.append(hours, 10.2)       # সময় ঠিক আছে, 15 km-এর মতোই
X = sm.add_constant(km)
fit = sm.OLS(hours, X).fit()
inf = OLSInfluence(fit)

n, p = X.shape
tab = pd.DataFrame({"km": km, "hours": hours, "hat": inf.hat_matrix_diag,
                    "stud_resid": inf.resid_studentized_external,
                    "cooks_d": inf.cooks_distance[0]})
print(tab.sort_values("cooks_d", ascending=False).head(4).round(3))
print(f"rules: hat > 2p/n = {2 * p / n:.3f}, Cook's D > 4/n = {4 / n:.3f}")

keep = np.arange(n) != tab.cooks_d.idxmax()   # শুধু সবচেয়ে প্রভাবশালী সারি বাদ
fit2 = sm.OLS(hours[keep], X[keep]).fit()
print("slope with the row    =", round(fit.params[1], 4))
print("slope without the row =", round(fit2.params[1], 4), "(true 0.25)")

grid = np.array([0, 160])
plt.figure(figsize=(6, 3.6)).subplots_adjust(bottom=0.16)
plt.scatter(km, hours, s=16, label="orders")
plt.scatter([150], [10.2], s=70, facecolors="none", edgecolors="r",
            label="km typed as 150")
plt.plot(grid, fit.params[0] + fit.params[1] * grid, "r--", label="OLS, all rows")
plt.plot(grid, fit2.params[0] + fit2.params[1] * grid, "g-", label="OLS, row dropped")
plt.xlabel("distance from hub (km)"); plt.ylabel("delivery time (hours)")
plt.legend(fontsize=8); plt.show()

       km  hours    hat  stud_resid  cooks_d
30  150.0   10.2  0.700      -7.553   22.669
24    8.0    3.7  0.071      -2.138    0.155
13    6.7    5.0  0.075      -1.737    0.114
12   56.9   20.2  0.055       1.358    0.053
rules: hat > 2p/n = 0.129, Cook's D > 4/n = 0.129
slope with the row    = 0.0597
slope without the row = 0.2561 (true 0.25)

৩ · Robust loss আবার — Huber আর M-estimator

ভুল সারি খুঁজে খুঁজে বাদ দেওয়ার বিকল্প হলো এমন loss নেওয়া, যা নিজেই বড় ভুলের টান সীমিত রাখে। পর্ব ০২-এ Huber loss দেখেছেন: residual r ছোট (∣r∣ ≤ δ) হলে r²/২, বড় হলে δ(∣r∣ − δ/২)। এখানে প্রশ্নটি অন্য: এই loss regression-এ কীভাবে কাজ করে, আর কেন কাজ করে?

যেকোনো loss ρ-এর জন্য Σρ(y_i − x_iᵀβ) কমানোর সমীকরণ হলো Σψ(r_i)x_i = ০, যেখানে ψ = ρ′ হলো প্রতিটি বিন্দুর রেখার ওপর “টান”। এমন অনুমানকে বলে M-estimator (maximum-likelihood ধাঁচের)। Squared loss-এ ψ(r) = r: টান residual-এর সমানুপাতিক, সীমাহীন। Huber-এ ψ(r) = r যতক্ষণ ∣r∣ ≤ δ, তারপর ±δ-তে থেমে যায়। সমীকরণটিকে Σw_ir_ix_i = ০ লেখা যায়, যেখানে ওজন w_i = ψ(r_i)/r_i = min(১, δ/∣r_i∣)। অর্থাৎ Huber regression আসলে একটি weighted least squares, যার ওজন residual দেখে বদলায়: পর্ব ২০-এর IRLS-এর মতোই বারবার ওজন হালনাগাদ করে চালানো যায়।

δ কীভাবে বাছবেন? residual-এর একক ঘণ্টা হলে δ-ও ঘণ্টায়, আর ডেটার scale বদলালে δ-ও বদলাতে হয়। তাই scikit-learn-এর HuberRegressor একটি scale σ-ও একসঙ্গে অনুমান করে, আর সীমা ধরে ∣r/σ∣ > ε. Default ε = ১.৩৫; Gaussian ভুলে এটি OLS-এর প্রায় ৯৫% efficiency রাখে, অর্থাৎ ভুল পরিষ্কার হলে সামান্যই দাম দিতে হয়। আরও দুটি খুঁটিনাটি: এর alpha default ০.০০০১, একটি ছোট L2 শাস্তি (পর্ব ১৭-এর ridge); আর এর loss-এর সূত্র পাঠ্যবইয়ের Huber-এর দ্বিগুণ (z² আর ২ε∣z∣ − ε²), যা সর্বনিম্ন বিন্দু বদলায় না।

Huber-এর একটি বড় সীমা আছে। ওজন residual দেখে কমে, কিন্তু অংশ ২-এর C বিন্দু রেখাকে টেনে নিজের residual ছোট করে ফেলে। তাই x-এ দূরের ভুলের বিরুদ্ধে Huber প্রায় OLS-এর মতোই দুর্বল। এই দুর্বলতার মাপের নাম breakdown point: ডেটার সবচেয়ে কম কত ভাগ ইচ্ছামতো নষ্ট করলে অনুমানকে যত খুশি দূরে সরানো যায়। গড় আর OLS-এর breakdown point ০ (একটি বিন্দুই যথেষ্ট), median-এর ৫০%. Huber-এর মতো M-estimator y-এর outlier-এ ভালো, কিন্তু leverage-সহ ভুলে তার breakdown point-ও ০. পরের দুটি পদ্ধতি অন্য পথে হাঁটে।

৪ · RANSAC আর Theil–Sen — ভোট আর median দিয়ে রেখা

RANSAC (RANdom SAmple Consensus; Fischler আর Bolles, ১৯৮১) এসেছে computer vision থেকে, যেখানে ছবির মিলের অর্ধেকও ভুল হতে পারে। কৌশলটি নির্বাচনের মতো: এলোমেলোভাবে ন্যূনতম কয়েকটি বিন্দু (সরলরেখায় ২টি) নিয়ে একটি রেখা বানান, তারপর গুনুন কতগুলো বিন্দু সেই রেখা থেকে একটি সীমার (threshold) ভেতরে আছে। এই কাজ বহুবার করে যে রেখার পক্ষে সবচেয়ে বেশি “ভোট” (inlier), তাকে রাখুন, আর শেষে শুধু সেই inlier-দের নিয়ে OLS চালান। একটি এলোমেলো জোড়া পরিষ্কার হওয়ার সম্ভাবনা যথেষ্ট থাকলে, বহু চেষ্টার মধ্যে অন্তত একটি পরিষ্কার জোড়া প্রায় নিশ্চিতভাবে আসে। দুর্বলতা: threshold আপনাকে বলে দিতে হয়, আর ফল এলোমেলো চেষ্টার ওপর নির্ভর করে, তাই random_state দিতে হয়। scikit-learn-এ threshold না দিলে default হলো y-এর median absolute deviation, residual-এর নয়: y-এর ছড়ানো বেশি হলে সীমাও ঢিলে হয়ে যায়।

Theil–Sen (Henri Theil ১৯৫০, Pranab Kumar Sen ১৯৬৮) আরও সরল: প্রতিটি জোড়া বিন্দুর ঢাল বের করুন, তারপর তাদের median নিন। একটি ভুল বিন্দু শুধু তার নিজের জোড়াগুলো নষ্ট করে, আর median সেগুলো এক পাশে ঠেলে দেয়। এর breakdown point প্রায় ২৯.৩%. নিচের stepper পাঁচটি বিন্দুতে পুরো হিসাব করে।

n টি বিন্দুতে জোড়া n(n − 1)/২টি, তাই বড় ডেটায় সব জোড়া গোনা ধীর। scikit-learn-এর TheilSenRegressor অনেক feature-এর জন্য সাধারণ রূপ চালায়: ছোট ছোট উপসেটে least squares, তারপর সব সমাধানের spatial median, আর উপসেটের সংখ্যা max_subpopulation (default ১০,০০০) ছাড়ালে এলোমেলো নমুনা। তাই এক feature-এও উত্তর ক্লাসিক “ঢালের median”-এর হুবহু সমান না-ও হতে পারে। এবার চারটি পদ্ধতি পাশাপাশি: ২০০টি অর্ডার, যার দূরের এলাকার ২০টিতে (১০%) ২৫ থেকে ৫০ ঘণ্টার বাড়তি ভুল।

import numpy as np
import pandas as pd
from sklearn.linear_model import (LinearRegression, HuberRegressor,
                                  RANSACRegressor, TheilSenRegressor)
from bazar_delivery import gross_errors, TRUE_SLOPE

km, hours, bad = gross_errors()          # 200 সারি, 20টিতে বড় ভুল
X = km.reshape(-1, 1)
print("rows:", len(km), "| planted errors:", bad.size)

models = {
    "OLS": LinearRegression(),
    "Huber": HuberRegressor(epsilon=1.35),
    "RANSAC": RANSACRegressor(residual_threshold=5.0, random_state=42),
    "TheilSen": TheilSenRegressor(random_state=42),
}
rows = []
for name, m in models.items():
    m.fit(X, hours)
    est = m.estimator_ if name == "RANSAC" else m
    slope, icpt = est.coef_[0], est.intercept_
    rows.append([name, icpt, slope, abs(slope - TRUE_SLOPE)])
res = pd.DataFrame(rows, columns=["model", "intercept", "slope", "slope_err"])
print(res.round(4).to_string(index=False))

hub = models["Huber"]
w_out = hub.outliers_                    # |r| > epsilon·sigma হলে True
print("Huber: scale sigma =", round(hub.scale_, 3),
      "| flagged =", int(w_out.sum()), "| planted among flagged =",
      int(np.isin(np.flatnonzero(w_out), bad).sum()))
ran = models["RANSAC"].inlier_mask_
print("RANSAC: inliers =", int(ran.sum()), "| planted kept as inliers =",
      int(ran[bad].sum()))

rows: 200 | planted errors: 20
   model  intercept  slope  slope_err
     OLS     2.1837 0.4828     0.2328
   Huber     5.9591 0.2599     0.0099
  RANSAC     6.6923 0.2269     0.0231
TheilSen     5.4706 0.2651     0.0151
Huber: scale sigma = 1.7 | flagged = 71 | planted among flagged = 20
RANSAC: inliers = 175 | planted kept as inliers = 0

৫ · গড় যথেষ্ট নয় — heteroscedasticity আর quantile regression

ঢাকার কাছের অর্ডার প্রায় সবসময় একই সময়ে পৌঁছায়; দূরের জেলায় কখনো ১২ ঘণ্টায়, কখনো দুদিনে, কারণ ফেরি, বৃষ্টি আর রাস্তা। ভুলের ছড়ানো x-এর সঙ্গে বদলানোকে বলে heteroscedasticity. OLS তখনও গড়ের একটি ভালো রেখা দিতে পারে, কিন্তু “গড় ± ধ্রুব প্রস্থ” দিয়ে বানানো পরিসর কাছে অযথা চওড়া আর দূরে বিপজ্জনকভাবে সরু। Quantile regression (Roger Koenker আর Gilbert Bassett, ১৯৭৮) প্রতিটি τ-এর জন্য আলাদা রেখা দেয়: Q_τ(y ∣ x) = xᵀβ_τ, অর্থাৎ দূরত্ব x হলে সময়ের τ-quantile. ছড়ানো বাড়লে ১০ আর ৯০তম quantile-এর রেখার ঢাল আলাদা হয়, রেখাগুলো পাখার মতো খোলে।

Loss হলো pinball loss: r = y − q̂ হলে r ≥ ০-এ τ·r, আর r < ০-এ (τ − ১)·r. কম বলা (আসল সময় অনুমানের চেয়ে বেশি) আর বেশি বলার দাম অসমান। নিচের stepper একটি residual দিয়ে শুরু করে দেখায়, এই অসমতা কেন ঠিক ৯০তম percentile-এ গিয়ে থামে।

Regression-এ q̂-এর জায়গায় xᵀβ বসে, আর Σ pinball(y_i − x_iᵀβ) কমাতে হয়। Loss টুকরো টুকরো রৈখিক আর convex, তাই সমস্যাটি একটি linear programming: প্রতিটি residual-কে দুটি অ-ঋণাত্মক অংশে (ওপরে u_i, নিচে v_i) ভাগ করে τΣu + (১ − τ)Σv কমানো, শর্ত y = Xβ + u − v. LP-র সমাধান একটি কোণে থাকে, তাই অন্তত p টি বিন্দু রেখার ঠিক ওপরে পড়ে। scikit-learn-এর QuantileRegressor এই LP সরাসরি SciPy-র HiGHS solver দিয়ে সমাধান করে (solver="highs", scikit-learn ১.৪ থেকে default). statsmodels-এর QuantReg iteratively reweighted least squares দিয়ে একই সর্বনিম্ন বিন্দুর কাছে পৌঁছায় আর standard error-ও দেয়। একটি ফাঁদ: QuantileRegressor-এর alpha default ১.০, একটি L1 শাস্তি (lasso-র মতো); শুধু quantile চাইলে alpha=0 দিতেই হবে।

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from sklearn.linear_model import QuantileRegressor
from bazar_delivery import fan

km, hours = fan(n=400)                   # দূরে গেলে ছড়ানো বাড়ে
X = sm.add_constant(km)
taus = [0.1, 0.5, 0.9]
fits = {t: sm.QuantReg(hours, X).fit(q=t) for t in taus}
for t in taus:
    b0, b1 = fits[t].params
    below = np.mean(hours <= b0 + b1 * km)
    print(f"tau={t}: intercept={b0:6.3f} slope={b1:.4f} share below={below:.3f}")
ols = sm.OLS(hours, X).fit().params
print(f"OLS mean: intercept={ols[0]:6.3f} slope={ols[1]:.4f}")

# যাচাই: scikit-learn, প্রথম 120টি সারিতে, LP solver "highs"
qr = QuantileRegressor(quantile=0.9, alpha=0.0, solver="highs")
qr.fit(km[:120].reshape(-1, 1), hours[:120])
sm90 = sm.QuantReg(hours[:120], X[:120]).fit(q=0.9).params
print("n=120, tau=0.9  sklearn:", np.round([qr.intercept_, qr.coef_[0]], 4),
      " statsmodels:", np.round(sm90, 4))

grid = np.linspace(0, 80, 2)
plt.figure(figsize=(6, 3.6)).subplots_adjust(bottom=0.16)
plt.scatter(km, hours, s=8, alpha=0.5, label="orders")
for t, c in zip(taus, ["tab:blue", "k", "tab:red"]):
    b0, b1 = fits[t].params
    plt.plot(grid, b0 + b1 * grid, color=c, label=f"quantile {t}")
plt.plot(grid, ols[0] + ols[1] * grid, "g--", label="OLS mean")
plt.xlabel("distance from hub (km)"); plt.ylabel("delivery time (hours)")
plt.legend(fontsize=8); plt.show()

tau=0.1: intercept= 5.490 slope=0.1655 share below=0.102
tau=0.5: intercept= 5.600 slope=0.2485 share below=0.500
tau=0.9: intercept= 6.814 slope=0.3586 share below=0.897
OLS mean: intercept= 5.558 slope=0.2682
n=120, tau=0.9  sklearn: [6.4908 0.388 ]  statsmodels: [6.4908 0.388 ]

৬ · Prediction interval, coverage আর quantile-এর কাটাকাটি

১০ আর ৯০তম quantile-এর রেখা মিলিয়ে একটি ৮০% prediction interval পাওয়া যায়: নতুন একটি অর্ডারের সময় প্রায় ৮০% ক্ষেত্রে এর ভেতরে পড়ার কথা। এটি গড়ের confidence interval নয়: CI বলে গড় রেখা কোথায় থাকতে পারে, prediction interval বলে একটি অর্ডারের সময় কোথায় পড়তে পারে, তাই এটি অনেক চওড়া। OLS থেকেও normal বণ্টন ধরে একটি ৮০% interval বানানো যায়, ŷ ± ১.২৮১৬·σ̂ (statsmodels-এর get_prediction দিয়ে এর পুরো রূপ পর্ব ২৩, Linear model গভীরে — basis function, interaction আর regularization path-এর বিষয়)। কিন্তু σ̂ সব দূরত্বে একই।

Interval-এর আসল পরীক্ষা coverage: training-এ কখনো না দেখা নতুন ডেটায় কত ভাগ মান interval-এর ভেতরে পড়ে। গড় coverage যথেষ্ট নয়; কাছে আর দূরে আলাদা করে দেখতে হয়, কারণ গড় ৮০% হতে পারে কাছে ৯৯% আর দূরে ৬০%-এর গড় হিসেবে। নিচের ব্লক আলাদা seed-এর ৪০০টি নতুন অর্ডারে দুটি interval মাপে, তারপর একটি দ্বিতীয় সমস্যা খোঁজে: quantile রেখার কাটাকাটি।

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.linear_model import QuantileRegressor
from bazar_delivery import fan

km, hours = fan(n=400, seed=24)              # train
km_t, hours_t = fan(n=400, seed=99)          # held-out, অন্য seed
X, Xt = sm.add_constant(km), sm.add_constant(km_t)
lo = Xt @ sm.QuantReg(hours, X).fit(q=0.1).params
hi = Xt @ sm.QuantReg(hours, X).fit(q=0.9).params
ols = sm.OLS(hours, X).fit()
mid, s = Xt @ ols.params, np.sqrt(ols.scale)
olo, ohi = mid - 1.2816 * s, mid + 1.2816 * s    # normal ধরে 80% interval
band = pd.cut(km_t, [0, 27, 54, 81], labels=["near", "mid", "far"])
inside = lambda a, b: (hours_t >= a) & (hours_t <= b)
cov = pd.DataFrame({"quantile_10_90": inside(lo, hi), "ols_normal_80": inside(olo, ohi)})
tab = cov.groupby(band, observed=True).mean()
tab.loc["all"] = cov.mean()
print(tab.round(3))
print("mean width  quantile:", round((hi - lo).mean(), 2), "| OLS:", round(2 * 1.2816 * s, 2))

rng = np.random.default_rng(42)             # ছোট নমুনা: 100 বার, প্রতিবার 25 সারি
grid = np.linspace(2, 80, 40).reshape(-1, 1)
crossed = 0
for _ in range(100):
    i = rng.choice(400, 25, replace=False)
    q = np.column_stack([QuantileRegressor(quantile=t, alpha=0, solver="highs")
                         .fit(km[i, None], hours[i]).predict(grid)
                         for t in (0.1, 0.5, 0.9)])
    crossed += bool(np.any(np.diff(q, axis=1) < 0))
print("samples with crossing inside 2-80 km:", crossed, "of 100")
q_fix = np.sort(q, axis=1)                   # rearrangement: প্রতি km-এ সাজানো
cross = lambda a: bool(np.any(np.diff(a, axis=1) < 0))
print("last sample, crossing before / after sorting:", cross(q), "/", cross(q_fix))

      quantile_10_90  ols_normal_80
near           0.818          0.991
mid            0.774          0.877
far            0.792          0.660
all            0.792          0.830
mean width  quantile: 9.65 | OLS: 10.89
samples with crossing inside 2-80 km: 48 of 100
last sample, crossing before / after sorting: True / False

এই কাটাকাটিকে বলে quantile crossing. কারণ সরল: প্রতিটি τ-এর রেখা আলাদাভাবে fit হয়, কেউ অন্যদের কথা জানে না। নমুনা ছোট হলে ঢালগুলো এলোমেলো হয়, আর দুটি সরলরেখার ঢাল আলাদা হলে কোথাও না কোথাও তারা কাটবেই; ডেটার পরিসরের বাইরে (extrapolation) এটি প্রায় নিশ্চিত। প্রতিকার তিনটি। (১) Rearrangement: প্রতিটি x-এ অনুমানগুলো সাজিয়ে নেওয়া; Chernozhukov, Fernández-Val আর Galichon ২০১০ সালে দেখান যে এতে অনুমান খারাপ হয় না, বরং সত্যের কাছে আসে। (২) Quantile-গুলো একসঙ্গে fit করা, ক্রমের শর্তসহ। (৩) Model-কে পরিসরের ভেতরে রাখা আর যথেষ্ট ডেটা নেওয়া।

৭ · কোন পরিস্থিতিতে কোন পদ্ধতি

কোনো একটি পদ্ধতি সবসময় সেরা নয়। প্রশ্নটি আগে ঠিক করুন: আপনি গড় চান, median চান, নাকি পুরো পরিসর? ডেটার সমস্যা কি y-এ, x-এ, নাকি ছড়ানোয়?

কয়েকটি কাজের নিয়ম। প্রথম, ভুল সারি চিনতে পারলে সংশোধনই সেরা প্রতিকার। দ্বিতীয়, OLS আর একটি robust পদ্ধতি পাশাপাশি চালান: দুটির উত্তর অনেক আলাদা হলে সেটিই সতর্কসংকেত, আর কারণ খুঁজুন। তৃতীয়, ছড়ানো বদলায় কি না residual বনাম fitted ছবিতে দেখুন; ফানেলের আকার পেলে interval-এর জন্য quantile regression বা log রূপান্তর ভাবুন, আর গড়ের standard error-এর জন্য robust SE (পর্ব ২০-এর sandwich SE). চতুর্থ, ব্যবসার খরচ অসমান হলে (দেরিতে পৌঁছানো আগে পৌঁছানোর চেয়ে খারাপ), গড় নয়, সেই খরচের অনুপাত থেকে τ বাছুন: দেরির দাম আগের দামের ৯ গুণ হলে τ = ৯/(৯ + ১) = ০.৯. পঞ্চম, Linear model-এর বাইরেও এই loss বাঁচে: gradient boosting-এ loss="quantile" বা loss="huber" দেওয়া যায়, যা পর্ব ৩৫ (Gradient boosting — residual থেকে functional gradient descent)-এ দেখবেন। আর outlier নিজেই যখন খোঁজার বিষয় (জালিয়াতি, যন্ত্রের ত্রুটি), তখন তা পর্ব ৪৩ (Anomaly detection — isolation forest, LOF আর one-class SVM)-এর কাজ।

৮ · AI-এর যুগে — robust আর quantile কোড যাচাই

“Outlier সরিয়ে model বানাও” বা “ডেলিভারির সময়ের range দাও” বললে AI assistant প্রায়ই চলনসই কোড দেয়, কিন্তু এই বিষয়ের ভুলগুলো চুপচাপ থাকে। নিচের ভুলগুলো বারবার দেখা যায়।

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

• ডেটা আগে ভাগ করেছি, test set তালা দেওয়া; outlier খোঁজা আর বাদ দেওয়া শুধু training অংশে
• প্রভাবশালী সারি মোছার আগে উৎস দেখেছি, কারণ লিখেছি, আর মোছার আগে-পরে ফল দুটোই রিপোর্টে রেখেছি
• OLS আর অন্তত একটি robust পদ্ধতি পাশাপাশি চালিয়েছি, কয়েকটি seed বা fold-এ
• scikit-learn-এর default (alpha, epsilon, residual_threshold, solver) নিজের version-এর docs-এ মিলিয়েছি
• ছোট একটি উদাহরণে hat value বা pinball loss হাতে হিসাব করে কোডের সঙ্গে মিলিয়েছি
• Interval-এর coverage held-out ডেটায়, x-এর ভাগে ভাগে মেপেছি, আর quantile crossing খুঁজেছি
• warning চাপা দিইনি, আর কোনো আসল গ্রাহকের ঠিকানা বা ফোন নম্বর AI-তে পেস্ট করিনি