ML দিয়ে forecasting — lag feature, backtesting আর অনেক series একসঙ্গে

AI-এর ভিত্তি কোর্সের পর্ব ২৬-এ pandas-এ তারিখ, resample, shift আর rolling শিখেছেন; পর্ব ৫৬-এ দেখেছেন সময়ের ডেটা এলোমেলো ভাগ করা কেন ভুল, আর TimeSeriesSplit-এর প্রথম পরিচয় পেয়েছেন; পর্ব ৬৩-এ lag feature এসেছিল একটি টীকা হিসেবে। আগের পর্ব ৫১ (Time series forecasting — decomposition, ARIMA আর exponential smoothing) একটি series-এর নিজের ছন্দ শেখার ক্লাসিক পথ দেখায়, সঙ্গে naive baseline আর MASE. এই পর্ব অন্য পথে হাঁটে: forecasting-কে সাধারণ regression বানিয়ে tree model দিয়ে শেখানো, আর তার সব ফাঁদ হাতে-কলমে ধরা। চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, pandas, scikit-learn আর matplotlib); প্রথমবার scikit-learn নামতে প্রায় দশ সেকেন্ড লাগে।

১ · সময়ের ধারা থেকে টেবিল — forecasting-কে supervised learning বানানো

Forecasting মানে অতীতের ধারা দেখে ভবিষ্যতের মান অনুমান। ক্লাসিক model (ARIMA, ETS) একটি series-এর নিজের গতিবিধি শেখে। ML পথে আমরা time series-টিকে একটি টেবিলে ভাঙি: প্রতিটি (তারিখ, জেলা) একটি সারি, লক্ষ্য সেই দিনের অর্ডার, আর কলামে এমন তথ্য, যা অনুমানের মুহূর্তে হাতে ছিল। এই মুহূর্তের নাম forecast origin, আর origin থেকে কত দিন সামনের অনুমান, তার নাম horizon h. টেবিল তৈরি হয়ে গেলে যেকোনো regression model চলে; এই পর্বে HistGradientBoosting, যা পর্ব ৩৬-এর LightGBM-ধাঁচের histogram boosting.

এই রূপান্তরের তিনটি লাভ। এক, একটি model একসঙ্গে অনেক series শিখতে পারে (অংশ ৫)। দুই, বাইরের তথ্য (ঈদ, দামের ছাড়, প্রচারণা) শুধু একটি কলাম যোগ করার ব্যাপার। তিন, tree model নিজে থেকেই অ-রৈখিক সম্পর্ক আর মিথস্ক্রিয়া ধরে, যেমন “ঈদের আগের শুক্রবার”। দামও তিনটি। Tree model প্রশিক্ষণের পরিসরের বাইরে যেতে পারে না: অর্ডার বছরে ২০% বাড়লে model নিজে সেই ঢাল আঁকে না, তাই স্তর ধরে রাখার কাজ lag-এর। সময়ের হিসাব আপনাকে নিজে রাখতে হয়, আর ভুল হলে leakage চুপচাপ ঢোকে (অংশ ২)। আর কয়েক দিন সামনের অনুমান কীভাবে হবে, তা আলাদা সিদ্ধান্ত (অংশ ৪)।

২ · Lag, rolling আর ক্যালেন্ডার feature — leakage ছাড়া

প্রতিটি কলামের জন্য একটিই প্রশ্ন: origin-এ দাঁড়িয়ে এই সংখ্যাটি কি সত্যিই জানা ছিল? এমন কলামকে বলে point-in-time feature. তিন ধরনের feature এই পরীক্ষায় পাস করে, যদি ঠিকভাবে বানানো হয়:

• Lag feature: k দিন আগের মান, pandas-এ groupby("district").shift(k). groupby বাদ দিলে এক জেলার প্রথম সারিতে আগের জেলার শেষ মান ঢুকে পড়ে (অনুশীলন ৩)। আগামীকালের অনুমানে (h = ১) k ≥ ১ নিরাপদ; h দিন সামনের অনুমানে k ≥ h হতে হবে।
• Rolling পরিসংখ্যান: আগের সাত দিনের গড়, যা দৈনিক এলোমেলো ওঠানামা মসৃণ করে স্তর দেখায়। নিয়ম: আগে shift, পরে rolling, উল্টো নয়। rolling(7).mean() সরাসরি লিখলে জানালার শেষ ঘর আজকের দিন।
• ক্যালেন্ডার feature: সপ্তাহের বার, ঈদ আসতে কত দিন, ঈদ গেছে কত দিন। তারিখ থেকেই হিসাব হয়, তাই ভবিষ্যতের দিনের জন্যও আগে জানা। এগুলো seasonality ধরে। আবহাওয়া এমন নয়: আগামীকালের বৃষ্টি আজ জানা নেই, তাই সেটি নিতে হলে আবহাওয়ার পূর্বাভাস নিতে হয়, আসল মাপ নয়।

এবার পুরো ডেটা। ছয় জেলার দুই বছরের দৈনিক অর্ডার নিচের module seed দিয়ে বানায়: জেলাভেদে আলাদা স্তর আর বৃদ্ধি, শুক্র-শনিবারে বেশি, ঈদের আগের দুই সপ্তাহে ভিড় আর ঈদের দিন ও পরের দুদিন কম, সঙ্গে কয়েক দিন টিকে থাকা এলোমেলো ধাক্কা। ঈদের তারিখগুলো আনুমানিক, আর প্রভাবের মাপ পুরোপুরি বানানো; কোনো আসল দোকানের বিক্রি নয়। চারটি ব্লকই এই module থেকে ডেটা নেয়।

# বাজারঘর (কাল্পনিক দোকান): ছয় জেলার দৈনিক অর্ডার, seed দিয়ে বানানো
import numpy as np
import pandas as pd

DISTRICTS = ["Dhaka", "Chattogram", "Sylhet", "Rajshahi", "Khulna", "Barishal"]
BASE = [220, 140, 70, 60, 55, 35]          # শুরুর দিকের গড় দৈনিক অর্ডার
# ঈদের আনুমানিক তারিখ; প্রভাবের আকার বানানো, কোনো আসল বিক্রির নয়
EID = pd.to_datetime(["2024-06-17", "2025-03-31", "2025-06-07",
                      "2026-03-20", "2026-05-27", "2027-03-10"])

def eid_features(dates):
    d = pd.DatetimeIndex(dates).values[:, None]
    diff = (EID.values[None, :] - d) / np.timedelta64(1, "D")
    ahead = np.where(diff >= 0, diff, 999).min(axis=1)   # ঈদ আসতে কত দিন
    since = np.where(diff <= 0, -diff, 999).min(axis=1)  # ঈদ গেছে কত দিন
    return np.minimum(ahead, 30), np.minimum(since, 10)

def make_orders(start="2024-07-01", days=730, seed=52):
    rng = np.random.default_rng(seed)
    dates = pd.date_range(start, periods=days, freq="D")
    t, dow = np.arange(days) / 365, dates.day_of_week.values
    ahead, since = eid_features(dates)
    eid = np.where(ahead <= 14, 1 + 0.06 * (14 - ahead), 1.0)  # ঈদের আগে ভিড়
    eid = np.where(since <= 2, 0.45, eid)       # ঈদের দিন আর পরের দুদিন কম
    parts = []
    for name, base in zip(DISTRICTS, BASE):
        fri, sat = rng.uniform(0.15, 0.40), rng.uniform(0.0, 0.15)
        week = 1 + fri * (dow == 4) + sat * (dow == 5) - 0.08 * (dow == 1)
        growth = 1 + rng.uniform(0.10, 0.35) * t
        mu = base * growth * week * eid * (1 + 0.08 * np.sin(2 * np.pi * t))
        shock = np.zeros(days)                  # টানা কয়েক দিন থাকা ধাক্কা (AR(1))
        for i in range(1, days):
            shock[i] = 0.8 * shock[i - 1] + rng.normal(0, 0.09)
        y = rng.poisson(mu * np.exp(shock))
        parts.append(pd.DataFrame({"date": dates, "district": name, "orders": y}))
    return pd.concat(parts, ignore_index=True)

FEATS = ["lag1", "lag7", "lag14", "roll7", "dow", "eid_ahead", "eid_since"]

def add_features(df):
    df = df.sort_values(["district", "date"]).reset_index(drop=True)
    g = df.groupby("district")["orders"]
    for k in (1, 7, 14):
        df[f"lag{k}"] = g.shift(k)
    df["roll7"] = g.transform(lambda s: s.shift(1).rolling(7).mean())  # আগে shift
    df["dow"] = df["date"].dt.day_of_week
    df["eid_ahead"], df["eid_since"] = eid_features(df["date"])
    df["district_id"] = df["district"].map(DISTRICTS.index)
    return df.dropna().reset_index(drop=True)

প্রথম ব্লক feature বানায়, ১৫ মে-র সারিটি ছাপে (স্টেপারের সঙ্গে মিলিয়ে দেখুন), তারপর একজন বিশ্লেষকের সাধারণ ভুল নকল করে: “trend” নামে একটি rolling(7).mean() কলাম, shift ছাড়া। শেষে শেষ আট সপ্তাহের test-এ চারটি অনুমানের MAE: seasonal naive (গত সপ্তাহের একই বার), naive (গতকাল), সৎ HistGradientBoosting, আর trend-যুক্ত HistGradientBoosting.

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor
from bazar_orders import make_orders, add_features, FEATS

raw = make_orders()
df = add_features(raw)
print(raw.shape, "->", df.shape, "| rows with a missing lag dropped")
g = df.groupby("district")["orders"]
df["trend7"] = g.transform(lambda s: s.rolling(7).mean())   # shift ভুলে গেছি!
one = df[(df.district == "Barishal") & (df.date == "2026-05-15")]
print(one[["orders"] + FEATS + ["trend7"]].round(2).to_string(index=False))
back = 7 * df.trend7 - 7 * df.roll7 + df.lag7               # leak থেকে আজকের অর্ডার
print("max |rebuilt - orders| =", round((back - df.orders).abs().max(), 6))

cut = df.date.max() - pd.Timedelta(days=56)                 # শেষ ৮ সপ্তাহ test
tr, te = df[df.date <= cut], df[df.date > cut]
mae = lambda p: np.mean(np.abs(te.orders.to_numpy() - p))
print(f"seasonal naive (lag7)   MAE = {mae(te.lag7.to_numpy()):5.2f}")
print(f"naive (lag1)            MAE = {mae(te.lag1.to_numpy()):5.2f}")
for name, cols in [("HGB honest", FEATS), ("HGB + trend7 (leak)", FEATS + ["trend7"])]:
    m = HistGradientBoostingRegressor(random_state=0).fit(tr[cols], tr.orders)
    print(f"{name:23s} MAE = {mae(m.predict(te[cols])):5.2f}")

(4380, 3) -> (4296, 11) | rows with a missing lag dropped
 orders  lag1  lag7  lag14  roll7  dow  eid_ahead  eid_since  trend7
     75  60.0  78.0   61.0  57.14    4       12.0       10.0   56.71
max |rebuilt - orders| = 0.0
seasonal naive (lag7)   MAE = 44.48
naive (lag1)            MAE = 27.86
HGB honest              MAE = 19.35
HGB + trend7 (leak)     MAE = 15.59

৩ · Backtesting — rolling origin, expanding আর sliding window

শেষ আট সপ্তাহের একটি time-based split সৎ, কিন্তু একটি মাত্র নমুনা। সেই আট সপ্তাহে ঈদ পড়লে এক রকম ফল, না পড়লে আরেক রকম। Backtesting বা rolling-origin evaluation তাই অতীতের কয়েকটি তারিখে দাঁড়ায়: প্রতিটি origin-এ শুধু তার আগের ডেটায় শেখায়, পরের কয়েক সপ্তাহে মাপে, তারপর origin সামনে সরায়। এটি সময়ের ডেটার cross-validation (AI-এর ভিত্তি কোর্সের পর্ব ৬৬-এর KFold-এর বদলে)। Train অংশ কেমন হবে, তার দুটি ধরন আছে: expanding window, যা শুরু থেকে origin পর্যন্ত সব নেয়, আর sliding window, যা শুধু শেষ নির্দিষ্ট কয়েক দিন নেয়।

দ্বিতীয় ব্লক এই ভাগ নিজে লেখে, gap = ০ আর ৭ দুই ক্ষেত্রেই scikit-learn-এর TimeSeriesSplit-এর সঙ্গে মেলায়, তারপর পাঁচটি fold-এ seasonal naive আর HistGradientBoosting তুলনা করে। একটি খুঁটিনাটি গুরুত্বপূর্ণ: ছয়টি জেলার সারি একটির নিচে আরেকটি সাজানো, তাই ভাগ হয় তারিখের তালিকায়, সারির নম্বরে নয়। সারির নম্বরে TimeSeriesSplit চালালে “train” পায় ঢাকা আর চট্টগ্রাম, “test” পায় বরিশাল, আর সময়ের ক্রম ভেঙে যায়।

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import TimeSeriesSplit
from sklearn.ensemble import HistGradientBoostingRegressor
from bazar_orders import make_orders, add_features, FEATS
df = add_features(make_orders())
days = np.sort(df.date.unique())              # ভাগ হবে তারিখে, সারিতে নয়
def rolling_origin(n, n_splits, test_size, gap=0):
    for k in range(n_splits, 0, -1):          # প্রথম fold সবচেয়ে পুরোনো
        end = n - (k - 1) * test_size
        yield np.arange(0, end - test_size - gap), np.arange(end - test_size, end)
for gap in (0, 7):
    mine = list(rolling_origin(len(days), 5, 28, gap))
    skl = list(TimeSeriesSplit(n_splits=5, test_size=28, gap=gap).split(days))
    same = all((a == c).all() and (b == d).all() for (a, b), (c, d) in zip(mine, skl))
    print(f"gap={gap}: identical to TimeSeriesSplit: {same} | last fold: "
          f"train 0..{mine[-1][0][-1]}, test {mine[-1][1][0]}..{mine[-1][1][-1]}")
rows, folds = [], list(rolling_origin(len(days), 5, 28))   # h = 1, তাই gap = 0
for k, (tr_i, te_i) in enumerate(folds, 1):
    tr, te = df[df.date.isin(days[tr_i])], df[df.date.isin(days[te_i])]
    m = HistGradientBoostingRegressor(random_state=0).fit(tr[FEATS], tr.orders)
    err = lambda p: np.mean(np.abs(te.orders.to_numpy() - p))
    rows.append([k, str(days[te_i[0]])[:10], err(te.lag7.to_numpy()),
                 err(m.predict(te[FEATS]))])
res = pd.DataFrame(rows, columns=["fold", "test_from", "snaive_MAE", "HGB_MAE"])
print(res.round(2).to_string(index=False))
print("mean MAE:", res[["snaive_MAE", "HGB_MAE"]].mean().round(2).to_dict())
fig, ax = plt.subplots(figsize=(6, 2.6)); fig.subplots_adjust(bottom=0.2)
for k, (tr_i, te_i) in enumerate(folds, 1):
    ax.barh(k, len(tr_i), color="#93c5fd"); ax.barh(k, 28, left=te_i[0], color="#f97316")
ax.set_xlabel("day index (train = blue, test = orange)"); ax.set_ylabel("fold")
ax.invert_yaxis(); plt.show()

gap=0: identical to TimeSeriesSplit: True | last fold: train 0..687, test 688..715
gap=7: identical to TimeSeriesSplit: True | last fold: train 0..680, test 688..715
 fold  test_from  snaive_MAE  HGB_MAE
    1 2026-02-11       22.04    13.68
    2 2026-03-11       59.16    21.19
    3 2026-04-08       22.19    14.88
    4 2026-05-06       62.43    22.27
    5 2026-06-03       26.52    16.24
mean MAE: {'snaive_MAE': 38.47, 'HGB_MAE': 17.65}

Gap কখন লাগে? এখানে h = ১, আর প্রতিটি সারি লক্ষ্য-তারিখ দিয়ে চিহ্নিত, তাই gap = ০ সঠিক। কিন্তু ধরুন সারি চিহ্নিত origin দিয়ে, আর লক্ষ্য h দিন পরের অর্ডার (বা “সামনের সাত দিনের মোট”)। তখন train-এর শেষ h − ১টি সারির লক্ষ্য test-এর সময়ের ভেতরে পড়ে: origin-এ সেগুলো এখনো ঘটেনি। gap = h − 1 সেই সারিগুলো ফেলে দেয়। ডেটা দেরিতে চূড়ান্ত হলেও (ফেরত আর বাতিল দুদিন পরে হিসাবে আসে) একই কাজ করতে হয়।

একই backtest-এ ক্লাসিক model — ডেস্কটপে সত্যিকারের run

পর্ব ৫১-এর ক্লাসিক model কি এই ছয় জেলায় আরও ভালো? statsmodels-এ প্রতি জেলায় প্রতিটি fold-এ SARIMA fit পাতায় ধীর, তাই একটি ডেস্কটপ প্রকল্পে Nixtla-র statsforecast চালানো হয়েছে, একই module আর একই পাঁচটি fold দিয়ে: প্রতি ২৮ দিনে নতুন করে fit (refit=28), আর মাঝের দিনগুলোতে আগের দিনের আসল মান দেখে পরের দিনের অনুমান, অর্থাৎ HistGradientBoosting-এর মতোই h = ১.

$ uv init --vcs none --author-from none --python 3.14 bazarghar-forecast
Initialized project `bazarghar-forecast` at `~/bazarghar-forecast`
$ cd bazarghar-forecast
$ uv add statsforecast
…
 + statsforecast==2.1.1
 + statsmodels==0.15.0
…
$ uv pip list | grep -E "^(numpy|pandas|statsforecast) "
numpy             2.5.4
pandas            2.3.3
statsforecast     2.1.1
$ uv run python sf_backtest.py

import sys
import time
import numpy as np
import pandas as pd
import statsforecast
from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA, AutoETS
from bazar_orders import make_orders

raw = make_orders()
data = raw.rename(columns={"district": "unique_id", "date": "ds", "orders": "y"})
sf = StatsForecast(models=[AutoETS(season_length=7),
                           AutoARIMA(season_length=7)], freq="D")
t0 = time.perf_counter()
cv = sf.cross_validation(df=data, h=1, step_size=1, n_windows=140, refit=28)
secs = time.perf_counter() - t0
cv = cv.sort_values(["unique_id", "ds"]).reset_index(drop=True)
cv["fold"] = (cv.groupby("unique_id").cumcount() // 28) + 1   # 28 দিনের ৫টি fold
models = ["AutoETS", "AutoARIMA"]
mae = cv.groupby("fold")[models].apply(
    lambda g: (g.sub(cv.loc[g.index, "y"], axis=0)).abs().mean())
print(mae.astype(float).round(2).to_string())
print("first test day per fold:", cv.groupby("fold").ds.min().dt.date.astype(str).tolist())
print("mean:", mae.astype(float).mean().round(2).to_dict())
print(f"rows: {len(cv)} | fit + forecast: {secs:.1f} s")
print("python", sys.version.split()[0], "| statsforecast", statsforecast.__version__,
      "| pandas", pd.__version__, "| numpy", np.__version__)

      AutoETS  AutoARIMA
fold                    
1       13.44      15.65
2       28.98      28.14
3       15.01      16.42
4       30.45      33.12
5       14.07      19.51
first test day per fold: ['2026-02-11', '2026-03-11', '2026-04-08', '2026-05-06', '2026-06-03']
mean: {'AutoETS': 20.39, 'AutoARIMA': 22.57}
rows: 840 | fit + forecast: 124.6 s
python 3.14.4 | statsforecast 2.1.1 | pandas 2.3.3 | numpy 2.5.4

৪ · কয়েক দিন সামনে — recursive বনাম direct

ঈদের আগে hub-এ মাল যায় এক বা দুই সপ্তাহ আগে, তাই আগামীকালের অনুমান যথেষ্ট নয়; দরকার ১ থেকে ১৪ দিন সামনের অনুমান। সমস্যা: lag1 দিয়ে শেখা model পরশুর অনুমানে আগামীকালের lag1 চায়, যা origin-এ অজানা। দুটি প্রচলিত পথ আছে।

• Recursive (iterated): একটি h = ১ model. আগামীকাল অনুমান করে সেই সংখ্যাকে আসল মানের জায়গায় বসান, তারপর পরের দিন, এভাবে ১৪ বার। Model একটি, কিন্তু প্রতিটি ধাপের ভুল পরের ধাপের feature-এ ঢোকে।
• Direct: প্রতিটি h-এর জন্য আলাদা model, যার feature শুধু origin পর্যন্ত জানা: h = ১০-এর model-এর সবচেয়ে কাছের lag হলো lag10, আর “একই বার”-এর lag হলো lag14. নিজের অনুমান ফেরত খায় না, কিন্তু ১৪টি model, আর দূরের h-এর model সাম্প্রতিক তথ্য কম পায়।

তৃতীয় ব্লক দুটি পথই চালায়। একটি origin পর্যন্ত (শেষ ১৫৫ দিনের আগে) শেখানো হয়, তারপর আর refit নয়; পরের ১৪১টি দিনের প্রতিটিকে origin ধরে ১ থেকে ১৪ দিন সামনের অনুমান, ছয় জেলায়। Feature-এর function X(P, t, h) একটিই: h দিলে শুধু t − h পর্যন্ত জানা মান নেয়। Recursive পথে P-তে origin-এর পরের সব মান মুছে দিয়ে নিজের অনুমান বসানো হয়; direct পথে প্রতিটি h-এর আলাদা model আসল ইতিহাস থেকে সরাসরি পড়ে। HistGradientBoosting প্রথমবার ১৪টি model শিখতে কয়েক সেকেন্ড নিতে পারে।

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import HistGradientBoostingRegressor as HGB
from bazar_orders import make_orders, eid_features
raw = make_orders()
Y = raw.pivot(index="date", columns="district", values="orders").to_numpy(float)
days = np.sort(raw.date.unique()); H, n_d = 14, Y.shape[1]
cal = [pd.DatetimeIndex(days).day_of_week.to_numpy(), *eid_features(days)]
def X(P, t, h):    # target দিন t-এর feature, শুধু t − h পর্যন্ত জানা তথ্য দিয়ে
    w, j = 7 * int(np.ceil(h / 7)), np.arange(len(t))   # w: একই বারের জানা দিন
    cols = [P[j, t - h], P[j, t - w], P[j, t - w - 7],
            np.mean([P[j, t - h - k] for k in range(7)], axis=0)]
    cols += [np.tile(v[t][:, None], (1, n_d)) for v in cal]
    return np.stack(cols + [np.tile(np.arange(n_d), (len(t), 1))], -1).reshape(-1, 8)
T0, t_tr = len(Y) - 155, np.arange(28, len(Y) - 154)   # T0-র পরে আর refit নয়
YY = np.broadcast_to(Y, (len(t_tr),) + Y.shape)
direct = {h: HGB(random_state=0, categorical_features=[7]).fit(X(YY, t_tr, h),
          Y[t_tr].ravel()) for h in range(1, H + 1)}       # প্রতি h-এ একটি model
o = np.arange(T0, len(Y) - H); J = np.arange(len(o))       # প্রতিদিন একটি origin
P = np.repeat(Y[None], len(o), 0); err = np.zeros((3, H))
P[np.arange(len(Y))[None, :] > o[:, None]] = np.nan       # origin-এর পরের সব মুছে দেওয়া
for h in range(1, H + 1):
    t, truth = o + h, Y[o + h].ravel()
    P[J, t] = direct[1].predict(X(P, t, 1)).reshape(-1, n_d)   # নিজের অনুমান ফেরত
    d = direct[h].predict(X(np.broadcast_to(Y, P.shape), t, h))
    s = Y[t - 7 * int(np.ceil(h / 7))].ravel()                  # seasonal naive
    err[:, h - 1] = [np.abs(truth - p).mean() for p in (P[J, t].ravel(), d, s)]
tab = pd.DataFrame(err.T, index=pd.Index(range(1, H + 1), name="h"),
                   columns=["recursive", "direct", "snaive"])
print(len(o), "origins x", n_d, "districts\n", tab.iloc[[0, 1, 3, 6, 9, 13]].round(1))
plt.figure(figsize=(6, 3)).subplots_adjust(bottom=0.17)
plt.plot(tab.index, tab.to_numpy(), marker="o", label=list(tab.columns))
plt.xlabel("horizon h (days ahead)"); plt.ylabel("MAE (orders)"); plt.legend(); plt.show()

141 origins x 6 districts
     recursive  direct  snaive
h                            
1        17.3    17.3    38.4
2        20.3    20.0    38.4
4        22.7    22.6    38.5
7        25.3    23.5    38.4
10       26.6    25.1    42.9
14       27.4    26.0    42.8

৫ · অনেক series একসঙ্গে — global বনাম local model

বাজারঘরের ছয়টি জেলা, আর বাস্তবে হাজারো পণ্য। দুটি পথ। Local model: প্রতিটি series-এর আলাদা model, শুধু নিজের ইতিহাস থেকে; ক্লাসিক ARIMA আর ETS এভাবেই চলে। Global model: সব series-এর সারি একটি টেবিলে, series-এর পরিচয় একটি কলাম, আর একটি model. Global model-এর সুবিধা: শুক্রবারের উল্লম্ফন বা ঈদের ভিড়ের মতো সবার মধ্যে মিল থাকা ধরন সে সব জেলার ডেটা থেকে একসঙ্গে শেখে, আর ছোট বা নতুন series বড়দের থেকে ধার নেয়। ঝুঁকি: বড় series (ঢাকা) loss-এ বেশি জায়গা নেয়, আর আকারে খুব আলাদা series-এ নিজের স্তর আলাদা রাখতে হয়। এখানে district_id একটি category কলাম, যা HistGradientBoosting-এর native categorical support দিয়ে split-এ জেলাগুলোকে দলে ভাগ করতে পারে; অনেকে প্রতিটি series-কে তার সাম্প্রতিক গড় দিয়ে ভাগ করে আকারও সমান করেন।

চতুর্থ ব্লকে একটি বাস্তব পরিস্থিতি: বরিশাল hub খুলেছে মাত্র ৯০ দিন আগে, তাই তার ইতিহাস ছোট। একই শেষ আট সপ্তাহে global আর local model তুলনা, তারপর global model দিয়ে ১০, ৫০ আর ৯০তম quantile-এর অনুমান (অংশ ৬)।

import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor as HGB
from bazar_orders import make_orders, add_features, FEATS

df = add_features(make_orders())
cut = df.date.max() - pd.Timedelta(days=56)              # শেষ ৮ সপ্তাহ test
new = (df.district == "Barishal") & (df.date <= cut - pd.Timedelta(days=90))
df = df[~new]                                            # বরিশাল hub নতুন: ৯০ দিনের ইতিহাস
tr, te = df[df.date <= cut], df[df.date > cut]
G = FEATS + ["district_id"]

glob = HGB(categorical_features=["district_id"], random_state=0).fit(tr[G], tr.orders)
te = te.assign(global_=glob.predict(te[G]), local=np.nan)
for d, part in tr.groupby("district"):                   # জেলাপ্রতি আলাদা model
    m = HGB(random_state=0).fit(part[FEATS], part.orders)
    te.loc[te.district == d, "local"] = m.predict(te.loc[te.district == d, FEATS])
ae = te[["global_", "local"]].sub(te.orders, axis=0).abs()
tab = ae.groupby(te.district, sort=False).mean()
tab["train_rows"] = tr.groupby("district", sort=False).size()
print(tab.round(1)); print("all:", ae.mean().round(2).to_dict())

q = {a: HGB(loss="quantile", quantile=a, categorical_features=["district_id"],
            random_state=0).fit(tr[G], tr.orders).predict(te[G]) for a in (0.1, 0.5, 0.9)}
inside = (te.orders >= q[0.1]) & (te.orders <= q[0.9])
print("coverage of 10-90 band:", round(inside.mean(), 3),
      "| below q0.1:", round((te.orders < q[0.1]).mean(), 3),
      "| above q0.9:", round((te.orders > q[0.9]).mean(), 3))
print(inside.groupby(te.district, sort=False).mean().round(2).to_dict())
row = te.district.eq("Dhaka").to_numpy().nonzero()[0][-1]
print("Dhaka, last day: q10/q50/q90 =", [round(q[a][row]) for a in q], "actual =", te.orders.iloc[row])

            global_  local  train_rows
district                              
Barishal        7.5   10.5          90
Chattogram     30.6   37.2         660
Dhaka          43.5   44.0         660
Khulna         10.6   10.6         660
Rajshahi       11.8   14.5         660
Sylhet         13.4   13.6         660
all: {'global_': 19.56, 'local': 21.73}
coverage of 10-90 band: 0.729 | below q0.1: 0.098 | above q0.9: 0.173
{'Barishal': 0.77, 'Chattogram': 0.68, 'Dhaka': 0.66, 'Khulna': 0.77, 'Rajshahi': 0.7, 'Sylhet': 0.8}
Dhaka, last day: q10/q50/q90 = [261, 307, 331] actual = 328

স্তরের অনুমান মেলানো। বাজারঘরের series একটি স্তরবিন্যাসে সাজানো: দেশ → বিভাগ → জেলা → পণ্য। প্রতিটি স্তরে আলাদা model বানালে সংখ্যাগুলো মেলে না: ছয় জেলার অনুমানের যোগফল ৭৬০, অথচ দেশের model বলছে ৮০০. Hierarchical forecasting-এ এই অমিল দূর করার কাজকে বলে reconciliation. Bottom-up: শুধু নিচের স্তর অনুমান করে যোগ করা (৭৬০); সহজ, কিন্তু নিচের series শব্দ-ভরা। Top-down: দেশের ৮০০-কে ইতিহাসের অনুপাতে জেলায় ভাগ; মসৃণ, কিন্তু জেলার নিজস্ব ধরন হারায়। আরও আধুনিক একটি পথ সব স্তরের অনুমান একসঙ্গে নিয়ে ভুলের covariance দিয়ে ওজন করে মেলায়, যার পরিচিত রূপ MinT (Wickramasuriya, Athanasopoulos আর Hyndman, ২০১৯)।

৬ · Quantile forecast আর মজুদের সিদ্ধান্ত

Hub-এর ম্যানেজার একটি সংখ্যা চান না, চান “কতটা রাখলে বেশিরভাগ দিন কম পড়বে না”। পর্ব ২৪-এ দেখেছেন, pinball loss একটি quantile শেখায়; HistGradientBoosting-এ loss="quantile" আর quantile=τ দিলে সেই loss-এ gradient boosting চলে। এভাবে পাওয়া quantile forecast-এর ১০ আর ৯০তম মান মিলিয়ে একটি ৮০% prediction interval.

ব্লক ৪-এর শেষ অংশ এই interval-কে শেষ আট সপ্তাহে যাচাই করে। Band ধরেছে ৭২.৯%, লক্ষ্য ৮০%. নিচে পড়েছে ৯.৮% (লক্ষ্য ১০%, ঠিক আছে), কিন্তু ওপরে বেরিয়েছে ১৭.৩% (লক্ষ্য ১০%): ওপরের লেজ ছোট করে ধরা। কারণ দুটি। অর্ডার বছরে বাড়ছে, আর tree model প্রশিক্ষণের সর্বোচ্চ মানের ওপরে যেতে পারে না; আর test-এ ঈদের ভিড় আছে। জেলাভেদে coverage ০.৬৬ (ঢাকা) থেকে ০.৮০ (সিলেট)। শিক্ষা: interval-এর নাম “৮০%” হলেই সে ৮০% ধরে না; backtest-এ মাপুন, জেলা ধরে ধরে। প্রতিকার: backtest-এর ভুল দেখে band চওড়া করা, লক্ষ্যকে স্তর দিয়ে ভাগ করে আকারহীন করা, বা conformal পদ্ধতি, যা পর্ব ৫৮ (Conformal prediction — নিশ্চয়তাসহ অনুমানের পরিসর)-এ আসবে।

কোন quantile মজুদ করবেন? এটি একটি সিদ্ধান্তের প্রশ্ন (AI-এর ভিত্তি কোর্সের পর্ব ৫৪-এর expected cost-এর হিসাব)। এক ইউনিট কম রাখলে হারানো লাভ আর ক্রেতার বিরক্তি মিলিয়ে দাম c_u, বেশি রাখলে গুদাম আর ফেরত পাঠানোর দাম c_o. ক্লাসিক newsvendor problem-এর উত্তর: চাহিদার τ = c_u/(c_u + c_o) quantile পরিমাণ রাখুন। c_u = ৳৮০ আর c_o = ৳২০ হলে τ = ৮০/১০০ = ০.৮; অর্থাৎ ৮০তম percentile-এর model দরকার, গড়ের নয়। ব্লকের শেষ লাইনে ঢাকার শেষ দিনে q১০/q৫০/q৯০ = ২৬১/৩০৭/৩৩১, আসল ৩২৮: median অনুযায়ী মাল রাখলে সেদিন ২১টি অর্ডার কম পড়ত।

৭ · Forecasting-এর library আর প্রতিযোগিতার শিক্ষা

এই পর্বের সব কিছু pandas আর scikit-learn দিয়ে নিজে লেখা, যাতে প্রতিটি ধাপ চোখে দেখা যায়। কাজের জায়গায় কয়েকটি library এই হিসাব গুছিয়ে দেয়। নিচের version লেখার সময়ে PyPI-তে দেখা; API বদলায়, তাই নিজের version-এর docs দেখুন।

Prophet নিয়ে একটি সতর্কতা: এটি সহজ বলে জনপ্রিয়, কিন্তু কয়েকটি প্রকাশিত তুলনায় দেখা গেছে, অনেক ডেটায় এটি seasonal naive বা ETS-এর মতো সরল পদ্ধতিকে নিয়মিত হারায় না। যেকোনো library-র “auto” model-কে নিজের backtest-এ seasonal naive-এর পাশে মেপে তবেই বিশ্বাস করুন।

প্রতিযোগিতা থেকে। Spyros Makridakis-এর M প্রতিযোগিতাগুলো কয়েক দশক ধরে forecasting পদ্ধতির খোলা পরীক্ষা। ২০১৮ সালের M4-এ ছিল এক লাখ series; জয়ী পদ্ধতি ছিল Slawek Smyl-এর exponential smoothing আর recurrent neural network-এর সংকর, আর কয়েকটি পদ্ধতির মিশ্রণ (combination) ভালো করেছিল, যেখানে একা ML পদ্ধতিগুলো সাধারণত পিছিয়ে ছিল বলে আয়োজকেরা লিখেছিলেন। ২০২০ সালের M5-এ একটি বড় খুচরা বিক্রেতার হাজারো পণ্যের দৈনিক বিক্রি, স্তরবিন্যাসসহ; শীর্ষ সমাধানগুলোর বেশিরভাগ LightGBM-ভিত্তিক global model, lag আর ক্যালেন্ডার feature দিয়ে, ঠিক এই পর্বের পথে। M5-এর একটি আলাদা অংশে quantile অনুমানও চাওয়া হয়েছিল। দুটি ফলাফল একসঙ্গে পড়ুন: অনেক সম্পর্কিত series আর বাইরের তথ্য থাকলে global ML model শক্তিশালী; একটি-দুটি ছোট series-এ ক্লাসিক model আর মিশ্রণ এখনো কঠিন প্রতিদ্বন্দ্বী।

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

“আগামী মাসের বিক্রির forecast model বানাও” বললে AI assistant দ্রুত চলনসই কোড দেয়, আর তার score প্রায়ই সুন্দর। সমস্যা হলো, সময়ের ডেটার ভুলগুলো score-কে খারাপ করে না, ভালো করে। নিচের ভুলগুলো বারবার দেখা যায়।

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

• একটি তারিখ বেছে তার প্রতিটি feature হাতে হিসাব করেছি; কোনোটিতে origin-এর পরের তথ্য নেই
• ভাগ তারিখে করেছি, সারিতে নয়; shuffle কোথাও নেই; শেষ অংশ চূড়ান্ত test হিসেবে তালা দেওয়া
• সব fold-এর MAE seasonal naive-এর পাশে দেখেছি, শুধু গড় বা শুধু শেষ fold নয়
• Multi-step অনুমানে h ধরে ভুল ছেপেছি, আর recursive backtest-এ আসল lag ঢোকেনি তা যাচাই করেছি
• অবিশ্বাস্য রকম ভালো score পেলে আগে leakage খুঁজেছি, একটি করে feature সরিয়ে
• Interval-এর coverage backtest-এ, series ধরে মেপেছি; scikit-learn আর pandas-এর parameter নিজের version-এর docs-এ মিলিয়েছি
• কোনো আসল ক্রেতা বা দোকানের বিক্রির ডেটা অনুমতি ছাড়া AI-তে পেস্ট করিনি