Bayesian linear regression — posterior আর predictive অনিশ্চয়তা

AI-এর ভিত্তি কোর্সের পর্ব ৫০-এ prior, likelihood, evidence আর posterior দেখেছেন, Beta–Binomial দিয়ে একটি সংখ্যার posterior আর credible interval-ও। এই কোর্সের পর্ব ১৮-এ MAP মানে prior যোগ করা MLE, পর্ব ১৭-এ ridge মানে Gaussian prior-এর MAP, আর পর্ব ১৯-এ বহুমাত্রিক Gaussian. কিন্তু MAP শেষ পর্যন্ত একটি বিন্দু। এই পর্বের প্রশ্ন: বাজারঘর ক্রেতাকে বলবে “আপনার পার্সেল ৪৫ থেকে ৬০ মিনিটে আসবে”; সেই পরিসর কত চওড়া হওয়া উচিত, আর ডেটা যেখানে কম সেখানে model কি নিজের অজ্ঞতা স্বীকার করবে? চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, SciPy, scikit-learn 1.8.0, statsmodels); প্রথমবার scikit-learn নামতে প্রায় ১০ সেকেন্ড লাগে।

১ · একটি রেখা বনাম রেখার বণ্টন — model আর prior

Model. প্রতিটি অর্ডারের input x (যেমন দূরত্ব) থেকে কয়েকটি feature বানাই, φ(x) = (1, x) বা (1, x, x²)। এগুলোকে basis function বলে (পর্ব ২৩)। তারপর ধরে নিই

y = wᵀφ(x) + ε,   ε ~ N(0, β⁻¹)

এখানে w হলো M-টি weight-এর vector, ε হলো noise. β হলো noise-এর precision, অর্থাৎ variance-এর উল্টো: β = 1/σ². σ = ৬ মিনিট হলে β = 1/36. N-টি অর্ডারের feature-গুলো সারি হিসেবে সাজালে পাই N × M design matrix Φ. Noise স্বাধীন ধরলে likelihood হলো p(y | w) = N(y | Φw, β⁻¹I), যার log হলো −(β/2)‖y − Φw‖² + ধ্রুবক। এটুকু সাধারণ least squares-এর গল্প।

Prior. নতুন অংশ: w-কেও একটি random variable ধরি, ডেটা দেখার আগে

w ~ N(0, α⁻¹I)

α হলো weight-এর prior precision. α = 1/400 মানে প্রতিটি weight-এর prior standard deviation ২০: “ইন্টারসেপ্ট বা প্রতি কিমির মিনিট সম্ভবত −৪০ থেকে +৪০-এর মধ্যে”, খুব ঢিলেঢালা একটি বিশ্বাস। α বড় হলে prior weight-কে শূন্যের দিকে কড়া টানে। α আর β, দুটোই weight নয়, weight-এর বণ্টন নিয়ন্ত্রণ করে; তাই এদের hyperparameter বলে।

লক্ষ্য. Bayes-এর নিয়মে posterior p(w | y) ∝ p(y | w) p(w). MLE বা MAP এই বণ্টনের একটি বিন্দু (চূড়া) নেয়; Bayesian পদ্ধতি পুরো বণ্টন রাখে, আর নতুন অর্ডারের অনুমান করার সময় সব সম্ভাব্য w-এর ওপর গড় করে। সেই গড় থেকেই আসে “৪৫ থেকে ৬০ মিনিট” ধরনের পরিসর।

২ · Completing the square — posterior-এর গড় আর covariance

Prior আর likelihood দুটোই w-এর দ্বিঘাত রাশির exponential. গুণ করলে log-এ যোগ হয়:

log p(w | y) = −(β/2)‖y − Φw‖² − (α/2)wᵀw + ধ্রুবক

প্রথম পদ খুলে লিখি: ‖y − Φw‖² = yᵀy − 2wᵀΦᵀy + wᵀΦᵀΦw. w-এর ওপর নির্ভর করে এমন পদগুলো সাজালে

log p(w | y) = −½ wᵀ(αI + βΦᵀΦ)w + wᵀ(βΦᵀy) + ধ্রুবক

এবার একটি Gaussian N(m, S)-এর log density খুলে দেখুন: −½(w − m)ᵀS⁻¹(w − m) = −½wᵀS⁻¹w + wᵀS⁻¹m + ধ্রুবক। দুটি রাশি পাশাপাশি রেখে দ্বিঘাত পদ আর একঘাত পদ মেলালেই উত্তর। এই কৌশলের নাম completing the square, স্কুলের ax² + bx-কে a(x + b/2a)²-এ আনার বহুমাত্রিক রূপ।

Prior আর posterior একই পরিবারের (দুটিই Gaussian) বলে Gaussian prior এখানে conjugate prior. ফলে নতুন ডেটা এলে আজকের posterior-কে কালকের prior ধরে একই সূত্র আবার খাটানো যায়: S_new⁻¹ = S_old⁻¹ + βφφᵀ. একবারে সব ডেটা দিলে আর একটি একটি করে দিলে শেষ ফল হুবহু এক। এটিই Bayesian updating. সংখ্যার দিক থেকে S_N⁻¹ সবসময় symmetric positive definite (α > 0 থাকায়), তাই inverse না নিয়ে Cholesky দিয়ে সমাধান করাই নিরাপদ (পর্ব ১০)।

এক মাত্রায় (একটিই weight, ইন্টারসেপ্ট নেই, y = wx + ε) সূত্র দুটি এত ছোট যে খাতায় করা যায়: precision = α + βΣx², আর গড় = βΣxy ÷ precision. নিচের stepper তিনটি বিন্দু একে একে যোগ করে।

৩ · একটি একটি করে অর্ডার — posterior কীভাবে সরু হয়

এবার বাজারঘরের ২০টি কাল্পনিক অর্ডার (seed দিয়ে বানানো): দূরত্ব ১ থেকে ১০ কিমি, আসল সম্পর্ক সময় = ২০ + ৫ × দূরত্ব, সঙ্গে σ = ৬ মিনিটের noise. আপাতত ধরে নিচ্ছি β = 1/36 জানা (অভিজ্ঞতা থেকে), আর α = 1/400. নিচের ব্লক ১, ২, ৫ আর ২০টি অর্ডার দেখার পর posterior হিসাব করে। ওপরের সারিতে weight-এর জগতে (w0 = ইন্টারসেপ্ট, w1 = প্রতি কিমিতে মিনিট) posterior-এর contour, লাল + আসল (২০, ৫)। নিচের সারিতে posterior থেকে আঁকা ছয়টি রেখা আর এ পর্যন্ত দেখা অর্ডার। Inverse-এর জন্য scipy.linalg, কারণ ব্রাউজারের NumPy অনেক linear algebra ব্যর্থতা নীরবে NaN দিয়ে ঢেকে দেয়।

import numpy as np
import scipy.linalg as sla
import matplotlib.pyplot as plt

rng = np.random.default_rng(45)
x = rng.uniform(1, 10, 20)                      # দূরত্ব (km)
y = 20 + 5 * x + rng.normal(0, 6, 20)          # ডেলিভারি সময় (মিনিট)
alpha, beta = 1 / 20**2, 1 / 6**2              # prior precision, noise precision

def posterior(xs, ys):
    Phi = np.column_stack([np.ones_like(xs), xs])
    A = alpha * np.eye(2) + beta * Phi.T @ Phi  # posterior precision
    S = sla.inv(A)                              # posterior covariance
    m = beta * S @ Phi.T @ ys                   # posterior mean
    return m, S

W0, W1 = np.meshgrid(np.linspace(-20, 60, 160), np.linspace(-4, 14, 160))
grid = np.dstack([W0, W1])
fig, ax = plt.subplots(2, 4, figsize=(11, 5.2))
for j, n in enumerate([1, 2, 5, 20]):
    m, S = posterior(x[:n], y[:n])
    sd = np.sqrt(np.diag(S)); r = S[0, 1] / (sd[0] * sd[1])
    print(f"n={n:2d}  mean=({m[0]:6.2f}, {m[1]:5.2f})  sd=({sd[0]:5.2f}, {sd[1]:4.2f})  corr={r:+.2f}")
    d = grid - m
    q = np.einsum("...i,ij,...j->...", d, sla.inv(S), d)
    ax[0, j].contour(W0, W1, np.exp(-q / 2), levels=[0.05, 0.3, 0.7])
    ax[0, j].plot(20, 5, "r+", ms=10)
    ax[0, j].set_title(f"posterior, n = {n}"); ax[0, j].set_xlabel("w0 (intercept)")
    for w in rng.multivariate_normal(m, S, 6):
        ax[1, j].plot([0, 10], [w[0], w[0] + 10 * w[1]], color="tab:blue", alpha=0.5)
    ax[1, j].plot(x[:n], y[:n], "ko", ms=4); ax[1, j].set_ylim(0, 90)
    ax[1, j].set_xlabel("distance (km)")
ax[0, 0].set_ylabel("w1 (min per km)"); ax[1, 0].set_ylabel("delivery time (min)")
fig.subplots_adjust(hspace=0.45, wspace=0.3)
plt.show()

n= 1  mean=(  1.26,  7.79)  sd=(19.74, 3.34)  corr=-0.96
n= 2  mean=(  3.97,  7.91)  sd=(19.49, 3.34)  corr=-0.98
n= 5  mean=( 13.90,  5.44)  sd=(12.56, 1.86)  corr=-0.98
n=20  mean=( 15.77,  5.23)  sd=( 5.07, 0.72)  corr=-0.96

৪ · Posterior mean মানেই ridge — আর predictive বণ্টন

m_N = β(αI + βΦᵀΦ)⁻¹Φᵀy. বন্ধনীর ভেতর থেকে β বের করে আনলে β কেটে যায়:

m_N = (ΦᵀΦ + λI)⁻¹Φᵀy,   λ = α/β

এটি হুবহু ridge regression-এর সমাধান (পর্ব ১৭-এ λ = σ²/τ² হিসেবে দেখেছেন; τ² = 1/α prior-এর variance, σ² = 1/β)। অর্থাৎ ridge-এর penalty ইচ্ছেমতো একটি সংখ্যা নয়: “noise-এর variance ভাগ prior-এর variance”. Noise বেশি হলে ডেটা কম বিশ্বাসযোগ্য, তাই prior-এর টান বেশি। একটি সূক্ষ্ম পার্থক্য: এখানে prior ইন্টারসেপ্টেও বসেছে, তাই তুলনা করতে হবে fit_intercept=False দিয়ে, ১-এর কলামটি Φ-এর ভেতরে রেখে। সাধারণ ridge ইন্টারসেপ্টকে শাস্তি দেয় না।

নতুন অর্ডারের অনুমান. নতুন দূরত্ব x*-এর feature φ* = φ(x*)। সময় y* = wᵀφ* + ε, যেখানে w আসে posterior থেকে আর ε নতুন noise. Gaussian-এর রৈখিক রূপান্তর আবার Gaussian, তাই সব w-এর ওপর integral করে পাই posterior predictive distribution:

p(y* | x*, ডেটা) = N( m_Nᵀφ*,  1/β + φ*ᵀS_Nφ* )

Variance-এর প্রথম অংশ 1/β হলো aleatoric uncertainty: ডেটার নিজের এলোমেলোভাব, w পুরো জানা থাকলেও থাকত। দ্বিতীয় অংশ φ*ᵀS_Nφ* হলো epistemic uncertainty: আমরা w ঠিক জানি না বলে যে অনিশ্চয়তা, বেশি ডেটায় কমে। “Plug-in” অনুমান শুধু w = m_N বসিয়ে দেয় আর variance ধরে 1/β, দ্বিতীয় অংশ বাদ। নিচের ব্লক আগে ridge-এর সঙ্গে মিল যাচাই করে, তারপর ০ থেকে ২৫ কিমি পর্যন্ত predictive পরিসর হিসাব করে আর statsmodels-এর OLS prediction interval-এর পাশে রাখে (পর্ব ২৩)।

import numpy as np
import scipy.linalg as sla
from sklearn.linear_model import Ridge
import statsmodels.api as sm

rng = np.random.default_rng(45)
x = rng.uniform(1, 10, 20)
y = 20 + 5 * x + rng.normal(0, 6, 20)
alpha, beta = 1 / 20**2, 1 / 6**2
Phi = np.column_stack([np.ones_like(x), x])
S = sla.inv(alpha * np.eye(2) + beta * Phi.T @ Phi)
m = beta * S @ Phi.T @ y

lam = alpha / beta                              # 36/400 = 0.09
ridge = Ridge(alpha=lam, fit_intercept=False).fit(Phi, y)
print("posterior mean :", m.round(4))
print("ridge (λ=α/β)  :", ridge.coef_.round(4), " λ =", round(lam, 4))
print("সর্বোচ্চ পার্থক্য:", f"{np.abs(m - ridge.coef_).max():.1e}")

ols = sm.OLS(y, Phi).fit()
print("OLS-এর অনুমিত σ:", round(np.sqrt(ols.scale), 2), "(আসল 6)")
xt = np.array([0.0, 2.0, 6.0, 10.0, 15.0, 25.0])
Pt = np.column_stack([np.ones_like(xt), xt])
mean = Pt @ m
epi = np.einsum("ij,jk,ik->i", Pt, S, Pt)      # φᵀ S φ : weight-এর অনিশ্চয়তা
tot = 1 / beta + epi                            # + noise variance
pi = ols.get_prediction(Pt).summary_frame(alpha=0.05)
print("\n  km   mean  sd_epi  sd_tot  Bayes 95%        OLS 95% PI")
for i, v in enumerate(xt):
    lo, hi = mean[i] - 1.96 * np.sqrt(tot[i]), mean[i] + 1.96 * np.sqrt(tot[i])
    print(f"{v:4.0f} {mean[i]:6.1f} {np.sqrt(epi[i]):6.2f} {np.sqrt(tot[i]):6.2f}"
          f"  [{lo:5.1f}, {hi:5.1f}]  [{pi.obs_ci_lower[i]:5.1f}, {pi.obs_ci_upper[i]:5.1f}]")
print("plug-in (শুধু noise): ±", round(1.96 * 6, 1), "সব দূরত্বে একই")

posterior mean : [15.7713  5.2304]
ridge (λ=α/β)  : [15.7713  5.2304]  λ = 0.09
সর্বোচ্চ পার্থক্য: 6.9e-14
OLS-এর অনুমিত σ: 5.22 (আসল 6)

  km   mean  sd_epi  sd_tot  Bayes 95%        OLS 95% PI
   0   15.8   5.07   7.86  [  0.4,  31.2]  [  2.2,  31.4]
   2   26.2   3.70   7.05  [ 12.4,  40.1]  [ 14.0,  40.0]
   6   47.2   1.46   6.17  [ 35.1,  59.3]  [ 36.0,  58.6]
  10   68.1   2.67   6.57  [ 55.2,  80.9]  [ 55.7,  79.7]
  15   94.2   6.06   8.53  [ 77.5, 110.9]  [ 77.3, 108.9]
  25  146.5  13.18  14.49  [118.1, 174.9]  [116.9, 171.2]
plug-in (শুধু noise): ± 11.8 সব দূরত্বে একই

৫ · Predictive variance ভেঙে দেখা — কোন অংশ কখন বড়

φ = (1, x) হলে φᵀSφ = S₀₀ + 2x S₀₁ + x² S₁₁, x-এর একটি দ্বিঘাত রাশি, যার সর্বনিম্ন বিন্দু ডেটার গড়ের কাছে। তাই পট্টির আকার ঘণ্টাঘড়ির মতো। নিচের stepper ব্লকের S_N দিয়ে ৬ কিমি আর ২৫ কিমিতে হিসাবটি খুলে দেখায়।

Frequentist interval-এর সঙ্গে সম্পর্ক. OLS-এর prediction interval-এর সূত্রও s²(1 + φ*ᵀ(ΦᵀΦ)⁻¹φ*), একই দুই অংশ। Prior খুব ঢিলে (α → ০) আর σ অজানা ধরে তার ওপরও একটি “তথ্যহীন” prior দিলে Bayesian predictive হুবহু Student-t prediction interval-এ নেমে আসে। তাই পার্থক্য সংখ্যায় নয়, বাক্যে: frequentist বলে “এই পদ্ধতিতে বানানো পরিসর বারবার ব্যবহার করলে ৯৫% বার নতুন অর্ডার ধরবে”; Bayesian বলে “এই model আর prior মেনে নিলে, এই অর্ডারটি এই পরিসরে থাকার সম্ভাবনা ৯৫%”. Prior কড়া হলে বা ডেটা কম হলে সংখ্যাও আলাদা হয়, আর তখনই prior-এর পছন্দ যাচাই করা জরুরি। দুই ক্ষেত্রেই আসল পরীক্ষা একই: আলাদা রাখা অর্ডারে পরিসরটি সত্যিই প্রায় ৯৫% বার ধরে কি না (coverage), যেমন পর্ব ২৬-এ সম্ভাবনার calibration দেখেছেন।

৬ · Basis function আর evidence — α, β ডেটা থেকে, ডিগ্রিও

এতক্ষণ α আর β হাতে দিয়েছি। বাস্তবে এগুলো জানা থাকে না। একটি পথ cross-validation (AI-এর ভিত্তি কোর্সের পর্ব ৬৬)। Bayesian পথ আরেকটি প্রশ্ন করে: এই α, β আর এই basis নিয়ে, দেখা ডেটা আসার সম্ভাবনা কত? উত্তর w-এর ওপর integral করে পাওয়া evidence বা marginal likelihood, p(y | α, β) = ∫ p(y | w, β) p(w | α) dw. Gaussian বলে integral বন্ধ রূপে আসে:

ln p(y | α, β) = (M/2)ln α + (N/2)ln β − E(m_N) − ½ ln|S_N⁻¹| − (N/2)ln 2π

যেখানে E(m) = (β/2)‖y − Φm‖² + (α/2)mᵀm, ঠিক ridge-এর লক্ষ্য-ফাংশন, আর |S_N⁻¹| হলো posterior precision matrix-এর determinant. Evidence সর্বোচ্চ করে α, β বাছাকে বলে empirical Bayes বা type-II maximum likelihood: weight-এর ওপর integral, hyperparameter-এর ওপর optimization. ডেরিভেটিভ শূন্য ধরলে MacKay-র (১৯৯২) fixed-point নিয়ম আসে। βΦᵀΦ-এর eigenvalue λ_i হলে

γ = Σ_i λ_i/(α + λ_i),   α ← γ / m_Nᵀm_N,   1/β ← ‖y − Φm_N‖² / (N − γ)

γ হলো effective number of parameters: যে weight-এর দিকে ডেটা prior-এর চেয়ে বেশি কথা বলে (λ_i ≫ α), সেটি প্রায় ১ গোনে, বাকিগুলো প্রায় ০ (পর্ব ১৭-এর effective degrees of freedom-এর সঙ্গে মিলিয়ে দেখুন)। দ্বিতীয় সূত্রটি noise variance-এর অনুমান, শুধু N-এর বদলে N − γ দিয়ে ভাগ, যেমন OLS-এ N − M দিয়ে।

নিচের ব্লকে নতুন ৩০টি অর্ডার, যেখানে আসল সম্পর্ক বাঁকা: সময় = ১২ + ৭x − ০.৩x² (দূরে গেলে রাইডার হাইওয়ে পান, প্রতি কিমির সময় কমে)। ডিগ্রি ০ থেকে ৭ পর্যন্ত polynomial basis (polynomial features), x-কে আগে প্রায় [−১, ১]-এ এনে। প্রতিটি ডিগ্রিতে নিজে লেখা fixed-point দিয়ে α, β আর evidence, পাশে scikit-learn-এর BayesianRidge একই কাজ করে। RBF basis (radial basis function) দিলেও একই কোড চলে, শুধু Φ বদলায়।

import numpy as np
import scipy.linalg as sla
from sklearn.linear_model import BayesianRidge

rng = np.random.default_rng(7)
x = rng.uniform(0.5, 12, 30)                    # দূরত্ব (km)
y = 12 + 7 * x - 0.3 * x**2 + rng.normal(0, 4, 30)
u = (x - 6) / 6                                 # প্রায় [-1, 1]-এ আনা
N = len(y)

def evidence(Phi, a=1.0, b=1.0):
    ev = np.linalg.eigvalsh(Phi.T @ Phi)        # ΦᵀΦ-এর eigenvalue
    for _ in range(500):                        # MacKay-র fixed-point
        A = a * np.eye(Phi.shape[1]) + b * Phi.T @ Phi
        m = b * sla.solve(A, Phi.T @ y, assume_a="pos")
        gamma = np.sum(b * ev / (a + b * ev))   # কার্যকর parameter-সংখ্যা
        sse = np.sum((y - Phi @ m) ** 2)
        a, b = gamma / (m @ m), (N - gamma) / sse
    A = a * np.eye(Phi.shape[1]) + b * Phi.T @ Phi
    m = b * sla.solve(A, Phi.T @ y, assume_a="pos")
    sse = np.sum((y - Phi @ m) ** 2)
    logev = (Phi.shape[1] * np.log(a) + N * np.log(b) - b * sse - a * m @ m
             - np.linalg.slogdet(A)[1] - N * np.log(2 * np.pi)) / 2
    return logev, a, b, np.sqrt(sse / N), gamma

print("deg   logev     α      β     rmse   γ  │ sk_score  sk_λ    sk_α")
for d in range(0, 8):
    Phi = np.vander(u, d + 1, increasing=True)  # 1, u, u², …
    le, a, b, rmse, g = evidence(Phi)
    br = BayesianRidge(fit_intercept=False, compute_score=True, tol=1e-10,
                       max_iter=1000).fit(Phi, y)
    print(f"{d:3d} {le:8.2f} {a:.4f} {b:.4f} {rmse:5.2f} {g:4.1f} │"
          f" {br.scores_[-1]:8.2f} {br.lambda_:.4f} {br.alpha_:.4f}")

deg   logev     α      β     rmse   γ  │ sk_score  sk_λ    sk_α
  0  -120.31 0.0006 0.0068 11.90  1.0 │  -120.31 0.0006 0.0068
  1   -93.53 0.0010 0.0525  4.22  2.0 │   -93.53 0.0010 0.0525
  2   -90.56 0.0013 0.0764  3.43  3.0 │   -90.56 0.0013 0.0764
  3   -91.46 0.0018 0.0783  3.33  3.9 │   -91.46 0.0018 0.0783
  4   -91.73 0.0018 0.0833  3.18  4.8 │   -91.73 0.0018 0.0833
  5   -92.33 0.0020 0.0822  3.18  5.1 │   -92.33 0.0020 0.0822
  6   -92.70 0.0021 0.0821  3.17  5.2 │   -92.70 0.0021 0.0821
  7   -93.10 0.0022 0.0813  3.17  5.4 │   -93.10 0.0022 0.0813

Seed বদলালে কী হয় (সত্যিকারের run). একই generator seed ০ থেকে ৯-এ ডেস্কটপে চালানো হয়েছে (লেখার সময়ে Python 3.14.4, NumPy 2.5.4, scikit-learn 1.8.0; একটি local run, সাধারণ ফল নয়); প্রতিটি সারিতে evidence যে ডিগ্রি বেছেছে আর সেরা দুটি log evidence:

$ python seeds.py
seed 0: best degree 7 | top-2 log evidence -92.33 -92.38
seed 1: best degree 2 | top-2 log evidence -91.21 -92.69
seed 2: best degree 2 | top-2 log evidence -93.04 -93.75
seed 3: best degree 2 | top-2 log evidence -94.05 -95.37
seed 4: best degree 2 | top-2 log evidence -88.55 -89.92
seed 5: best degree 2 | top-2 log evidence -94.89 -96.35
seed 6: best degree 4 | top-2 log evidence -89.20 -89.68
seed 7: best degree 2 | top-2 log evidence -90.56 -91.46
seed 8: best degree 3 | top-2 log evidence -98.81 -99.18
seed 9: best degree 2 | top-2 log evidence -88.47 -89.15

দশটির সাতটিতে আসল ডিগ্রি ২ জিতেছে। বাকি তিনটিতে (seed ০, ৬, ৮) জিতেছে ৭, ৪ আর ৩, কিন্তু প্রথম আর দ্বিতীয়ের ব্যবধান মাত্র ০.০৫, ০.৪৮ আর ০.৩৭ log একক, অর্থাৎ Bayes factor ২-এর কম। Evidence যখন জোর দিয়ে কিছু বলে না, তখন সরল model নেওয়া, বা কয়েকটি ডিগ্রির অনুমান evidence-এর ওজনে মেলানো (Bayesian model averaging), নিরাপদ পথ।

Evidence কেন বেশি নমনীয় model-কে শাস্তি দেয়, যখন training ভুল তার পক্ষে? log evidence-কে দুই ভাগে লেখা যায়: fit পদ (N/2)ln β − (β/2)‖y − Φm_N‖² − (N/2)ln 2π, আর বাকিটা, (M/2)ln α − (α/2)m_Nᵀm_N − ½ln|S_N⁻¹|. দ্বিতীয়টিকে বলে Occam factor-এর log: মোটামুটি “posterior-এর আয়তন ভাগ prior-এর আয়তন”-এর log, প্রতিটি নতুন weight-এর জন্য ঋণাত্মক একটি দাম।

৭ · ARD — প্রতিটি weight-এর নিজস্ব prior

BayesianRidge সব weight-কে একই α দেয়। Automatic relevance determination (ARD) প্রতিটি weight w_j-কে আলাদা precision α_j দেয়: w_j ~ N(0, α_j⁻¹)। তারপর evidence সর্বোচ্চ করে প্রতিটি α_j শেখে, update α_j ← γ_j/m_j², যেখানে γ_j = 1 − α_j(S_N)_jj. কোনো feature ডেটা ব্যাখ্যায় কাজে না লাগলে তার α_j অনেক বড় হয়, prior weight-কে শূন্যে চেপে ধরে, আর Occam দাম বাঁচে। scikit-learn-এর ARDRegression-এ α_j হলো lambda_[j], আর threshold_lambda (default ১০,০০০) ছাড়ালে weight পুরোপুরি বাদ দেয়। নিচের ৬০টি অর্ডারে তিনটি আসল feature (দূরত্ব, ওজন, বৃষ্টি) আর তিনটি অপ্রাসঙ্গিক (অর্ডার নম্বরের শেষ দুই অঙ্ক, দাম, রেটিং)। সব feature আগে standardize করা, কারণ prior-এর precision feature-এর মাপের ওপর নির্ভর করে।

import numpy as np
import pandas as pd
from sklearn.linear_model import ARDRegression, BayesianRidge, LinearRegression
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(11)
n = 60
X = pd.DataFrame({
    "distance_km": rng.uniform(1, 12, n),
    "weight_kg": rng.uniform(0.2, 8, n),
    "rain": rng.integers(0, 2, n).astype(float),
    "order_id_mod": rng.integers(0, 100, n).astype(float),   # অপ্রাসঙ্গিক
    "price_tk": rng.uniform(200, 3000, n),                   # অপ্রাসঙ্গিক
    "rating": rng.integers(1, 6, n).astype(float),           # অপ্রাসঙ্গিক
})
y = (15 + 5 * X.distance_km + 1.5 * X.weight_kg + 8 * X.rain
     + rng.normal(0, 4, n))
Z = StandardScaler().fit_transform(X)            # সব feature এক মাপে

ols = LinearRegression().fit(Z, y)
br = BayesianRidge().fit(Z, y)
ard = ARDRegression().fit(Z, y)
pd.set_option("display.width", 80)
tab = pd.DataFrame({"OLS": ols.coef_, "BayesRidge": br.coef_,
                    "ARD": ard.coef_, "ARD λ_j": ard.lambda_,
                    "ARD sd": np.sqrt(np.diag(ard.sigma_))},
                   index=X.columns)
print(tab.round(3).to_string())
print("ARD noise sd ≈", round(1 / np.sqrt(ard.alpha_), 2),
      "| BayesianRidge-এর একক λ =", round(br.lambda_, 4))

                 OLS  BayesRidge     ARD   ARD λ_j  ARD sd
distance_km   16.556      16.447  16.432     0.004   0.549
weight_kg      3.945       3.912   3.908     0.064   0.548
rain           3.759       3.759   3.750     0.070   0.548
order_id_mod   0.087       0.087   0.000  1270.122   0.028
price_tk      -0.482      -0.457  -0.002   681.149   0.038
rating         0.291       0.290   0.001  1158.056   0.029
ARD noise sd ≈ 4.18 | BayesianRidge-এর একক λ = 0.0199

৮ · Gaussian-এর বাইরে — Laplace approximation, সীমা আর পরের পথ

Bayesian logistic regression. ফেরত (return) হবে কি না, এমন হ্যাঁ-না লক্ষ্যে likelihood Bernoulli, আর Gaussian prior-এর সঙ্গে গুণ করলে posterior আর Gaussian থাকে না; বন্ধ সূত্রও নেই। সবচেয়ে পুরোনো সহজ পথ Laplace approximation: (১) MAP বিন্দু w* খুঁজুন (L2-শাস্তিযুক্ত logistic regression, যেমন LogisticRegression(C=1/α)); (২) সেখানে −log posterior-এর Hessian H = αI + Σ_n σ_n(1 − σ_n)φ_nφ_nᵀ হিসাব করুন, যেখানে σ_n হলো অর্ডার n-এর অনুমিত সম্ভাবনা; (৩) posterior ≈ N(w*, H⁻¹)। অর্থাৎ চূড়ার চারপাশে log posterior-কে একটি দ্বিঘাত রাশি দিয়ে বদলানো, পর্ব ১৮-এর Fisher information-এর ভাবনার মতো। তারপর নতুন অর্ডারের সম্ভাবনায় অনিশ্চয়তা আনতে একটি পরিচিত আনুমানিক সূত্র: p ≈ σ(μ / √(1 + πs²/8)), যেখানে μ = w*ᵀφ আর s² = φᵀH⁻¹φ. s² বড় হলে logit শূন্যের দিকে নামে, অর্থাৎ ডেটা থেকে দূরে model কম আত্মবিশ্বাসী সম্ভাবনা দেয়, যা প্রায়ই calibration-এ সাহায্য করে।

সীমাগুলো সৎভাবে.

• Model ভুল হলে অনিশ্চয়তাও ভুল। Epistemic অংশ “এই basis-এর ভেতরে কোন w” নিয়ে অনিশ্চয়তা; আসল সম্পর্ক বাঁকা হলে সরলরেখার posterior আত্মবিশ্বাসের সঙ্গে ভুল জায়গায় সরু হতে পারে।
• Noise সব জায়গায় সমান ধরা হয়। ঈদের আগে বা বৃষ্টিতে ছড়ানো বেশি হলে (heteroscedastic noise) একটি β যথেষ্ট নয়; পর্ব ২০ বা পর্ব ২৪-এর মতো পথ লাগে।
• Prior-এর মাপ feature-এর মাপের ওপর নির্ভর করে। Standardize না করলে “একই α” আসলে একেক feature-এ একেক শক্তির prior.
• খরচ. S_N একটি M × M matrix; M হাজার ছাড়ালে Cholesky O(M³)। N অনেক বড় হলে সমস্যা নেই (ΦᵀΦ একবারই বানাতে হয়)।

পরের পথ. weight-এর বদলে সরাসরি function-এর ওপর prior দিলে, আর basis-এর সংখ্যা অসীম করলে, আসে পর্ব ৪৬ (Gaussian process — kernel দিয়ে অনিশ্চয়তাসহ অনুমান); পর্ব ৩০-এর kernel ridge-এর গড়ের সঙ্গে তখন variance-ও পাওয়া যায়। Posterior বন্ধ সূত্রে না এলে নমুনা নিয়ে কাজ চালানো হয় পর্ব ৪৭ (MCMC — Metropolis-Hastings আর Gibbs sampling)-এ, আর আনুমানিক বণ্টন দিয়ে পর্ব ৪৮ (Variational inference-এর ধারণা — ELBO আর approximation)-এ। Predictive অনিশ্চয়তা কাজে লাগিয়ে পরের পরীক্ষা কোথায় করবেন, তা পর্ব ৪৯ (Bayesian optimization — অল্প চেষ্টায় hyperparameter খোঁজা)-এর বিষয়।

৯ · AI-এর যুগে — Bayesian regression-এর কোড আর দাবি যাচাই

AI assistant Bayesian regression-এর derivation আর কোড দ্রুত লিখে দেয়, আর সেগুলো দেখতে খুব গাণিতিক বলে ভুল ধরা কঠিন। নিচের ভুলগুলো বারবার দেখা যায়।

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

• এক মাত্রার তিন-বিন্দুর উদাহরণে posterior হাতে আর কোডে মিলিয়েছি
• Posterior mean-কে Ridge(alpha=α/β, fit_intercept=False)-এর সঙ্গে মিলিয়েছি
• scikit-learn-এর alpha_/lambda_-এর অর্থ নিজের version-এর docs-এ দেখেছি
• Predictive sd-তে noise আর epistemic দুই অংশই আছে, আর ডেটা থেকে দূরে পরিসর চওড়া হচ্ছে
• আলাদা রাখা ডেটায় coverage মেপেছি, কয়েকটি seed-এ, একটি সরল baseline-এর পাশে
• Feature standardize করে prior দিয়েছি, আর ARD-এর ফলকে কারণ-কার্য বলে দাবি করিনি
• আসল ক্রেতার ঠিকানা, ফোন নম্বর বা অর্ডারের ব্যক্তিগত তথ্য AI-তে পেস্ট করিনি