Exponential family আর GLM — Poisson, binomial আর link function

AI-এর ভিত্তি কোর্সের পর্ব ৩৬-এ Bernoulli, binomial আর Poisson বণ্টন দেখেছেন, পর্ব ৫৮-এ linear regression আর পর্ব ৫৯-এ logistic regression. এই কোর্সের পর্ব ১১-এ logistic-এর Hessian XᵀSX বের করেছিলেন, আর পর্ব ১৫-এ শুনেছিলেন যে পরিসংখ্যানের library এই Newton-এর ধাপকে IRLS নামে চালায়। এই পর্ব সেই সুতোগুলো এক জায়গায় বাঁধে: linear, logistic আর Poisson regression আসলে একটিই model-এর তিন রূপ, যার নাম generalized linear model. Likelihood আর MLE-এর মূলনীতি পর্ব ১৮ (Maximum likelihood আর MAP — model fit করার মূলনীতি)-এর বিষয়; এখানে শুধু যতটুকু লাগে ততটুকু। চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, SciPy, statsmodels আর scikit-learn); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে, scikit-learn আর statsmodels-এ প্রায় দশ সেকেন্ড।

১ · গোনা সংখ্যা আর হ্যাঁ/না — সাধারণ regression কোথায় আটকায়

Linear regression ধরে নেয় y = Xβ + ভুল, যেখানে ভুল সব জায়গায় একই মাপের normal বণ্টন। গোনা সংখ্যার (count) লক্ষ্যে তিনটি অনুমানই ভাঙে: অনুমান ঋণাত্মক হতে পারে, ভুলের ছড়ানো গড়ের সঙ্গে বাড়ে (heteroscedasticity), আর feature-এর প্রভাব যোগ নয়, গুণ। হ্যাঁ/না লক্ষ্যে অবস্থা আরও খারাপ: সম্ভাবনা ০ আর ১-এর বাইরে চলে যায়। সমাধান তিনটি আলাদা প্রশ্নে ভাগ করা: y কোন বণ্টন থেকে আসে? feature-গুলো কীভাবে একত্র হয়? সেই যোগফল গড়ের সঙ্গে কীভাবে জোড়া? এই ভাগের নামই generalized linear model (GLM).

টেবিলের প্রথম তিনটি বণ্টন একটি পরিবারের সদস্য। সেই পরিবারটি আগে চিনে নিলে বাকি সব যন্ত্র, অর্থাৎ fit করার পদ্ধতি, residual আর যাচাই, একবারেই সবার জন্য লেখা যায়।

২ · Exponential family — এক ছাঁচে Bernoulli, Poisson আর Gaussian

একটি বণ্টন exponential family-র সদস্য, যদি তার সম্ভাবনা (বা density) এই রূপে লেখা যায়:

p(y | θ) = h(y) · exp( θ · T(y) − A(θ) )

চারটি অংশ। θ হলো natural parameter: বণ্টনের এমন একটি parameter, যা y-এর সঙ্গে শুধু গুণ হয়ে সূচকে বসে। T(y) হলো sufficient statistic: ডেটার যে সারাংশ জানলেই θ সম্পর্কে ডেটার সব তথ্য পাওয়া যায়; এই পর্বের তিনটি বণ্টনেই T(y) = y. A(θ) হলো log-partition function: সম্ভাবনার যোগফল ১ রাখার জন্য যে log-সংশোধন লাগে। আর h(y) θ-মুক্ত একটি ভিত্তি অংশ। Poisson-কে এই রূপে আনার কাজটি ধাপে ধাপে দেখুন।

A′(θ) = গড় কেন হয়? সম্ভাবনার যোগফল ১, তাই e^A(θ) = Σ h(y)e^θy. দুই দিকে θ-এর সাপেক্ষে derivative নিলে A′(θ)e^A(θ) = Σ y·h(y)e^θy, অর্থাৎ A′(θ) = Σ y·p(y) = E[y]. আরেকবার derivative নিলে A″(θ) = Var(y). Variance কখনো ঋণাত্মক নয়, তাই A″ ≥ ০: A(θ) একটি convex function. এই ছোট তথ্যটিই পরে নিশ্চিত করবে যে canonical link-এর GLM-এ log-likelihood-এর একটিই চূড়া।

শেষ কলামের V(μ)-কে বলে variance function: গড় μ হলে variance কত। Gaussian-এ σ² আলাদা একটি সংখ্যা, যাকে বলে dispersion parameter φ; সাধারণ রূপে Var(y) = φ·V(μ). Bernoulli আর Poisson-এ φ = ১ ধরা থাকে, কারণ variance পুরোপুরি গড় দিয়েই ঠিক হয়। নিচের ব্লক তিনটি বণ্টনকে ছাঁচের রূপে লিখে SciPy-র নিজস্ব logpmf/logpdf-এর সঙ্গে মেলায়, আর central difference দিয়ে A′ ও A″ মেপে গড় ও variance-এর সঙ্গে তুলনা করে।

import numpy as np
from scipy import stats
from scipy.special import gammaln, expit

# p(x) = h(x)·exp(θ·T(x) − A(θ)); তিনটিতেই T(x) = x
fams = {
    "Bernoulli p=0.3": (np.log(0.3 / 0.7),            # θ = log-odds
                        lambda t: np.log1p(np.exp(t)),  # A(θ)
                        lambda x: 0 * x,                # log h(x)
                        np.array([0, 1]), stats.bernoulli(0.3).logpmf),
    "Poisson lam=4":   (np.log(4.0), np.exp,
                        lambda x: -gammaln(x + 1),
                        np.arange(0, 12), stats.poisson(4).logpmf),
    "Gaussian mu=2":   (2.0, lambda t: t ** 2 / 2,      # σ = 1 জানা
                        lambda x: -x ** 2 / 2 - 0.5 * np.log(2 * np.pi),
                        np.linspace(-2, 6, 9), stats.norm(2, 1).logpdf),
}
h = 1e-4
for name, (t, A, logh, x, ref) in fams.items():
    logp = logh(x) + t * x - A(t)                       # exponential-family রূপ
    gap = np.max(np.abs(logp - ref(x)))                 # scipy-র সঙ্গে পার্থক্য
    dA = (A(t + h) - A(t - h)) / (2 * h)                # A′(θ)
    d2A = (A(t + h) - 2 * A(t) + A(t - h)) / h ** 2     # A″(θ)
    print(f"{name:16s} theta={t:+.4f} A={A(t):.4f} max|gap|={gap:.0e}")
    print(f"{'':16s} A'(theta)={dA:.4f}  A''(theta)={d2A:.4f}")

print("Bernoulli: expit(theta) =", round(expit(np.log(0.3 / 0.7)), 4))
print("Poisson:   exp(theta)   =", round(np.exp(np.log(4.0)), 4))

Bernoulli p=0.3  theta=-0.8473 A=0.3567 max|gap|=2e-16
                 A'(theta)=0.3000  A''(theta)=0.2100
Poisson lam=4    theta=+1.3863 A=4.0000 max|gap|=0e+00
                 A'(theta)=4.0000  A''(theta)=4.0000
Gaussian mu=2    theta=+2.0000 A=2.0000 max|gap|=2e-15
                 A'(theta)=2.0000  A''(theta)=1.0000
Bernoulli: expit(theta) = 0.3
Poisson:   exp(theta)   = 4.0

৩ · GLM-এর তিন অংশ — random component, linear predictor আর link

John Nelder আর Robert Wedderburn ১৯৭২ সালে একটি প্রবন্ধে দেখান, linear, logistic, probit, Poisson আর আরও কয়েকটি regression একই কাঠামোয় লেখা যায় আর একই algorithm দিয়ে fit করা যায়। কাঠামোর তিনটি অংশ:

Random component ঠিক করে ভুলের আকৃতি, linear predictor ঠিক করে feature-গুলো কীভাবে মেশে, আর link function দুটিকে জোড়ে। Poisson-এ log link নিলে μ = e^η সবসময় ধনাত্মক, আর η-তে যোগ মানে μ-তে গুণ: চায়ের দোকানের “২০% কম” ঠিক এভাবেই আসে। নিচের ছবিতে বাজারঘরের একটি জেলার একটি সারি পুরো পথ পেরোয়।

Canonical link আর logistic regression-এর পরিচয়

যে link η-কে সরাসরি natural parameter θ বানায়, অর্থাৎ g(μ) = θ, তাকে বলে canonical link. অংশ ২-এর টেবিল থেকে পড়ে নিন: Gaussian-এর canonical link identity, Poisson-এর log, আর Bernoulli-র logit, log(μ/(1 − μ)). তাই logistic regression হলো Bernoulli random component + linear predictor + logit link-এর GLM; AI-এর ভিত্তি কোর্সের পর্ব ৫৯-এর sigmoid আসলে logit-এর উল্টো। জেলার ১০০টি অর্ডারের মধ্যে কয়টি ফেরত এল, এমন ডেটায় binomial GLM একই logit link নেয়, শুধু প্রতিটি সারির ওজন তার অর্ডারের সংখ্যা।

Canonical link-এর দুটি সুবিধা আছে। প্রথমত, log-likelihood β-তে concave, তাই চূড়া একটিই (ডেটা যথেষ্ট আর X-এর column-গুলো স্বাধীন হলে)। দ্বিতীয়ত, MLE-এর সমীকরণ খুব সরল হয়: Xᵀ(y − μ) = ০. Intercept থাকলে এর প্রথম সারি বলে Σy_i = Σμ̂_i: fit হওয়া model মোট অর্ডার হুবহু মেলায়। তবে canonical link বাধ্যতামূলক নয়। Binomial-এ probit বা complementary log-log, Gamma-তে canonical inverse link-এর বদলে log link খুব প্রচলিত, কারণ ব্যাখ্যা সহজ আর μ ধনাত্মক থাকার নিশ্চয়তা আছে। Link একটি model-অনুমান; ব্যাখ্যা আর যাচাইয়ের ফল দেখে বাছতে হয়।

৪ · Poisson regression আর exposure — জেলাভিত্তিক অর্ডারের model

এই পর্বের ডেটা: বাজারঘরের ৬৪টি জেলার এক মাসের অর্ডার, সব seed দিয়ে বানানো, কোনো আসল দোকানের নয়। প্রতিটি সারিতে মোট অর্ডার, মাসে কত দিন ডেলিভারি খোলা ছিল (নতুন শাখা বা বন্যার কারণে ৫ থেকে ৩০ দিন), সেই মাসে ছাড়ের ক্যাম্পেইন চলেছে কি না, আর ঢাকার hub থেকে দূরত্ব (১০০ km এককে)। ডেটা বানানোর আসল সূত্রে প্রতি দিনের হার e^{১.২ + ০.৪·promo − ০.৬·dist}; model সেটি জানে না, ডেটা থেকে খুঁজে বের করবে। চারটি ব্লকই নিচের module থেকে ডেটা নেয়।

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

def make_orders(seed=20, nb_alpha=None):
    rng = np.random.default_rng(seed)
    n = 64                                     # ৬৪টি জেলা, প্রতিটি এক সারি
    promo = rng.integers(0, 2, n)              # সেই মাসে ছাড়ের ক্যাম্পেইন (1 = হ্যাঁ)
    dist = rng.uniform(0.1, 4.0, n).round(2)   # ঢাকার hub থেকে দূরত্ব, 100 km এককে
    days = rng.integers(5, 31, n)              # মাসে কত দিন ডেলিভারি খোলা ছিল
    rate = np.exp(1.2 + 0.4 * promo - 0.6 * dist)   # প্রতি দিনের আসল গড় অর্ডার
    mu = rate * days
    if nb_alpha:                               # জেলাভেদে অদেখা পার্থক্য: overdispersion
        mu = mu * rng.gamma(1 / nb_alpha, nb_alpha, n)
    return pd.DataFrame({"orders": rng.poisson(mu), "days": days,
                         "promo": promo, "dist": dist})

Poisson regression-এর model:

orders_i ~ Poisson(μ_i),   log μ_i = log(days_i) + β₀ + β₁·promo_i + β₂·dist_i

log(days)-কে বলে offset, আর days হলো exposure: ঘটনা ঘটার “সুযোগ” কতটুকু ছিল। দুই দিকে log(days) বিয়োগ করলে log(μ_i/days_i) = Xβ, অর্থাৎ model আসলে প্রতি দিনের হার নিয়ে, অথচ ডেটা থাকে গোনা সংখ্যায়। তিনটি ভুল বিকল্পের সঙ্গে তুলনা করুন। (ক) days বাদ দিলে ৯ দিন খোলা জেলাকে ৩০ দিনের জেলার সঙ্গে সরাসরি তুলনা করা হয়, আর exposure-এর পার্থক্য promo বা dist-এর ঘাড়ে চাপে। (খ) days-কে সাধারণ feature বানালে model তার coefficient নিজে fit করে; সেটি ১ না হলে “দ্বিগুণ দিন খোলা মানে দ্বিগুণ অর্ডার” কথাটিই ভাঙে। (গ) orders/days-এ সাধারণ regression করলে ৫ দিনের জেলার অনিশ্চিত হার আর ৩০ দিনের জেলার স্থির হার সমান ওজন পায়। Offset-সহ Poisson এই তিন সমস্যাই এড়ায়।

Fit করা হয় maximum likelihood দিয়ে। log y! β-মুক্ত বলে বাদ দিলে log-likelihood:

ℓ(β) = Σ_i [ y_iη_i − e^η_i ],   ∇ℓ = Xᵀ(y − μ),   ∇²ℓ = −XᵀWX, W = diag(μ_i)

Gradient-এর রূপটি logistic regression-এর Xᵀ(y − p)-এর হুবহু ভাই, কারণ দুটিই canonical link. Hessian −XᵀWX সবসময় negative semidefinite, তাই ℓ concave আর পর্ব ১৫-এর Newton-এর পদ্ধতি এখানে খুব ভালো চলে। পরের অংশে সেই Newton-ই দেখবেন, পরিসংখ্যানবিদদের প্রিয় পোশাকে।

৫ · IRLS — Newton-এর ধাপ weighted least squares-এর রূপে

Newton-এর ধাপ β_নতুন = β + (XᵀWX)⁻¹Xᵀ(y − μ). এখানে Xᵀ(y − μ)-কে XᵀW·W⁻¹(y − μ) লিখে একটু সাজালে:

β_নতুন = (XᵀWX)⁻¹ XᵀW z,   z = η + (y − μ)/μ

এটি ঠিক একটি weighted least squares-এর সমাধান: “লক্ষ্য” z-কে X দিয়ে regression, প্রতিটি সারির ওজন W_ii = μ_i. z-কে বলে working response: বর্তমান η থেকে শুরু করে link-এর ঢাল ধরে y-এর দিকে একটি সরলরৈখিক সংশোধন। প্রতি ধাপে μ বদলায়, তাই ওজনও বদলায়: iteratively reweighted least squares, সংক্ষেপে IRLS. Offset থাকলে z থেকে offset বাদ দিয়ে regression করা হয়, কারণ offset-এর coefficient fit হয় না। যেকোনো link আর বণ্টনে সাধারণ সূত্র W_ii = ১/(V(μ_i)·g′(μ_i)²) আর z = η + (y − μ)·g′(μ). Canonical link-এ এটি হুবহু Newton; অন্য link-এ এটি Hessian-এর জায়গায় তার প্রত্যাশিত মান (Fisher information) ব্যবহার করে, যাকে বলে Fisher scoring. নিচের stepper চারটি জেলার ছোট ডেটায় একটি ধাপ হাতে চালায়।

এবার পুরো ৬৪ জেলায় তিনটি পথ পাশাপাশি: হাতে লেখা IRLS (β = ০ থেকে, ইচ্ছা করে দূর থেকে), statsmodels-এর GLM (default পদ্ধতি IRLS), আর scikit-learn-এর PoissonRegressor. শেষেরটিতে offset দেওয়ার সরাসরি parameter নেই; তার বদলে হার y/days-কে লক্ষ্য আর days-কে sample_weight করলে log-likelihood-এর β-নির্ভর অংশ হুবহু একই হয়, তাই উত্তরও একই। আর alpha=0 দিতেই হবে: scikit-learn-এ default alpha=1.0, অর্থাৎ L2 শাস্তি (পর্ব ১৭-এর ridge), যা coefficient-কে শূন্যের দিকে টানে।

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.linear_model import PoissonRegressor
from bazar_district_orders import make_orders

df = make_orders()                                   # ৬৪ জেলা, এক মাস
X = np.column_stack([np.ones(64), df.promo, df.dist])
y, off = df.orders.to_numpy(float), np.log(df.days.to_numpy(float))

b = np.zeros(3)                                      # IRLS, শুরু β = 0
for it in range(1, 30):
    eta = X @ b + off                                # linear predictor + offset
    mu = np.exp(eta)                                 # inverse link
    z = (eta - off) + (y - mu) / mu                  # working response
    b_new = np.linalg.solve(X.T @ (mu[:, None] * X), X.T @ (mu * z))
    print(f"iter {it}: beta = {np.round(b_new, 6)}")
    if np.max(np.abs(b_new - b)) < 1e-10: break
    b = b_new

res = sm.GLM(y, X, family=sm.families.Poisson(), offset=off).fit()
skl = PoissonRegressor(alpha=0, solver="newton-cholesky", tol=1e-12)
skl.fit(df[["promo", "dist"]], y / df.days, sample_weight=df.days)

names = ["const", "promo", "dist"]
print(pd.DataFrame({"IRLS": b, "statsmodels": res.params,
      "sklearn": np.r_[skl.intercept_, skl.coef_]}, index=names).round(4))
ci = np.exp(res.conf_int())                          # 95% CI, তারপর exp
print(pd.DataFrame({"SE": res.bse, "rate_ratio": np.exp(res.params),
      "low": ci[:, 0], "high": ci[:, 1]}, index=names).round(4))
print("statsmodels IRLS iterations:", res.fit_history["iteration"])

iter 1: beta = [ 1.985814  0.600538 -0.859639]
iter 2: beta = [ 1.486857  0.433004 -0.701539]
iter 3: beta = [ 1.314721  0.331727 -0.622724]
iter 4: beta = [ 1.297698  0.317158 -0.614058]
iter 5: beta = [ 1.297526  0.316965 -0.613968]
iter 6: beta = [ 1.297526  0.316965 -0.613968]
iter 7: beta = [ 1.297526  0.316965 -0.613968]
         IRLS  statsmodels  sklearn
const  1.2975       1.2975   1.2975
promo  0.3170       0.3170   0.3170
dist  -0.6140      -0.6140  -0.6140
           SE  rate_ratio     low    high
const  0.0545      3.6602  3.2896  4.0727
promo  0.0483      1.3730  1.2489  1.5093
dist   0.0273      0.5412  0.5130  0.5709
statsmodels IRLS iterations: 5

৬ · Coefficient পড়া — rate ratio আর তার অনিশ্চয়তা

Log link-এ η-তে β যোগ হওয়া মানে μ-তে e^β গুণ হওয়া। তাই Poisson regression-এর coefficient সরাসরি পড়া হয় না; পড়া হয় e^β, যাকে বলে rate ratio (চিকিৎসা-গবেষণায় incidence rate ratio, IRR): অন্য সব feature স্থির রেখে এই feature এক একক বাড়লে প্রতি একক exposure-এ ঘটনার হার কত গুণ হয়। Logistic regression-এর odds ratio-র সঙ্গে এর মিল স্পষ্ট, তবে odds নয়, সরাসরি হার।

দুটি বাড়তি সতর্কতা। প্রথমত, এই CI-গুলো Wald interval, অর্থাৎ “coefficient ± ১.৯৬ × SE”; নমুনা ছোট হলে বা অনুমানগুলো সন্দেহজনক হলে bootstrap দিয়ে যাচাই করা ভালো, পর্ব ২১ (Bootstrap আর resampling — model-এর অনিশ্চয়তা মাপা)-এ যা আসবে। দ্বিতীয়ত, intercept-এর e^β0 = ৩.৬৬ মানে “dist = ০ আর ক্যাম্পেইন নেই এমন জেলার দৈনিক হার”; ডেটায় dist-এর সর্বনিম্ন ০.১৭, তাই এটি সামান্য বাইরে টেনে আনা অনুমান। Feature-কে কেন্দ্রে আনলে (গড় বাদ দিলে) intercept পড়া সহজ হয়।

৭ · গোনা সংখ্যায় OLS কেন ভুল পথে যায়

অনেকে প্রথমে সবচেয়ে পরিচিত পথ নেন: orders-কে লক্ষ্য করে সাধারণ least squares (OLS), আর days-কে একটি feature. নিচের ব্লক একই ডেটায় সেটি আর Poisson GLM পাশাপাশি চালায়, তারপর residual-এর ছবি আঁকে। ছবির লেখা ইংরেজিতে, কারণ পাতার matplotlib বাংলা অক্ষর ঠিকমতো আঁকে না।

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from bazar_district_orders import make_orders

df = make_orders()
X = sm.add_constant(df[["promo", "dist", "days"]].astype(float))
ols = sm.OLS(df.orders, X).fit()                  # সাধারণ least squares, গোনা সংখ্যায়
Xp = sm.add_constant(df[["promo", "dist"]].astype(float))
poi = sm.GLM(df.orders, Xp, family=sm.families.Poisson(),
             offset=np.log(df.days)).fit()

f_ols, f_poi = ols.fittedvalues, poi.fittedvalues
print("OLS: negative predictions =", int((f_ols < 0).sum()), "| min =", round(f_ols.min(), 2))
# নতুন জেলা: hub থেকে 380 km, ক্যাম্পেইন নেই, 8 দিন খোলা
new = np.array([1.0, 0, 3.8, 8])
print("new district, OLS forecast =", round(float(new @ ols.params), 2))
print("new district, Poisson forecast =", round(float(np.exp(new[:3] @ poi.params) * 8), 2))
lo, hi = f_ols < f_ols.median(), f_ols >= f_ols.median()
print("OLS residual SD, low half =", round(ols.resid[lo].std(), 1),
      "| high half =", round(ols.resid[hi].std(), 1))
r_p = poi.resid_pearson
print("Pearson residual SD, low half =", round(r_p[f_poi < f_poi.median()].std(), 2),
      "| high half =", round(r_p[f_poi >= f_poi.median()].std(), 2))

fig, ax = plt.subplots(1, 2, figsize=(8, 3.4))
ax[0].scatter(f_ols, ols.resid, s=14); ax[0].axhline(0, color="k", lw=0.8)
ax[0].axvline(0, color="r", ls="--", lw=0.8)
ax[0].set(title="OLS on counts", xlabel="fitted orders", ylabel="raw residual")
ax[1].scatter(f_poi, r_p, s=14, color="tab:green"); ax[1].axhline(0, color="k", lw=0.8)
ax[1].set(title="Poisson GLM", xlabel="fitted orders", ylabel="Pearson residual")
fig.subplots_adjust(wspace=0.35, bottom=0.18); plt.show()

OLS: negative predictions = 7 | min = -12.85
new district, OLS forecast = -18.53
new district, Poisson forecast = 2.84
OLS residual SD, low half = 8.1 | high half = 11.4
Pearson residual SD, low half = 1.14 | high half = 0.95

আরেকটি জনপ্রিয় বিকল্প হলো log(y + ১)-এ OLS. এতে ঋণাত্মক অনুমান যায়, কিন্তু তিনটি সমস্যা থাকে: “+১” একটি ইচ্ছামতো সংখ্যা আর ছোট count-এ ফল তার ওপর অনেক নির্ভর করে; E[log y] থেকে ফিরে exp করলে গড় নয়, তার চেয়ে কম কিছু পাওয়া যায় (Jensen-এর অসমতা), তাই মোট অর্ডারের অনুমান কম আসে; আর heteroscedasticity থাকলে coefficient নিজেও বাঁকা হতে পারে। অর্থনীতির trade data নিয়ে Santos Silva আর Tenreyro-র ২০০৬ সালের পরিচিত প্রবন্ধ “The Log of Gravity” ঠিক এই যুক্তিতে log-linear OLS-এর বদলে Poisson-ভিত্তিক অনুমানের পক্ষে মত দেয়, এমনকি লক্ষ্য count না হলেও।

৮ · Deviance, residual আর overdispersion — negative binomial পর্যন্ত

Linear regression-এ fit কতটা ভালো, তা মাপে residual-এর বর্গের যোগফল। GLM-এ তার সাধারণ রূপ deviance:

D = 2 [ ℓ(saturated) − ℓ(model) ],   Poisson-এ D = 2 Σ [ y log(y/μ) − (y − μ) ]

Saturated model মানে প্রতিটি সারির জন্য আলাদা parameter, যা প্রতিটি y হুবহু মেলায়; অর্থাৎ সবচেয়ে ভালো সম্ভব likelihood. Deviance তাই “নিখুঁত মেলানো থেকে কতটা দূরে”। Gaussian-এ D হুবহু residual-এর বর্গের যোগফল (σ² = ১ ধরে), আর Bernoulli-তে D = ২ × মোট log loss; তাই AI-এর ভিত্তি কোর্সের পর্ব ৩৮-এর cross-entropy আসলে Bernoulli deviance-এর অর্ধেক। Nested দুটি model-এর deviance-এর পার্থক্য likelihood ratio test-এর পরিসংখ্যান; model তুলনার পুরো আলোচনা পর্ব ২২ (Model তুলনার পরিসংখ্যান — paired test, McNemar আর একাধিক তুলনা)-এ। ভালো fit-এ D মোটামুটি residual-এর degrees of freedom (n − p)-এর কাছাকাছি থাকার কথা, যদিও ছোট count-এ এই নিয়ম দুর্বল।

দ্বিতীয় যাচাই Pearson χ² = Σ (y − μ)²/V(μ), অর্থাৎ Pearson residual-এর বর্গের যোগফল। Poisson ঠিক হলে প্রতিটি পদের গড় প্রায় ১, তাই Pearson χ²/df ≈ ১. এর চেয়ে অনেক বেশি হলে variance গড়ের চেয়ে বড়: overdispersion. বাস্তবে এটিই নিয়ম, ব্যতিক্রম নয়। দুটি জেলার দূরত্ব আর ক্যাম্পেইন এক হলেও একটিতে বড় একটি অফিসপাড়া থাকতে পারে, যা feature-এ নেই; ঈদের আগে অর্ডার দলবেঁধে আসে; একজন পাইকার একাই ২০টি অর্ডার দেন। এমন অদেখা পার্থক্য μ-কে নিজেই এলোমেলো করে, আর তাতে variance বাড়ে। নিচের ব্লক দুটি ডেটা নেয়: আগের বিশুদ্ধ Poisson ডেটা, আর একই module-এ প্রতি জেলার μ-কে একটি এলোমেলো gamma গুণক দিয়ে নাড়ানো ডেটা, যা nb_alpha=0.5 দিয়ে বানানো হয়। তারপর দ্বিতীয়টিতে negative binomial model, যার variance μ + αμ².

import numpy as np
import statsmodels.api as sm
from scipy.special import xlogy
from bazar_district_orders import make_orders

def poisson_fit(df):
    X = sm.add_constant(df[["promo", "dist"]].astype(float))
    return X, sm.GLM(df.orders, X, family=sm.families.Poisson(),
                     offset=np.log(df.days)).fit()

for name, df in [("plain Poisson data", make_orders()),
                 ("overdispersed data", make_orders(nb_alpha=0.5))]:
    X, res = poisson_fit(df)
    y, mu = df.orders, res.fittedvalues
    dev = 2 * np.sum(xlogy(y, y / mu) - (y - mu))       # Poisson deviance
    pear = np.sum((y - mu) ** 2 / mu)                  # Pearson χ²
    print(f"{name}: deviance {dev:.1f} (statsmodels {res.deviance:.1f}),"
          f" df {res.df_resid:.0f}")
    print(f"   Pearson chi2/df = {pear / res.df_resid:.2f}"
          f" (statsmodels {res.pearson_chi2 / res.df_resid:.2f})")

nb = sm.NegativeBinomial(df.orders, X, offset=np.log(df.days)).fit(disp=0)
print("NB converged:", nb.mle_retvals["converged"], "| alpha =", round(nb.params["alpha"], 3))
print("promo coef  Poisson", round(res.params["promo"], 3), "SE", round(res.bse["promo"], 3))
print("promo coef  NegBin ", round(nb.params["promo"], 3), "SE", round(nb.bse["promo"], 3))
print("Poisson CI for rate ratio:", np.exp(res.conf_int().loc["promo"]).round(2).tolist())
print("NegBin  CI for rate ratio:", np.exp(nb.conf_int().loc["promo"]).round(2).tolist())

plain Poisson data: deviance 73.7 (statsmodels 73.7), df 61
   Pearson chi2/df = 1.14 (statsmodels 1.14)
overdispersed data: deviance 727.7 (statsmodels 727.7), df 61
   Pearson chi2/df = 12.19 (statsmodels 12.19)
NB converged: True | alpha = 0.427
promo coef  Poisson 0.187 SE 0.05
promo coef  NegBin  0.415 SE 0.178
Poisson CI for rate ratio: [1.09, 1.33]
NegBin  CI for rate ratio: [1.07, 2.15]

Overdispersion-এর প্রতিকার কয়েকটি, প্রতিটির নিজের অনুমান আছে:

একটি ফাঁদ: sm.families.NegativeBinomial()-এ α না দিলে statsmodels নিজে α = ১ ধরে নেয় আর একটি ValueWarning দেখায়। GLM-এর NB family α-কে অনুমান করে না; α হয় নিজে দিতে হয়, নয়তো উপরের মতো sm.NegativeBinomial (discrete model) দিয়ে α-ও MLE করতে হয়।

৯ · AI-এর যুগে — AI-এর লেখা GLM কোড যাচাই

“অর্ডারের সংখ্যা predict করার model বানাও” বললে AI assistant প্রায়ই এমন কোড দেয় যা চলে, সংখ্যা ছাপে, আর দেখতে বিশ্বাসযোগ্য; কিন্তু GLM-এর ভুলগুলো চুপচাপ থাকে। নিচের ভুলগুলো বারবার দেখা যায়।

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

• লক্ষ্যের ধরন দেখে বণ্টন আর link বেছেছি, আর fitted মান ঋণাত্মক হয় না তা দেখেছি
• Exposure থাকলে offset দিয়েছি, আর ছোট একটি উদাহরণে হাতে η আর μ মিলিয়েছি
• দুটি library (বা হাতে লেখা IRLS) একই coefficient দেয় কি না মিলিয়েছি, আর scikit-learn-এ alpha কত তা দেখেছি
• Coefficient-কে rate ratio আর CI হিসেবে লিখেছি, আর পর্যবেক্ষণের ডেটা থেকে কারণ-ফল দাবি করিনি
• Pearson χ²/df আর residual-এর ছবি দেখেছি; overdispersion থাকলে SE ঠিক করেছি
• ডেটা আগে ভাগ করেছি, test set তালা দেওয়া; held-out Poisson deviance একটি baseline (যেমন শুধু গড় হার)-এর সঙ্গে কয়েকটি fold-এ তুলনা করেছি
• নিজের version-এর docs পড়েছি, warning চাপা দিইনি, আর কোনো আসল গ্রাহকের ডেটা AI-তে পেস্ট করিনি