Maximum likelihood আর MAP — model fit করার মূলনীতি

AI-এর ভিত্তি কোর্সের পর্ব ৩৬-এ Bernoulli, Poisson আর normal distribution চিনেছেন, পর্ব ৩৮-এর “আরও গভীরে” বাক্সে দেখেছেন cross-entropy কমানো আর likelihood বাড়ানো একই কথা, আর পর্ব ৫০-এ Beta prior দিয়ে ফেরতের হারের MLE আর MAP পাশাপাশি রেখেছেন। এই কোর্সের পর্ব ০২-এ loss function বেছেছেন, আর পর্ব ১৭-এ শাস্তিকে prior হিসেবে পড়েছেন। এই পর্ব সেই টুকরোগুলোকে একটি নীতির নিচে আনে: model fit করা মানে প্রায় সবসময় কোনো likelihood-কে সর্বোচ্চ করা, কখনো একটি prior-সহ। হাতিয়ার আগের মডিউলের: পর্ব ১১-এর gradient আর Hessian. চালানো যায় এমন চারটি ব্লক পাতাতেই চলে (ব্রাউজারের Python, NumPy, SciPy আর matplotlib); প্রথমবার প্যাকেজ নামতে কয়েক সেকেন্ড লাগে।

১ · Likelihood — ডেটা স্থির, parameter বদলায়

ধরুন x₁, …, x_n ডেটা, আর model বলে প্রতিটি x_i এসেছে একটি distribution p(x | θ) থেকে, যেখানে θ (থিটা) হলো অজানা parameter. পর্যবেক্ষণগুলো স্বাধীন আর একই distribution থেকে এলে (iid), পুরো ডেটার সম্ভাবনা গুণফল। একে θ-এর function হিসেবে পড়লে তার নাম likelihood:

L(θ) = ∏_i=1..n p(x_i | θ)     ℓ(θ) = log L(θ) = ∑_i=1..n log p(x_i | θ)

সূত্রটি সম্ভাবনার মতো দেখায়, কিন্তু পড়ার দিক আলাদা। সম্ভাবনায় θ স্থির আর x বদলায়, তাই সব x-এর ওপর যোগফল ১. Likelihood-এ x স্থির আর θ বদলায়, আর θ-এর ওপর L(θ)-এর যোগফল বা সমাকলন ১ হওয়ার কোনো বাধ্যবাধকতা নেই। তাই “θ = ০.৪ হওয়ার likelihood ০.০৩৫” মানে এই নয় যে θ = ০.৪ হওয়ার সম্ভাবনা ০.০৩৫; সেই দাবির জন্য prior লাগে (অংশ ৭)। Maximum likelihood estimate (MLE) হলো θ̂ = argmax_θ L(θ), অর্থাৎ যে θ দেখা ডেটাকে সবচেয়ে প্রত্যাশিত করে।

Log নেওয়ার তিনটি কারণ। প্রথম, log একটানা বাড়ে (monotone), তাই L আর ℓ-এর সর্বোচ্চ বিন্দু একই θ. দ্বিতীয়, গুণফল যোগফলে বদলায়, আর যোগফলের derivative নেওয়া সহজ। তৃতীয়টি ব্যবহারিক: অনেক ছোট সম্ভাবনার গুণফল দ্রুত float-এর সবচেয়ে ছোট মানের নিচে নেমে চুপচাপ ০ হয়ে যায় (underflow); অংশ ২-এর ব্লকে ২০০০টি সংখ্যার গুণফল ঠিক এভাবেই ০.০ হবে, যেখানে log-এর যোগফল একটি সাধারণ সংখ্যা। Optimizer কমানোর কাজে অভ্যস্ত, তাই বাস্তবে আমরা negative log-likelihood (NLL) −ℓ(θ) কমাই।

সবচেয়ে ছোট উদাহরণ: বাজারঘরের একটি নতুন পণ্যের পাঁচটি অর্ডার, ১ মানে ফেরত এসেছে, ০ মানে আসেনি। Model: প্রতিটি অর্ডার স্বাধীনভাবে θ সম্ভাবনায় ফেরত আসে, অর্থাৎ Bernoulli distribution, p(x | θ) = θ^x(1 − θ)^1−x. নিচের stepper পদে পদে log-likelihood গড়ে।

২ · MLE হাতে — Bernoulli, Poisson আর exponential

রেসিপি সবসময় একই: ℓ(θ) লিখুন, θ-এর সাপেক্ষে derivative নিন, শূন্য বসিয়ে সমাধান করুন, তারপর দ্বিতীয় derivative ঋণাত্মক কি না দেখে নিশ্চিত হোন যে বিন্দুটি সর্বোচ্চ। ℓ-এর প্রথম derivative-এর আলাদা নাম আছে, score function; “score = ০” সমীকরণকে বলে likelihood equation.

Bernoulli. n-টি অর্ডারে k-টি ফেরত হলে ℓ(θ) = k log θ + (n − k) log(1 − θ). Derivative: k/θ − (n − k)/(1 − θ). শূন্য বসালে k(1 − θ) = (n − k)θ, অর্থাৎ θ̂ = k/n. দ্বিতীয় derivative −k/θ² − (n − k)/(1 − θ)², যা সবসময় ঋণাত্মক, তাই এটি সর্বোচ্চ। (k = ০ বা k = n হলে সর্বোচ্চ সীমানায়, θ̂ = ০ বা ১; অংশ ৮-এ এর বিপদ দেখব।)

Poisson. একটি জেলায় দিনে কতটি অর্ডার আসে, এমন গণনার জন্য Poisson distribution: p(x | λ) = e^−λλ^x/x!, যেখানে λ (ল্যামডা) গড় হার। নিচের stepper পাঁচ দিনের ডেটায় derivative শূন্য বসানোর প্রতিটি ধাপ দেখায়।

Exponential. দুটি অর্ডারের মাঝের সময় বা একজন রাইডারের পরের কাজ পাওয়ার অপেক্ষা প্রায়ই exponential distribution দিয়ে ধরা হয়: p(x | λ) = λe^−λx, x ≥ ০. ℓ(λ) = n log λ − λ∑x_i, derivative n/λ − ∑x_i = ০, তাই λ̂ = n/∑x_i = 1/x̄. পাঁচটি অপেক্ষা ২, ৫, ৩, ৬ আর ৪ মিনিট হলে x̄ = ৪, আর λ̂ = ০.২৫ প্রতি মিনিটে। লক্ষ করুন, গড় অপেক্ষা 1/λ-এর MLE সরাসরি x̄ = ৪: parameter-কে নতুন রূপে লিখলে (λ থেকে 1/λ) MLE-ও একই রূপান্তরে যায়। এটি MLE-এর invariance ধর্ম, যা অংশ ৫-এ আবার আসবে।

নিচের ব্লক বাজারঘরের দুটি কাল্পনিক ডেটায় (৪০টি অর্ডারের ফেরত, আর এক জেলার ৩০ দিনের অর্ডারের সংখ্যা) likelihood আর log-likelihood-এর বক্ররেখা আঁকে, closed-form MLE লাল দাগে চিহ্নিত করে, grid-এ খোঁজা সর্বোচ্চ বিন্দুর সঙ্গে মেলায়, আর underflow-টি নিজের চোখে দেখায়।

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

rng = np.random.default_rng(42)            # কাল্পনিক বাজারঘর
ret = rng.binomial(1, 0.12, size=40)        # ৪০টি অর্ডার: 1 = ফেরত এসেছে
orders = rng.poisson(6.5, size=30)          # এক জেলায় ৩০ দিনের অর্ডারের সংখ্যা
k, n = ret.sum(), ret.size
th = np.linspace(0.005, 0.4, 400)
L_b = th**k * (1 - th)**(n - k)             # likelihood: θ-এর function
ll_b = k * np.log(th) + (n - k) * np.log(1 - th)
lam = np.linspace(4, 10, 400)
ll_p = stats.poisson.logpmf(orders[:, None], lam).sum(axis=0)

print(f"returns: k = {k} of n = {n}; MLE k/n = {k / n:.4f}; grid argmax = {th[ll_b.argmax()]:.4f}")
print(f"orders : sum = {orders.sum()}, n = {orders.size}; MLE mean = {orders.mean():.4f}; grid argmax = {lam[ll_p.argmax()]:.4f}")
print(f"largest likelihood value (returns): {L_b.max():.3e}")
print(f"Poisson likelihood at the MLE     : {np.exp(ll_p.max()):.3e}")
big = rng.poisson(6.5, size=2000)
print("product of 2000 pmf values:", np.prod(stats.poisson.pmf(big, 6.5)),
      "| sum of logpmf:", round(stats.poisson.logpmf(big, 6.5).sum(), 1))

fig, ax = plt.subplots(2, 2, figsize=(7, 5), layout="constrained")
for a, x, y, m, t in [(ax[0, 0], th, L_b, k / n, "Bernoulli likelihood"),
                      (ax[1, 0], th, ll_b, k / n, "Bernoulli log-likelihood"),
                      (ax[0, 1], lam, np.exp(ll_p - ll_p.max()), orders.mean(), "Poisson likelihood (scaled)"),
                      (ax[1, 1], lam, ll_p, orders.mean(), "Poisson log-likelihood")]:
    a.plot(x, y, lw=2); a.axvline(m, color="tab:red", ls="--", label=f"MLE = {m:.3f}")
    a.set_title(t, fontsize=10); a.legend(fontsize=8)
ax[1, 0].set_xlabel("return rate theta"); ax[1, 1].set_xlabel("orders per day lambda")
plt.show()

returns: k = 5 of n = 40; MLE k/n = 0.1250; grid argmax = 0.1248
orders : sum = 188, n = 30; MLE mean = 6.2667; grid argmax = 6.2707
largest likelihood value (returns): 2.850e-07
Poisson likelihood at the MLE     : 9.929e-30
product of 2000 pmf values: 0.0 | sum of logpmf: -4732.3

৩ · Gaussian — গড়, variance আর n দিয়ে ভাগের bias

ঢাকার ভেতরের ডেলিভারির সময় (ঘণ্টায়) একটি ঘণ্টাকৃতি normal distribution N(μ, σ²) দিয়ে ধরা যাক। দুটি parameter, তাই দুটি আংশিক derivative:

ℓ(μ, σ²) = −(n/2) log(2πσ²) − ∑(x_i − μ)² / (2σ²)

μ-এর সাপেক্ষে: ∑(x_i − μ)/σ² = ০, তাই μ̂ = x̄. σ²-এর সাপেক্ষে (σ²-কে একটি চলক v ধরে): −n/(2v) + ∑(x_i − μ)²/(2v²) = ০, তাই σ̂² = (1/n)∑(x_i − x̄)². লক্ষ করুন, μ̂ বের করতে σ² জানতে হয়নি; আর ℓ-এর প্রথম পদ দেখায় কেন σ শূন্যের দিকে যেতে পারে না: log σ² তখন −∞-এর দিকে যায়, আর −(n/2) log σ² পদটি ফুলে উঠলেও দ্বিতীয় পদ আরও দ্রুত −∞ হয়।

AI-এর ভিত্তি কোর্সের পর্ব ৩৭-এ শিখেছেন নমুনা থেকে variance অনুমানে n − 1 দিয়ে ভাগ করা হয় (Bessel’s correction). এখানে দেখা যাচ্ছে MLE ঠিক সেটি করে না: সে n দিয়ে ভাগ করে। কারণ: বিচ্যুতিগুলো মাপা হচ্ছে x̄ থেকে, যা একই ডেটা থেকে এমনভাবে বসানো যে ∑(x_i − x̄)² সম্ভাব্য সবচেয়ে ছোট হয়। আসল μ থেকে মাপলে বর্গের যোগফল গড়ে বড় হতো। হিসাব করলে E[σ̂²] = ((n − 1)/n)·σ², অর্থাৎ গড়ে একটু ছোট। একে বলে bias: অনুমানের গড় আর আসল মানের পার্থক্য, এখানে −σ²/n. n বাড়লে bias শূন্যে যায়, তাই অনুমানটি “biased কিন্তু consistent”. নিচের ব্লক ২০,০০০টি ছোট নমুনায় এটি সংখ্যায় দেখায়, আর একটি অপ্রত্যাশিত দিকও: বর্গ-ভুলের মাপে biased অনুমানটিই জেতে।

import numpy as np

rng = np.random.default_rng(7)             # কাল্পনিক: ঢাকার ভেতরে ডেলিভারির সময় (ঘণ্টা)
mu, sigma = 26.0, 6.0                      # আসল গড় আর sd, তাই আসল variance 36
reps = 20000                               # ২০,০০০টি আলাদা ছোট নমুনা

for n in (2, 5, 20, 100):
    x = rng.normal(mu, sigma, size=(reps, n))
    v_mle = x.var(axis=1, ddof=0)          # MLE: n দিয়ে ভাগ
    v_unb = x.var(axis=1, ddof=1)          # unbiased: n − 1 দিয়ে ভাগ
    print(f"n = {n:3d} | mean of mu_hat {x.mean(axis=1).mean():6.2f}"
          f" | E[MLE var] {v_mle.mean():5.2f} (theory {36 * (n - 1) / n:5.2f})"
          f" | E[unbiased] {v_unb.mean():5.2f}")

x = rng.normal(mu, sigma, size=(reps, 5))  # n = 5-এ বর্গ-ভুল (MSE) তুলনা
for name, d in (("MLE      (ddof=0)", 0), ("unbiased (ddof=1)", 1)):
    v = x.var(axis=1, ddof=d)
    print(f"{name}: bias {v.mean() - 36:+6.2f}  sd {v.std():5.2f}"
          f"  MSE {np.mean((v - 36) ** 2):6.1f}")

n =   2 | mean of mu_hat  25.95 | E[MLE var] 17.70 (theory 18.00) | E[unbiased] 35.39
n =   5 | mean of mu_hat  26.02 | E[MLE var] 28.84 (theory 28.80) | E[unbiased] 36.05
n =  20 | mean of mu_hat  26.01 | E[MLE var] 34.14 (theory 34.20) | E[unbiased] 35.94
n = 100 | mean of mu_hat  26.00 | E[MLE var] 35.63 (theory 35.64) | E[unbiased] 35.99
MLE      (ddof=0): bias  -7.30  sd 20.48  MSE  472.8
unbiased (ddof=1): bias  -0.13  sd 25.60  MSE  655.5

৪ · MLE মানে loss কমানো — cross-entropy, KL আর squared error

গড় NLL লিখুন: −(1/n)∑ log q_θ(x_i), যেখানে q_θ হলো model. ডেটার নিজস্ব বণ্টন p̂ (প্রতিটি দেখা মানে ১/n ওজন) ধরলে এটি ঠিক −∑ p̂(x) log q_θ(x), অর্থাৎ p̂ আর q_θ-এর cross-entropy H(p̂, q_θ). আর H(p̂, q) = H(p̂) + KL(p̂ ‖ q), যেখানে H(p̂) θ-এর ওপর নির্ভর করে না। তাই:

MLE = argmin_θ H(p̂, q_θ) = argmin_θ KL(p̂ ‖ q_θ)

অর্থাৎ MLE ডেটার বণ্টনের সবচেয়ে কাছের model খোঁজে, যেখানে “কাছ” মাপা হয় KL দিয়ে। n বড় হলে p̂ আসল বণ্টন p-এর কাছে যায় (law of large numbers), তাই MLE আসলে KL(p ‖ q_θ) কমায়। Model যদি সঠিক হয়, সেই সর্বনিম্ন KL শূন্য, আর আসল θ-তে পৌঁছায়; model ভুল হলে (misspecification) MLE পৌঁছায় KL-এ “সবচেয়ে কাছের ভুল model”-এ। এই দ্বিতীয় কথাটি ভাবনার প্রশ্ন ৩-এ আবার আসবে।

এবার supervised learning. Regression-এ model বলে y = f_w(x) + ε, ε ~ N(০, σ²). তাহলে একটি সারির NLL = (y − f_w(x))²/(2σ²) + ½ log(2πσ²). σ স্থির ধরলে w-এর জন্য বাকি থাকে শুধু বর্গ-ভুলের যোগফল। অর্থাৎ least squares হলো Gaussian noise ধরে MLE. ঐতিহাসিকভাবে যুক্তিটি উল্টো দিক থেকেও এসেছিল: Adrien-Marie Legendre ১৮০৫ সালে least squares প্রকাশ করেন, আর Carl Friedrich Gauss ১৮০৯ সালে দেখান কোন ভুলের বণ্টন ধরলে গড়ই সবচেয়ে সম্ভাব্য মান হয়; উত্তর ছিল আজকের normal distribution. পুরো টেবিল:

পর্ব ০২-এ loss বেছেছিলেন ব্যবসার প্রশ্ন দেখে; এই টেবিল দেখায় প্রতিটি পরিচিত loss-এর পেছনে একটি লুকানো noise-এর অনুমান আছে। MSE বাছা মানে “ভুলগুলো ঘণ্টাকৃতি আর সব জায়গায় সমান ছড়ানো” মেনে নেওয়া। ডেলিভারির সময়ের মতো ডান দিকে লম্বা লেজের ডেটায় এই অনুমান দুর্বল, তাই সেখানে অন্য বণ্টন (পরের অংশের gamma) বা অন্য loss ভালো। AI-এর ভিত্তি কোর্সের পর্ব ৫৯-এর logistic regression-এর cross-entropy loss ঠিক তৃতীয় সারি, আর Poisson-এর সারি পর্ব ২০ (Exponential family আর GLM — Poisson, binomial আর link function)-এর বিষয়।

৫ · Fisher information — চূড়ার বক্রতা থেকে standard error

MLE একটি সংখ্যা দেয়, কিন্তু সেই সংখ্যা কতটা বিশ্বাসযোগ্য? Stepper ২-এর শেষ ধাপের ইঙ্গিত: log-likelihood-এর চূড়া সরু হলে ডেটা θ-কে শক্তভাবে আটকে রাখে; চ্যাপ্টা হলে অনেক θ প্রায় সমান ভালো। চূড়ার বক্রতা মাপে দ্বিতীয় derivative, আর তার ঋণাত্মক রূপের নাম Fisher information:

I(θ) = −E[ d²ℓ/dθ² ]     J(θ̂) = −d²ℓ/dθ² |_{θ = θ̂}     SE(θ̂) ≈ 1/√J(θ̂)

I(θ) হলো প্রত্যাশিত (expected) information, আর J(θ̂) হলো দেখা ডেটায় হিসাব করা observed information. একাধিক parameter-এ J হয়ে যায় NLL-এর Hessian (পর্ব ১১), আর SE আসে J⁻¹-এর কর্ণের বর্গমূল থেকে। n-টি iid পর্যবেক্ষণে information যোগ হয়: I_n(θ) = n·I₁(θ). তাই SE ∝ 1/√n, যা পর্ব ০৪-এর “mean ± sd/√n”-এর সাধারণ রূপ। কয়েকটি পরিচিত মান, প্রতিটি পর্যবেক্ষণে:

একটি numerical: ৮০টি অর্ডারে ১২টি ফেরত। θ̂ = ০.১৫, SE = √(০.১৫ × ০.৮৫/৮০) = ০.০৩৯৯. বড় n-এ θ̂ মোটামুটি normal, তাই প্রায় ৯৫% confidence interval ০.১৫ ± ১.৯৬ × ০.০৩৯৯ = (০.০৭২, ০.২২৮). এই ধাঁচের পরিসরকে Wald interval বলে; θ̂ শূন্য বা একের কাছে আর n ছোট হলে এটি খারাপ আচরণ করে (যেমন k = ০-এ প্রস্থ শূন্য), তখন অন্য পদ্ধতি লাগে।

MLE-এর তিনটি তাত্ত্বিক গুণ, কথায়

• Consistency: model ঠিক থাকলে আর কিছু নিয়মিত শর্ত মানলে, n বাড়ার সঙ্গে θ̂ আসল θ-এর দিকে যায়। এমন অনুমানকে consistent estimator বলে। অংশ ৩-এর σ̂² biased কিন্তু consistent.
• Asymptotic normality আর efficiency: বড় n-এ θ̂ মোটামুটি N(θ, 1/I_n(θ)) (asymptotic normality), আর Cramér–Rao bound বলে কোনো unbiased অনুমানের variance 1/I_n(θ)-এর নিচে যেতে পারে না। তাই বড় n-এ MLE সর্বোচ্চ দক্ষ। ছোট n-এ এই প্রতিশ্রুতি নেই।
• Invariance: g(θ)-এর MLE হলো g(θ̂). λ̂ = 1/x̄ হলে গড় অপেক্ষার MLE x̄, আর ফেরতের odds θ/(1 − θ)-এর MLE θ̂/(1 − θ̂). Unbiasedness এই ধর্ম মানে না: σ̂² unbiased করলেও σ̂ (বর্গমূল) biased থাকে।

৬ · সূত্র নেই — scipy.optimize দিয়ে numerical MLE

জেলার বাইরের ডেলিভারির সময় ধনাত্মক আর ডান দিকে হেলানো: বেশিরভাগ ২ থেকে ৪ দিনে, কিছু ৮–৯ দিনে। Gaussian এখানে মানায় না (ঋণাত্মক সময় দেয়)। স্বাভাবিক পছন্দ gamma distribution, যার shape k আর scale s; গড় k·s, variance k·s². এর log-likelihood:

ℓ(k, s) = (k − 1)∑log x_i − ∑x_i/s − n·k·log s − n·log Γ(k)

s-এর সাপেক্ষে derivative শূন্য বসালে ŝ = x̄/k, অর্থাৎ k̂·ŝ = x̄ হুবহু। কিন্তু k-এর সমীকরণ log k − ψ(k) = log x̄ − (1/n)∑log x_i, যেখানে ψ হলো digamma function; এর কোনো closed form নেই। তখন সংখ্যা দিয়ে NLL কমাতে হয়, ঠিক পর্ব ১৫-এর BFGS বা Newton দিয়ে। ভালো অভ্যাসগুলো:

• NLL কমান, likelihood বাড়াবেন না: minimize কমায়, আর log-এর যোগফল underflow এড়ায়। stats.gamma.logpdf-এর মতো log রূপ সরাসরি ব্যবহার করুন, np.log(pdf) নয়।
• সীমাবদ্ধতা সরান: k, s > ০; optimizer-কে log k আর log s দিন, তাহলে যেকোনো বাস্তব সংখ্যাই বৈধ (reparameterization). বিকল্প হলো bounds-সহ L-BFGS-B (পর্ব ১৬)।
• ভালো শুরুর বিন্দু: method of moments: গড় আর variance মিলিয়ে k₀ = x̄²/v, s₀ = v/x̄.
• যাচাই: res.success দেখুন, library-র fit-এর সঙ্গে মেলান (stats.gamma.fit-এ floc=0 না দিলে সে তৃতীয় একটি location parameter-ও fit করে), আর SE-র জন্য সর্বোচ্চ বিন্দুতে Hessian নিজে হিসাব করুন। BFGS-এর res.hess_inv শুধু optimizer-এর চলার পথে জমানো একটি আনুমানিক মান, আর এখানে তা log-parameter-এর জগতে।

import numpy as np
from scipy import stats, optimize
np.set_printoptions(precision=4, suppress=True)

rng = np.random.default_rng(11)            # কাল্পনিক: জেলার বাইরে ৬০টি ডেলিভারির সময় (দিন)
x = rng.gamma(shape=2.5, scale=1.2, size=60)

def nll(p, x):                             # p = (log k, log s): ধনাত্মক রাখার কৌশল
    return -stats.gamma.logpdf(x, a=np.exp(p[0]), scale=np.exp(p[1])).sum()

def fit(x):
    m, v = x.mean(), x.var()               # শুরু: method of moments
    res = optimize.minimize(nll, np.log([m * m / v, v / m]), args=(x,), method="BFGS")
    return np.exp(res.x), res

theta, res = fit(x)
print("MLE (k, scale):", theta, "| converged:", res.success, "| NLL:", round(res.fun, 4))
print("scipy gamma.fit:", np.array(stats.gamma.fit(x, floc=0))[[0, 2]])

def f(t): return -stats.gamma.logpdf(x, a=t[0], scale=t[1]).sum()
h = 1e-4 * theta; H = np.zeros((2, 2))     # (k, scale)-তে NLL-এর Hessian, central difference
for i in range(2):
    for j in range(2):
        ei, ej = np.eye(2)[i] * h[i], np.eye(2)[j] * h[j]
        H[i, j] = (f(theta + ei + ej) - f(theta + ei - ej)
                   - f(theta - ei + ej) + f(theta - ei - ej)) / (4 * h[i] * h[j])
se = np.sqrt(np.diag(np.linalg.inv(H)))
print("observed information H:\n", H)
print("SE from inverse Hessian:", se)

boot = np.array([fit(rng.choice(x, x.size))[0] for _ in range(300)])
print("SE from 300 bootstrap refits:", boot.std(axis=0, ddof=1))
print("mean delivery k*scale:", round(theta.prod(), 3), "| sample mean:", round(x.mean(), 3))

MLE (k, scale): [2.7463 1.1727] | converged: True | NLL: 117.0368
scipy gamma.fit: [2.7463 1.1727]
observed information H:
 [[ 26.2958  51.1627]
 [ 51.1627 119.8146]]
SE from inverse Hessian: [0.4741 0.2221]
SE from 300 bootstrap refits: [0.4577 0.1717]
mean delivery k*scale: 3.221 | sample mean: 3.221

৭ · MAP — prior যোগ করা MLE

Bayes-এর নিয়মে posterior ∝ likelihood × prior. MAP (maximum a posteriori) হলো posterior-এর সর্বোচ্চ বিন্দু। Log নিলে:

θ̂_MAP = argmax_θ [ ℓ(θ) + log p(θ) ] = argmin_θ [ NLL(θ) − log p(θ) ]

অর্থাৎ MAP হলো MLE যেখানে loss-এর সঙ্গে একটি শাস্তি −log p(θ) যোগ হয়েছে। দুটি গুরুত্বপূর্ণ উদাহরণ:

Prior কতটা টানে? Beta-র সূত্রটি একটু সাজালে θ̂_MAP = w·(k/n) + (1 − w)·m, যেখানে m = (α − 1)/(α + β − 2) prior-এর mode আর w = n/(n + α + β − 2) ডেটার ওজন। নিচের stepper Beta(২, ১৮) prior নিয়ে (গড় ০.১, অর্থাৎ “আমাদের দোকানে সাধারণত দশে একটি ফেরত”) ডেটা বাড়ার সঙ্গে টানটি দেখায়।

৮ · MLE কখন ভুল পথে যায় — ছোট n আর পুরো আলাদা ডেটা

MLE-এর সুন্দর গুণগুলো “বড় n”-এর প্রতিশ্রুতি। ছোট n-এ বা ডেটার বিশেষ গঠনে এটি চরম উত্তর দেয়, অর্থাৎ overfitting-এর সবচেয়ে সরল রূপ। দুটি পরিচিত ক্ষেত্র:

• ছোট n, বিরল ঘটনা: ১০টি অর্ডারে শূন্য ফেরত হলে θ̂ = ০. পরে একটিও ফেরত এলে তার log loss −log ০ = ∞. Naive Bayes-এ না দেখা শব্দের একই সমস্যা pseudo-count দিয়ে সারানো হয়, যা Beta বা Dirichlet prior-এর MAP.
• পুরো আলাদা করা যায় এমন ডেটা: দুই শ্রেণিকে একটি রেখা পুরোপুরি ভাগ করলে (linearly separable), logistic regression-এর ওজনের দিক ঠিক রেখে মাপ বাড়ালে প্রতিটি সারির সম্ভাবনা ১-এর আরও কাছে যায়, আর likelihood কেবল বাড়তেই থাকে। সসীম MLE নেই (perfect separation); Albert আর Anderson ১৯৮৪ সালে এই অস্তিত্বের শর্তগুলো লিখেছিলেন।

নিচের ব্লক দুটিই সংখ্যায় দেখায়: (ক) ২০,০০০টি কাল্পনিক নতুন পণ্য, প্রতিটির মাত্র ১০টি অর্ডার, আসল হার ০.০৮; MLE আর Beta(২, ১৮)-এর MAP-এর ভুল তুলনা। (খ) ২০টি পুরো আলাদা বিন্দুতে gradient ascent, শাস্তি ছাড়া (MLE) আর Gaussian prior-সহ (MAP, λ = ১), ১০ থেকে ১০,০০০ ধাপে ওজনের মাপ।

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(5)             # (ক) নতুন পণ্য: আসল ফেরতের হার 0.08, মাত্র ১০টি অর্ডার
a, b, n, true = 2, 18, 10, 0.08            # prior Beta(2, 18): গড় 0.10
k = rng.binomial(n, true, size=20000)      # ২০,০০০টি কাল্পনিক পণ্য
mle, map_ = k / n, (k + a - 1) / (n + a + b - 2)
print(f"share with k = 0 (MLE says 'never returned'): {np.mean(k == 0):.3f}")
for name, est in (("MLE", mle), ("MAP", map_)):
    print(f"{name}: bias {est.mean() - true:+.4f}  RMSE {np.sqrt(np.mean((est - true) ** 2)):.4f}")

x = np.r_[rng.normal(-2, 0.6, 10), rng.normal(2, 0.6, 10)]   # (খ) দুই দল, পুরো আলাদা
X, y = np.c_[np.ones(20), x], np.r_[np.zeros(10), np.ones(10)]
print("separable:", x[:10].max() < x[10:].min())
ts = (10, 100, 1000, 10000)

def run(lam):                              # gradient ascent: (log-lik − (λ/2)·w1²)/n
    w, keep = np.zeros(2), []
    for t in range(1, ts[-1] + 1):
        p = 1 / (1 + np.exp(-X @ w))
        w += 0.5 * (X.T @ (y - p) - lam * w * np.r_[0, 1]) / 20
        if t in ts: keep.append(w[1])
    return keep

res = {lam: run(lam) for lam in (0.0, 1.0)}
for lam, ws in res.items():
    print(f"lambda = {lam}: " + "  ".join(f"t={t}: w1={v:.2f}" for t, v in zip(ts, ws)))
plt.figure(figsize=(6, 3.4), layout="constrained")
for lam, ws in res.items():
    plt.semilogx(ts, ws, "o-", lw=2, label=f"lambda = {lam}")
plt.xlabel("gradient steps (log scale)"); plt.ylabel("slope weight w1")
plt.title("separable data: MLE weight keeps growing, MAP stops"); plt.legend(); plt.show()

share with k = 0 (MLE says 'never returned'): 0.435
MLE: bias +0.0002  RMSE 0.0862
MAP: bias -0.0156  RMSE 0.0345
separable: True
lambda = 0.0: t=10: w1=1.46  t=100: w1=2.66  t=1000: w1=4.08  t=10000: w1=5.65
lambda = 1.0: t=10: w1=1.31  t=100: w1=1.55  t=1000: w1=1.55  t=10000: w1=1.55

ডেস্কটপে সত্যিকারের run: একই ধরনের ২০টি পুরো আলাদা বিন্দু (seed একই, তবে প্রথমে binomial টানা হয়নি বলে সংখ্যাগুলো আলাদা) statsmodels আর scikit-learn-এ দিলে কী হয়? লেখার সময়ে Python 3.14.4, statsmodels 0.15.0 আর scikit-learn 1.9.1 দিয়ে, একটি uv venv-এ চালানো:

import numpy as np
import statsmodels.api as sm
from sklearn.linear_model import LogisticRegression

rng = np.random.default_rng(5)
x = np.r_[rng.normal(-2, 0.6, 10), rng.normal(2, 0.6, 10)]
y = np.r_[np.zeros(10), np.ones(10)]

res = sm.Logit(y, sm.add_constant(x)).fit(disp=0)
print("statsmodels Logit params:", res.params.round(2), "bse:", res.bse.round(1))

for C in (1e10, 1.0):
    m = LogisticRegression(C=C, max_iter=10000).fit(x.reshape(-1, 1), y)
    print(f"sklearn C={C:g}: coef {m.coef_[0, 0]:.2f}, intercept {m.intercept_[0]:.2f}")

~/bazarghar-ml$ python separation.py
…/statsmodels/base/optimizer.py:536: PerfectSeparationWarning: Perfect separation or prediction detected, parameter may not be identified
  callback(newparams)
… (একই সতর্কবার্তা মোট ১৫ বার)
…/statsmodels/discrete/discrete_model.py:268: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  mlefit = super().fit(
statsmodels Logit params: [ 0.66 23.34] bse: [73615.1 74358.1]
sklearn C=1e+10: coef 7.43, intercept 0.23
sklearn C=1: coef 1.68, intercept 0.17

statsmodels ১৫ বার PerfectSeparationWarning দিয়েছে, তারপর জানিয়েছে optimization মেলেনি; ঢাল ২৩.৩৪, আর তার standard error প্রায় ৭৪,০০০. এটিই অংশ ৫-এর সূত্রের চরম রূপ: likelihood-এর চূড়া নেই, ঢালের দিকে বক্রতা প্রায় শূন্য, তাই 1/√J বিশাল। scikit-learn-এ C=1e10 (প্রায় শাস্তিহীন) ৭.৪৩ দিয়েছে কোনো সতর্কবার্তা ছাড়াই: solver-এর থামার সহনসীমা যেখানে পৌঁছেছে, সংখ্যাটি সেখানকার, “আসল” কোনো মান নয়। Default C=1 (Gaussian prior, τ² = ১) দিয়েছে ১.৬৮, একটি সসীম আর পুনরুৎপাদনযোগ্য MAP. শিক্ষা: শাস্তিহীন fit-এর বিশাল coefficient বা SE দেখলে প্রথমে separation সন্দেহ করুন।

আরও কয়েকটি জায়গায় MLE বিপদে পড়ে: Gaussian mixture-এ একটি উপাদান একটি বিন্দুর ওপর বসে σ → ০ করলে likelihood অসীম হয় (পর্ব ৩৯ (Gaussian mixture আর EM algorithm — নরম clustering)-এ সমাধান); parameter-এর সংখ্যা n-এর সঙ্গে বাড়লে MLE inconsistent হতে পারে (Neyman আর Scott-এর ১৯৪৮ সালের উদাহরণ, যেখানে variance-এর MLE কখনো আসল মানে পৌঁছায় না); আর model ভুল হলে SE-র সূত্রও ভুল হয়।

৯ · AI-এর যুগে — likelihood-এর কোড আর derivation যাচাই

AI assistant MLE-এর কোড দ্রুত লেখে, আর derivation-ও লিখে দেয় আত্মবিশ্বাসের সঙ্গে। ভুলগুলো সাধারণত চুপচাপ: চিহ্নের ভুল, ধ্রুবক বাদ পড়া, ভুল parameterization (scipy-র scale বনাম rate), বা optimizer না মিললেও ফল ছাপা। নিচের টেবিল এই বিষয়ের সাধারণ ভুলগুলো ধরে।

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

• Distribution-এর parameterization (rate না scale, কোনটি স্থির) নিজের scipy version-এর docs থেকে পড়েছি
• ছোট একটি ডেটায় closed form বা grid-এর সঙ্গে optimizer-এর উত্তর মিলিয়েছি
• res.success দেখেছি, আর শুরুর বিন্দু বদলেও একই উত্তর পেয়েছি
• SE সর্বোচ্চ বিন্দুর Hessian থেকে নিয়েছি, আর bootstrap দিয়ে যাচাই করেছি
• Model তুলনা held-out log-likelihood-এ করেছি, test set শুধু শেষে একবার
• শাস্তিহীন fit-এর বিশাল coefficient বা SE দেখে separation বা ছোট n পরীক্ষা করেছি
• AI-এর derivation একটি ছোট সংখ্যার উদাহরণে নিজে হাতে মিলিয়েছি, আর কোনো আসল গ্রাহকের ডেটা পেস্ট করিনি