MCMC — Metropolis-Hastings আর Gibbs sampling

AI-এর ভিত্তি কোর্সের পর্ব ৫০ শিখিয়েছে prior, likelihood আর posterior, আর Beta–Binomial-এ posterior থেকে নমুনা নেওয়া। সেখানে posterior-এর সূত্র হাতে ছিল। এই কোর্সের পর্ব ১৮ posterior-এর সবচেয়ে উঁচু বিন্দু (MAP) খুঁজেছে। কিন্তু একটি বিন্দু অনিশ্চয়তার কথা বলে না। এই পর্বের প্রশ্ন: সূত্র হাতে না থাকলেও পুরো posterior থেকে নমুনা কীভাবে পাব? চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python: NumPy, SciPy, matplotlib আর scikit-learn সহ)। প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে, scikit-learn-এ প্রায় ১০ সেকেন্ড।

১ · কেন নমুনা লাগে — posterior-এর integral আর অজানা হর

Bayes-এর সূত্রে posterior হলো

p(θ | D) = p(D | θ)·p(θ) ÷ Z, Z = ∫ p(D | θ)·p(θ) dθ

এখানে θ হলো parameter (যেমন দিনে ফেরতের গড় হার), D ডেটা, p(D | θ) likelihood, p(θ) prior, আর Z হলো normalising constant, যার আরেক নাম evidence. লব হিসাব করা সহজ: যেকোনো θ-তে দুটি সংখ্যা গুণ। হর Z কঠিন, কারণ সব সম্ভাব্য θ-এর ওপর integral লাগে। আমরা আসলে যা চাই, তা-ও integral:

E[f(θ) | D] = ∫ f(θ)·p(θ | D) dθ ≈ (১/S)·∑_s f(θ^(s))

f(θ) = θ নিলে posterior গড়। f হলো “θ > ৪.৫ হলে ১, নইলে ০” নিলে একটি সম্ভাবনা। Percentile থেকে credible interval. ডান দিকের গড় হলো Monte Carlo অনুমান: posterior থেকে S-টি নমুনা থাকলে সব প্রশ্নের উত্তর গড় নিয়ে মেলে, আর স্বাধীন নমুনায় ভুল কমে ১/√S হারে (AI-এর ভিত্তি কোর্সের পর্ব ৩৬)।

Grid কেন চলে না? এক মাত্রায় ১০০টি বিন্দুতে লব হিসাব করে যোগ করলেই Z মেলে। দুই মাত্রায় ১০০² = ১০,০০০টি বিন্দু, দশ মাত্রায় ১০০^১০ = ১০^২০টি। একটি ছোট logistic regression-এও ২০টি feature থাকতে পারে। তাই grid-এর পথ বন্ধ। আবার posterior থেকে সরাসরি স্বাধীন নমুনা নেওয়ার উপায়ও সাধারণত নেই। ব্যতিক্রম conjugate prior: Poisson likelihood আর Gamma prior হলে posterior আবার Gamma. অংশ ৪-এ এই বন্ধ উত্তরটিকেই আমরা মাপকাঠি বানাব।

২ · Markov chain দিয়ে নমুনা — stationary distribution আর detailed balance

AI-এর ভিত্তি কোর্সের পর্ব ৫৩-এ Markov chain এসেছে: পরের অবস্থা শুধু এখনকার অবস্থার ওপর নির্ভর করে, transition matrix P দিয়ে। আর stationary distribution π সেই বণ্টন, যা একটি ধাপের পরেও বদলায় না: πP = π. সাধারণ Markov chain-এ P জানা, π খুঁজতে হয়। MCMC উল্টো কাজ করে: π (posterior) জানা, এমন P বানাতে হবে যার stationary distribution ঠিক π.

সবচেয়ে সহজ শর্তটি হলো detailed balance:

π(i)·P(i → j) = π(j)·P(j → i) প্রতিটি জোড়া i, j-এর জন্য

বাঁ দিক হলো “i-তে থাকা আর i থেকে j-তে যাওয়া”-র প্রবাহ, ডান দিক উল্টো প্রবাহ। প্রতিটি জোড়ায় দুই দিকের প্রবাহ সমান হলে কোনো অবস্থায় সম্ভাবনা জমে না বা কমে না। গাণিতিকভাবে দুই দিকে i-এর ওপর যোগ করলে ∑_i π(i)P(i → j) = π(j)∑_i P(j → i) = π(j), অর্থাৎ πP = π. Detailed balance তাই stationarity-র জন্য যথেষ্ট (তবে আবশ্যক নয়)। আরও দুটি শর্ত লাগে, যাতে যেকোনো শুরু থেকে chain π-তে পৌঁছায়: সব অবস্থায় পৌঁছানো যায় (irreducible), আর chain একটি নির্দিষ্ট ছন্দে আটকে থাকে না (aperiodic). নিচের stepper দুই অবস্থার ছোট একটি chain-এ সব হিসাব হাতে করে।

৩ · Metropolis-Hastings — প্রস্তাব, গ্রহণ আর বাতিল

অবিচ্ছিন্ন θ-তে একই কৌশল। ধরা যাক π̃(θ) = p(D | θ)p(θ), অর্থাৎ Z ছাড়া posterior. বর্তমান মান θ থেকে প্রতিটি ধাপে:

1. একটি proposal distribution q(θ′ | θ) থেকে নতুন মান θ′ আঁকা।
2. গ্রহণের সম্ভাবনা হিসাব: α = min(১, π̃(θ′)·q(θ | θ′) ÷ [π̃(θ)·q(θ′ | θ)]).
3. ০ থেকে ১-এর একটি uniform সংখ্যা u আঁকা। u < α হলে θ′ গ্রহণ, নইলে θ-ই থাকে। দুই ক্ষেত্রেই chain-এ এই ধাপের মান লেখা হয়।

অনুপাতে Z উপরে-নিচে থাকে, তাই কেটে যায়। এই কারণেই MCMC অজানা হরের posterior-এ কাজ করে। q প্রতিসম হলে (যেমন θ′ = θ + σ_prop·z, z ~ N(০, ১)) q-এর অংশও কাটে, আর থাকে শুধু π̃(θ′)/π̃(θ). এটিই random-walk Metropolis. প্রস্তাব অসমমিত হলে q-এর অনুপাতটি হলো Hastings-এর সংশোধন। যেমন ধনাত্মক λ-এ গুণ আকারের লাফ (λ′ = λ·e^σz) দিলে সংশোধনটি λ′/λ. Metropolis-Hastings কেন detailed balance মানে, এক লাইনে দেখা যায়: θ ≠ θ′ হলে π(θ)q(θ′ | θ)α(θ → θ′) = min(π(θ)q(θ′ | θ), π(θ′)q(θ | θ′)), যা θ আর θ′-এ প্রতিসম।

এবার অংশ ৪-এর আসল লক্ষ্যে এক-দুটি ধাপ হাতে করা যাক। সেখানে log π̃(λ) = ১৪৮·ln λ − ৩০.৫·λ (সংখ্যাগুলো কোথা থেকে, অংশ ৪-এ আসবে)। বাস্তব কোডে সবসময় log-এ কাজ হয়, কারণ likelihood-এর গুণফল সহজেই underflow করে।

৪ · এক মাত্রায় হাতে-কলমে — ফেরতের হার, Gamma prior আর আসল উত্তর

বাজারঘরে ৩০ দিনে প্রতিদিন কতগুলো অর্ডার ফেরত এসেছে, সেই কাল্পনিক গণনা নিন। Model: প্রতিদিনের গণনা y_i ~ Poisson(λ), prior λ ~ Gamma(shape ২, rate ০.৫), যার গড় ৪. Gamma–Poisson conjugate, তাই posterior Gamma(২ + ∑y, ০.৫ + n). এখানে ∑y = ১৪৭, n = ৩০, তাই posterior Gamma(১৪৯, ৩০.৫), আর log π̃(λ) = ১৪৮·ln λ − ৩০.৫·λ. নিচের sampler এই সূত্র ব্যবহার করে না। সে শুধু log π̃ জানে; শেষে আসল উত্তরের সঙ্গে মেলানো হয়।

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng(47)            # কাল্পনিক বাজারঘর: ৩০ দিনের ফেরত অর্ডার
y = rng.poisson(4.2, size=30)
a0, b0 = 2.0, 0.5                          # Gamma(shape 2, rate 0.5) prior, গড় 4
S, n = y.sum(), len(y)
def log_post(lam):                         # Z ছাড়া log posterior
    return (a0 - 1 + S) * np.log(lam) - (b0 + n) * lam if lam > 0 else -np.inf
lam, step, draws, acc = 1.0, 0.6, [], 0    # শুরু ইচ্ছা করে দূরে
for t in range(6000):
    prop = lam + step * rng.normal()       # random-walk প্রস্তাব
    if np.log(rng.uniform()) < log_post(prop) - log_post(lam):
        lam, acc = prop, acc + 1           # গ্রহণ
    draws.append(lam)                      # বাতিল হলে পুরোনো মানই আবার
draws = np.array(draws[1000:])             # প্রথম 1000 বাদ (warm-up)
exact = stats.gamma(a0 + S, scale=1 / (b0 + n))   # conjugate উত্তর
print("returns per day:", y[:10], "... sum", S, "n", n)
print(f"acceptance rate: {acc / 6000:.2f}")
print(f"posterior mean  MCMC {draws.mean():.3f}  exact {exact.mean():.3f}")
print(f"posterior sd    MCMC {draws.std():.3f}  exact {exact.std():.3f}")
print("95% interval MCMC", np.percentile(draws, [2.5, 97.5]).round(2),
      " exact", np.array(exact.interval(0.95)).round(2))
print(f"P(lambda > 4.5): MCMC {np.mean(draws > 4.5):.3f}  exact {exact.sf(4.5):.3f}")
g = np.linspace(2.5, 6, 300)
fig, ax = plt.subplots(figsize=(6, 3.4), layout="constrained")
ax.hist(draws, bins=40, density=True, alpha=0.5, label="Metropolis draws")
ax.plot(g, exact.pdf(g), "k", label="exact Gamma posterior")
ax.set_xlabel("lambda (returns per day)"); ax.set_ylabel("density")
ax.set_title("Random-walk Metropolis vs conjugate posterior"); ax.legend(fontsize=8)
plt.show()

returns per day: [ 5  4  8  3 11  6  2  3  6  4] ... sum 147 n 30
acceptance rate: 0.60
posterior mean  MCMC 4.876  exact 4.885
posterior sd    MCMC 0.396  exact 0.400
95% interval MCMC [4.14 5.67]  exact [4.13 5.7 ]
P(lambda > 4.5): MCMC 0.824  exact 0.832

৫ · লাফের মাপ, autocorrelation আর effective sample size

Proposal-এর মাপ σ_prop একমাত্র “নব”, আর এটি দুই দিকেই ভুল হতে পারে। খুব ছোট লাফ প্রায় সবই গ্রহণ হয়, কিন্তু chain পিঁপড়ের মতো হাঁটে। খুব বড় লাফ প্রায় সবই লেজে পড়ে বাতিল হয়, chain একই জায়গায় দাঁড়িয়ে থাকে। মাঝের মাপটি খোঁজার জন্য দুটি মাপকাঠি লাগে।

Autocorrelation ρ_k হলো chain-এর মান আর k ধাপ পরের মানের correlation. ρ_k ধীরে কমলে পাশাপাশি নমুনাগুলো প্রায় একই তথ্য বহন করে। Effective sample size:

ESS = N ÷ (১ + ২∑_{k ≥ ১} ρ_k)

হরের রাশিটিকে বলে integrated autocorrelation time τ. স্বাধীন নমুনায় সব ρ_k = ০, তাই ESS = N. যদি ρ_k = ০.৮^k হয়, তাহলে ∑ρ_k = ০.৮/০.২ = ৪, τ = ৯, আর ১০,০০০টি নমুনা কাজ করে প্রায় ১,১১১টি স্বাধীন নমুনার মতো। গড়ের Monte Carlo ভুল তাই MCSE ≈ sd ÷ √ESS, √N নয়। বাস্তবে দূরের lag-এর ρ_k শব্দে ভরা, তাই যোগফল কোথাও কেটে দিতে হয়। নিচের ব্লক সবচেয়ে সরল নিয়ম নেয়: প্রথম ঋণাত্মক ρ_k পর্যন্ত যোগ। ArviZ-এর মতো টুল Geyer-এর আরও সাবধানী নিয়ম আর কয়েকটি chain একসঙ্গে ব্যবহার করে, তাই তাদের সংখ্যা একটু আলাদা হবে।

import math
import numpy as np
import matplotlib.pyplot as plt
A, B = 149.0, 30.5                         # Gamma(149, rate 30.5): আগের ব্লকের posterior
def log_post(lam):
    return (A - 1) * math.log(lam) - B * lam if lam > 0 else -math.inf
def metropolis(step, n=5000, lam=4.0, seed=1):
    rng = np.random.default_rng(seed)
    z, u, out, acc = rng.normal(size=n), rng.uniform(size=n), np.empty(n), 0
    for t in range(n):
        prop = lam + step * z[t]
        if math.log(u[t]) < log_post(prop) - log_post(lam):
            lam, acc = prop, acc + 1
        out[t] = lam
    return out, acc / n
def acf(x):                                # autocorrelation, FFT দিয়ে
    x = x - x.mean(); f = np.fft.rfft(x, 2 * len(x))
    r = np.fft.irfft(f * np.conj(f))[:len(x)]
    return r / r[0]
def ess(x):                                # প্রথম ঋণাত্মক lag পর্যন্ত যোগ
    rho = acf(x); k = np.argmax(rho < 0)
    return len(x) / (1 + 2 * rho[1:k].sum())
fig, axes = plt.subplots(3, 1, figsize=(6, 5), sharex=True, layout="constrained")
for ax, step in zip(axes, [0.03, 1.0, 8.0]):
    d, rate = metropolis(step)
    rho = acf(d)
    print(f"step {step:5}: accept {rate:.2f}  lag1 {rho[1]:.3f}  lag10 {rho[10]:.3f}"
          f"  ESS {ess(d):6.0f} of {len(d)}  mean {d.mean():.3f}")
    ax.plot(d, lw=0.6); ax.set_ylabel(f"step {step}")
axes[0].set_title("Trace plots: proposal too small / about right / too large")
axes[2].set_xlabel("iteration"); plt.show()
d, _ = metropolis(1.0, n=20000)
for k in (1, 5, 10):                       # thinning: প্রতি k-তম রাখা
    print(f"thin {k:2d}: kept {len(d[::k]):5d}  ESS {ess(d[::k]):6.0f}")

step  0.03: accept 0.97  lag1 0.998  lag10 0.977  ESS      6 of 5000  mean 4.734
step   1.0: accept 0.43  lag1 0.608  lag10 0.000  ESS   1125 of 5000  mean 4.864
step   8.0: accept 0.07  lag1 0.916  lag10 0.444  ESS    175 of 5000  mean 4.876
thin  1: kept 20000  ESS   4969
thin  5: kept  4000  ESS   3442
thin 10: kept  2000  ESS   2000

কোন গ্রহণের হার লক্ষ্য? Gaussian লক্ষ্যে random-walk Metropolis-এর তাত্ত্বিক ফল: এক মাত্রায় প্রায় ০.৪৪, আর অনেক মাত্রায় প্রায় ০.২৩৪ (Gelman, Roberts আর Gilks ১৯৯৬; Roberts, Gelman আর Gilks ১৯৯৭)। এগুলো আদর্শ অবস্থার নির্দেশিকা, আইন নয়। এক মাত্রায় ভালো লাফের মাপ মোটামুটি ২.৪ × posterior sd. এখানে ২.৪ × ০.৪ ≈ ১.০, ঠিক ব্লকের মাঝের মাপ। তাই অনেক sampler warm-up-এ গ্রহণের হার দেখে লাফের মাপ নিজেই মানিয়ে নেয়। Warm-up শেষ হলে মাপ স্থির করতে হয়; নইলে chain আর Markov chain থাকে না, আর নিশ্চয়তাগুলো ভাঙে।

৬ · কয়েকটি chain, warm-up আর R-hat — দুই মাত্রার logistic posterior

একটি chain কখন “যথেষ্ট” চলেছে, তা সেই chain দেখে বলা কঠিন। সমাধান: কয়েকটি chain, একে অপরের থেকে দূরের শুরু থেকে। সবগুলো যদি একই বণ্টন দেখায়, তবেই ভরসা। এই তুলনার মাপ R-hat (Gelman আর Rubin ১৯৯২)। m-টি chain, প্রতিটিতে n-টি মান। W = প্রতিটি chain-এর ভেতরের variance-এর গড়; B/n = chain-গুলোর গড়ের variance. তারপর

V̂ = ((n − ১)/n)·W + B/n, R̂ = √(V̂ ÷ W)

Chain-গুলো একই বণ্টনে থাকলে তাদের গড় প্রায় এক, B/n ছোট, আর R̂ ≈ ১. কোনো chain অন্য জায়গায় আটকে থাকলে B/n বড়, R̂ > ১. Vehtari ও সহকর্মীরা (২০২১) rank-normalised, split রূপ প্রস্তাব করেন, আর R̂ < ১.০১ ও ESS > ৪০০-কে নির্দেশিকা বলেন; ArviZ সেই রূপটি দেয়। নিচের ব্লক সহজ ক্লাসিক রূপটি নিজে হিসাব করে।

Model: বাজারঘরের ২০০টি কাল্পনিক অর্ডার, অর্ডারের মূল্য x (হাজার টাকায়), ফেরত হলে y = ১. Logistic regression: P(y = ১) = σ(a + b·x), prior a, b ~ N(০, ৫²)। চারটি chain চার কোণ থেকে শুরু, ৫,০০০ ধাপ, প্রথম ১,০০০টি warm-up. চারটি chain একসঙ্গে array-তে চলে, তাই loop একটিই। শেষে posterior গড় মেলানো হয় scikit-learn-এর penalty ছাড়া fit-এর (C = 1e10) MLE-র সঙ্গে।

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
rng = np.random.default_rng(7)             # কাল্পনিক বাজারঘর: অর্ডার মূল্য (হাজার টাকা), ফেরত?
x = rng.uniform(0.3, 5.0, 200)
y = (rng.uniform(size=200) < 1 / (1 + np.exp(-(-2.5 + 0.6 * x)))).astype(float)
X = np.column_stack([np.ones(200), x])
def log_post(th):                          # th: (chain, 2); prior N(0, 5²) দুটিতেই
    eta = X @ th.T                         # 200 × chain
    ll = (y[:, None] * eta - np.logaddexp(0, eta)).sum(axis=0)
    return ll - (th ** 2).sum(axis=1) / (2 * 25)
th = np.array([[-6.0, 2.0], [3.0, -1.5], [-6.0, -1.5], [3.0, 2.0]])  # ছড়ানো শুরু
lp, T, keep, acc = log_post(th), 5000, [], 0
for t in range(T):
    prop = th + rng.normal(size=th.shape) * [0.5, 0.15]
    lp_new = log_post(prop)
    ok = np.log(rng.uniform(size=4)) < lp_new - lp       # প্রতিটি chain আলাদা সিদ্ধান্ত
    th[ok], lp[ok] = prop[ok], lp_new[ok]
    keep.append(th.copy()); acc += ok.sum()
draws = np.array(keep)                     # (T, 4 chain, 2)
def rhat(c):                               # c: (n, m chain); Gelman-Rubin
    n = len(c); W = c.var(axis=0, ddof=1).mean()
    Bn = c.mean(axis=0).var(ddof=1)        # B/n
    return np.sqrt(((n - 1) / n * W + Bn) / W)
for name, d in [("first 100 (warm-up)", draws[:100]), ("after warm-up", draws[1000:])]:
    print(f"R-hat {name:20s}: a {rhat(d[:, :, 0]):.3f}  b {rhat(d[:, :, 1]):.3f}")
post = draws[1000:].reshape(-1, 2); mle = LogisticRegression(C=1e10).fit(x[:, None], y)
print(f"acceptance rate: {acc / (4 * T):.2f}")
print("posterior mean (a, b):", post.mean(axis=0).round(3), " sd:", post.std(axis=0).round(3))
print("MLE sklearn    (a, b):", np.round([mle.intercept_[0], mle.coef_[0, 0]], 3))
print("95% interval b:", np.percentile(post[:, 1], [2.5, 97.5]).round(3))
fig, ax = plt.subplots(figsize=(6, 3.2), layout="constrained")
ax.plot(draws[:, :, 1], lw=0.5); ax.axvspan(0, 1000, color="grey", alpha=0.15)
ax.set_xlabel("iteration (shaded: warm-up)"); ax.set_ylabel("slope b")
ax.set_title("Four chains from dispersed starts"); plt.show()

R-hat first 100 (warm-up) : a 1.317  b 1.548
R-hat after warm-up       : a 1.001  b 1.001
acceptance rate: 0.24
posterior mean (a, b): [-2.592  0.621]  sd: [0.431 0.131]
MLE sklearn    (a, b): [-2.563  0.613]
95% interval b: [0.377 0.892]

Warm-up বনাম burn-in. “Burn-in” পুরোনো নাম: শুরুর অংশ ফেলে দেওয়া। আধুনিক টুলে এই সময়ে লাফের মাপ আর দিকও মানিয়ে নেওয়া হয়, তাই নাম warm-up বা tuning. কতটা বাদ দেবেন, তার কোনো জাদু সংখ্যা নেই। Trace plot আর R-hat দেখে ঠিক করুন। নিচের stepper অংশ ৫ আর ৬-এর আসল chain থেকে বানানো ছোট sparkline দিয়ে trace plot পড়া শেখায়।

৭ · Gibbs sampling — একবারে একটি চলক

অনেক model-এ পুরো posterior জটিল, কিন্তু বাকি সব স্থির রাখলে একটি parameter-এর বণ্টন পরিচিত। এই বণ্টনকে বলে full conditional. Gibbs sampling পালাক্রমে প্রতিটি parameter তার full conditional থেকে আঁকে। এটি Metropolis-Hastings-এরই বিশেষ রূপ: প্রস্তাব যখন ঠিক full conditional, তখন α = ১, তাই কোনো বাতিল নেই আর লাফের মাপও বাছতে হয় না।

বাজারঘরের ডেলিভারি সময়। ৪০টি কাল্পনিক ডেলিভারির সময় y_i (ঘণ্টা), model y_i ~ N(μ, σ²)। Prior: μ ~ N(৪৮, ২০²), আর σ² ~ Inverse-Gamma(২, ২০০)। দুটি prior আলাদা, তাই যুক্ত posterior-এর বন্ধ রূপ সহজ নয়। কিন্তু full conditional দুটি পরিচিত:

ব্লকের দ্বিতীয় অংশ Gibbs-এর দুর্বলতা দেখায়। লক্ষ্য একটি দুই মাত্রার standard Gaussian, correlation ρ. এখানে a | b ~ N(ρb, ১ − ρ²), আর b | a একইভাবে। ρ = ০.৯৯ হলে প্রতিটি শর্তসাপেক্ষ বণ্টনের sd মাত্র √(১ − ০.৯৮০১) ≈ ০.১৪, অথচ লক্ষ্যের ছড়ানো ১. Chain তাই সরু উপত্যকায় অক্ষ-বরাবর ছোট ছোট সিঁড়ি-ধাপে চলে। তত্ত্ব বলে, a-এর lag-১ autocorrelation ঠিক ρ².

import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(11)            # কাল্পনিক বাজারঘর: ৪০টি ডেলিভারির সময় (ঘণ্টা)
y = rng.normal(52, 14, 40); n, ybar = len(y), y.mean()
m0, s0sq, a0, b0 = 48.0, 20.0 ** 2, 2.0, 200.0   # μ ~ N(48, 20²), σ² ~ Inv-Gamma(2, 200)
mu, sig2, out = 30.0, 100.0, []
for t in range(4000):
    v = 1 / (1 / s0sq + n / sig2)                       # μ | σ², y : Normal
    mu = rng.normal(v * (m0 / s0sq + n * ybar / sig2), np.sqrt(v))
    rate = b0 + 0.5 * np.sum((y - mu) ** 2)             # σ² | μ, y : Inv-Gamma
    sig2 = rate / rng.gamma(a0 + n / 2)                 # 1/Gamma = Inv-Gamma
    out.append((mu, np.sqrt(sig2)))
d = np.array(out[500:])
print(f"data: mean {ybar:.2f} h, sd {y.std(ddof=1):.2f} h")
print("posterior mean (mu, sigma):", d.mean(axis=0).round(2))
print("95% interval mu:", np.percentile(d[:, 0], [2.5, 97.5]).round(2))
def gibbs_2d(rho, T=2000):                 # bivariate Normal, corr rho
    a, b, path = 3.0, -3.0, []
    for t in range(T):
        a = rng.normal(rho * b, np.sqrt(1 - rho ** 2))  # a | b
        b = rng.normal(rho * a, np.sqrt(1 - rho ** 2))  # b | a
        path.append((a, b))
    return np.array(path)
fig, axes = plt.subplots(1, 2, figsize=(6.4, 3.2), layout="constrained")
for ax, rho in zip(axes, [0.5, 0.99]):
    p = gibbs_2d(rho); lag1 = np.corrcoef(p[:-1, 0], p[1:, 0])[0, 1]
    print(f"rho {rho}: lag-1 autocorr of a {lag1:.3f}, mean of a {p[:, 0].mean():.3f}")
    ax.plot(p[:40, 0], p[:40, 1], "o-", ms=2, lw=0.7)
    ax.set_xlim(-4, 4); ax.set_ylim(-4, 4); ax.set_title(f"Gibbs, rho = {rho}, 40 steps")
    ax.set_xlabel("a"); ax.set_ylabel("b")
plt.show()

data: mean 50.85 h, sd 12.00 h
posterior mean (mu, sigma): [50.8  12.03]
95% interval mu: [47.05 54.39]
rho 0.5: lag-1 autocorr of a 0.204, mean of a -0.010
rho 0.99: lag-1 autocorr of a 0.981, mean of a -0.177

প্রতিকার। Correlated চলকগুলো একসঙ্গে আঁকা (blocked Gibbs), parameter বদলে নেওয়া (যেমন x-কে centre করলে intercept আর slope-এর correlation কমে), বা gradient-ভিত্তিক পদ্ধতি (অংশ ৮)। Topic model-এর সাহিত্যে collapsed Gibbs sampling খুব পরিচিত (Griffiths আর Steyvers ২০০৪); পর্ব ৪৪-এ ব্যবহার করা scikit-learn-এর LDA অবশ্য variational পদ্ধতিতে চলে, যার ধারণা পর্ব ৪৮ (Variational inference-এর ধারণা — ELBO আর approximation)-এ আসবে।

৮ · HMC আর NUTS — ঢাল দেখে লম্বা লাফ, আর PyMC-র একটি সত্যিকারের run

Random walk-এর সমস্যা মাত্রার সঙ্গে বাড়ে। d মাত্রায় গ্রহণের হার ধরে রাখতে লাফ ছোট করতে হয়, আর এক প্রান্ত থেকে অন্য প্রান্তে যেতে মোটামুটি d-এর অনুপাতে বেশি ধাপ লাগে। Hamiltonian Monte Carlo (HMC) একটি ভৌত তুলনা ব্যবহার করে। θ-কে একটি কণার অবস্থান ধরুন, −log π̃(θ)-কে পাহাড়ি উপত্যকার উচ্চতা। কণাকে একটি এলোমেলো ধাক্কা (momentum) দিয়ে কিছুক্ষণ গড়াতে দিন। Gradient কণাকে উপত্যকার আকার মেনে বাঁকায়, তাই সে সরু, হেলানো উপত্যকা ধরেও অনেক দূর যায়। গতিপথ সংখ্যাগতভাবে (leapfrog ধাপে) হিসাব হয়। ভুল জমে, তাই শেষে Metropolis-এর একটি গ্রহণ-পরীক্ষা দিয়ে সংশোধন হয়। ফল: দূরের প্রস্তাব, তবু উচ্চ গ্রহণের হার।

HMC-তে দুটি নব আছে: leapfrog-এর ধাপের মাপ আর ধাপের সংখ্যা। Hoffman আর Gelman-এর No-U-Turn Sampler (NUTS; প্রবন্ধ ২০১৪ সালে Journal of Machine Learning Research-এ) দ্বিতীয়টি নিজে ঠিক করে: গতিপথ যখন নিজের দিকে ফিরতে শুরু করে, তখন থামে। প্রথমটি warm-up-এ মানিয়ে নেয়। Stan আর PyMC-র ডিফল্ট sampler এটিই। একটি সীমা মনে রাখুন: gradient লাগে বলে বিচ্ছিন্ন (discrete) parameter, যেমন mixture-এর দল-নম্বর, সরাসরি চলে না; সেগুলো যোগ করে বাদ দিতে (marginalise) হয়।

নিচে অংশ ৬-এর একই ডেটা আর একই prior, এবার PyMC দিয়ে। এটি ব্রাউজারে চলে না, তাই একটি আলাদা uv project-এ ডেস্কটপে চালানো হয়েছে (সত্যিকারের run; লেখার সময়ের সর্বশেষ stable version). একে বলে probabilistic programming: model লিখে দিলে inference টুল নিজে চালায়।

~$ uv init --vcs none --author-from none bazarghar-ml
~$ cd bazarghar-ml && uv add pymc
 + pymc==6.4.0
 + pytensor==3.3.7
 …
~/bazarghar-ml$ cat nuts_pymc.py

import sys
import numpy as np
import pymc as pm
import arviz as az
import pandas as pd
pd.set_option("display.width", 120); pd.set_option("display.max_columns", None)

rng = np.random.default_rng(7)             # same fictional data as the in-page block
x = rng.uniform(0.3, 5.0, 200)
y = (rng.uniform(size=200) < 1 / (1 + np.exp(-(-2.5 + 0.6 * x)))).astype(int)

print("python", sys.version.split()[0], "| pymc", pm.__version__, "| arviz", az.__version__)
with pm.Model() as model:
    a = pm.Normal("a", 0, 5)
    b = pm.Normal("b", 0, 5)
    pm.Bernoulli("returned", logit_p=a + b * x, observed=y)
    idata = pm.sample(draws=1000, tune=1000, chains=4, cores=1,
                      random_seed=47, progressbar=False)
print(az.summary(idata, var_names=["a", "b"], round_to=3))
print("divergences:", int(idata.sample_stats["diverging"].sum()))

~/bazarghar-ml$ uv run python nuts_pymc.py
g++ not available, if using conda: `conda install gxx`
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (4 chains in 1 job)
NUTS: [a, b]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.
python 3.14.4 | pymc 6.4.0 | arviz 1.3.0
    mean    sd  eti89_lb  eti89_ub  ess_bulk  ess_tail  r_hat  mcse_mean  mcse_sd
a -2.569  0.43    -3.262    -1.879   915.576   888.047  1.010      0.014    0.010
b  0.613  0.13     0.407     0.817   918.595   929.604  1.009      0.004    0.003
divergences: 0

পড়া যাক। প্রথম লাইন একটি সতর্কবার্তা: এই Windows যন্ত্রে C++ compiler (g++) নেই, তাই PyMC-র ভেতরের PyTensor দ্রুততম পথে চলে না। তবু ১,০০০ tune আর ১,০০০ draw-এর চারটি chain শেষ হয়েছে ২ সেকেন্ডে। গড় (−২.৫৬৯, ০.৬১৩), sd (০.৪৩০, ০.১৩০), নিজের sampler-এর (−২.৫৯২, ০.৬২১) আর (০.৪৩১, ০.১৩১)-এর খুব কাছে। ArviZ 1.x ডিফল্টে ৮৯% equal-tailed interval (eti89) দেখায়, ৯৫% নয়। r_hat ১.০১০ আর ১.০০৯, নির্দেশিকার একেবারে সীমায়; আরও draw নিলে কমার কথা। ESS: ৪,০০০টি draw থেকে bulk ESS প্রায় ৯১৬ আর ৯১৯. একই ডেস্কটপে ArviZ-এর az.ess দিয়ে অংশ ৬-এর random-walk chain-গুলো মাপলে ১৬,০০০টি draw থেকে প্রায় ৫৮৬ (এক local run). প্রতি draw-এ NUTS প্রায় ছয় গুণ বেশি তথ্য দিয়েছে, যদিও প্রতি draw-এ সে কয়েকবার gradient হিসাব করে। এখানে posterior-এ a আর b-এর correlation প্রায় −০.৯২, তাই NUTS-ও পুরো দক্ষ নয়; x-কে centre করলে দুটিরই উন্নতি হওয়ার কথা। Divergence ০, অর্থাৎ গতিপথ কোথাও ভেঙে পড়েনি।

৯ · AI-এর যুগে — MCMC কোড আর ফল যাচাই

AI assistant দ্রুত sampler লিখে দেয়, কিন্তু MCMC-র ভুল প্রায়ই কোনো error ছাড়াই ভুল সংখ্যা দেয়:

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

• একটি conjugate বা আগে জানা উত্তরের উদাহরণে sampler-এর গড়, sd আর পরিসর মিলিয়েছি
• সব হিসাব log-এ হয়েছে, আর বাতিল ধাপে পুরোনো মানটি আবার লেখা হয়েছে, তা কোডে দেখেছি
• অন্তত ৪টি ছড়ানো শুরু থেকে chain চালিয়েছি; R-hat, ESS আর trace plot দেখেছি
• প্রতিটি টুলের parameter-এর নাম আর ডিফল্ট নিজের version-এর docs-এ মিলিয়েছি
• Posterior-এর সংখ্যা একটি সরল baseline-এর (যেমন MLE বা bootstrap) সঙ্গে তুলনা করেছি, আর prior বদলে ফল কতটা নড়ে তা দেখেছি
• কোনো আসল গ্রাহকের ব্যক্তিগত ডেটা AI-তে পেস্ট করিনি