সমাধান — অধ্যায় ৪.২ · Point Estimation: The Method of Moments¶
অধ্যায় ফাইল:
part-4-inference/04-02-method-of-moments.md(§৭ অনুশীলনী)। সংখ্যাগত উত্তরnumpy/scipyদিয়ে verifiable (seed উল্লেখ থাকলে reproducible)। মূল বস্তু — population moment \(\mu_k'=\mathbb{E}[X^k]\), sample moment \(\hat\mu_k'=\frac1n\sum_i X_i^k\), MoM-নিয়ম "\(\mu_k'(\theta)=\hat\mu_k'\) সমাধান করো", এবং consistency/bias/variance বিশ্লেষণ। চলমান উদাহরণ: E1 Normal\((\mu,\sigma^2)\); E2 Exponential\((\lambda)\) → \(\hat\lambda=1/\bar X\); E3 Uniform\((0,\theta)\) → \(\hat\theta=2\bar X\); E4 Gamma\((\alpha,\beta)\)।
ক · ধারণাগত (conceptual)¶
সমাধান ১ (★)¶
MoM-এর ধাপগুলো: 1. population moment লেখো: মডেল থেকে তাত্ত্বিক moment বের করো প্যারামিটারের মাধ্যমে — \(\mu_k'(\theta)=\mathbb{E}_\theta[X^k]\)। (এক প্যারামিটারে সাধারণত শুধু \(\mu_1'\) লাগে; দুই প্যারামিটারে \(\mu_1'\) ও \(\mu_2'\)।) 2. sample moment হিসাব করো: data থেকে \(\hat\mu_k'=\frac1n\sum_i X_i^k\)। 3. সমীকরণ বসাও: তত্ত্ব = data, অর্থাৎ \(\mu_k'(\theta)=\hat\mu_k'\)। 4. সমাধান করো: এই সমীকরণ(গুলো) \(\theta\)-এর জন্য সমাধান করে \(\hat\theta_{\text{MoM}}\) পাও।
কয়টা moment? ঠিক যত প্যারামিটার তত moment — \(p\)টা অজানা প্যারামিটারের জন্য প্রথম \(p\)টা moment। কারণ \(p\)টা অজানা বের করতে \(p\)টা independent সমীকরণ লাগে। কম মেলালে (যেমন \(2\) প্যারামিটারে \(1\) সমীকরণ) সমাধান অনন্য নয় (under-determined); বেশি মেলালে (over-determined) সাধারণত একসাথে সব সমীকরণ মেলানো অসম্ভব। Figure 1-এ একটাই প্যারামিটার \(\lambda\), তাই একটাই সমীকরণ \(\bar x=1/\lambda\) — আর সেটাই যথেষ্ট ও অনন্য সমাধান \(\hat\lambda=1/\bar x\) দেয়। (E1 Normal-এ \(2\) প্যারামিটার, তাই \(2\) moment — সমাধান ৭ দ্রষ্টব্য।)
সমাধান ২ (★)¶
plug-in principle: তত্ত্ব আমাদের প্রায়ই বলে একটা পরিমাণ population-এর কোনো বৈশিষ্ট্যের (যেমন \(\mathbb{E}[X^k]\)) মাধ্যমে — কিন্তু সেই population-বৈশিষ্ট্য অজানা। plug-in idea: যেখানেই population-পরিমাণ দরকার, সেখানে তার sample-প্রতিরূপ "বসিয়ে দাও"। MoM ঠিক এর একটা উদাহরণ: population moment \(\mu_k'=\mathbb{E}[X^k]\)-এর জায়গায় sample moment \(\hat\mu_k'=\frac1n\sum_i X_i^k\) বসিয়ে moment-সমীকরণ সমাধান করি। (একই নীতি পরে empirical CDF, bootstrap ইত্যাদিতেও ফিরে আসবে।)
কেন যুক্তিসঙ্গত? কারণ Law of Large Numbers (3.3): \(X_1,\dots,X_n\) independent ও identically distributed হলে
অর্থাৎ বড় নমুনায় sample moment population moment-এর খুব কাছাকাছি — তাই "বসিয়ে দেওয়ার" ভুল ছোট, এবং \(n\to\infty\)-এ মিলিয়ে যায়। এটাই MoM-এর consistency-র মূল ভিত্তি (Figure 2-তে দৃশ্যমান)।
সমাধান ৩ (★★)¶
কেন প্রথম moment হুবহু মেলে, দ্বিতীয়টা নয়: MoM-এ আমরা estimator বেছেছি শুধু প্রথম moment-সমীকরণ সমাধান করে — \(\bar x=1/\hat\lambda\), অর্থাৎ \(\hat\lambda=1/\bar x\)। তাই গঠনগতভাবে (by construction) পাওয়া \(\hat\lambda\)-তে model-এর প্রথম moment \(\mathbb{E}[X]=1/\hat\lambda\) ঠিক \(\bar x\)-এর সমান হতেই বাধ্য — এটাই তো সমীকরণটা ছিল। কিন্তু দ্বিতীয় moment \(\mathbb{E}[X^2]=2/\hat\lambda^2\) আমরা estimator বানাতে ব্যবহার করিনি; তাই sample-এর দ্বিতীয় moment \(\hat\mu_2'\)-এর সাথে এর মেলার কোনো নিশ্চয়তা নেই — শুধু consistency-র কারণে কাছাকাছি আসে (Figure 1-এ ৫.২৫ বনাম ৫.৪২)।
সাধারণ সত্য: MoM যে moment-গুলো মেলায়, ঠিক সেগুলোই হুবহু মেলে; বাকি (না-মেলানো) moment-গুলো শুধু আনুমানিকভাবে মেলে। এটাই MoM-এর সীমাবদ্ধ "fit" — এটা পুরো distribution-কে নয়, শুধু নির্বাচিত কয়েকটা moment-কে নিখুঁত করে।
সমাধান ৪ (★★)¶
না, MoM সবসময় বৈধ প্যারামিটার দেয় না। MoM শুধু moment-সমীকরণ(গুলো) বীজগাণিতিকভাবে সমাধান করে — সেই সমাধান প্যারামিটার-স্থানের (parameter space) ভেতরে পড়বে, তার কোনো অন্তর্নিহিত নিশ্চয়তা নেই।
উদাহরণ ও কারণ: - E3 Uniform\((0,\theta)\): \(\hat\theta=2\bar X\) কোনো নমুনায় সবচেয়ে বড় observation \(X_{(n)}\)-এর চেয়ে ছোট হতে পারে। কিন্তু মডেল বলে সব \(X_i\le\theta\), তাই \(\theta\ge X_{(n)}\) হতেই হবে — \(2\bar X<X_{(n)}\) হলে সেই estimate মডেলের সাথে অসঙ্গত (সমাধান ৬খ দ্রষ্টব্য)। - আরও জটিল মডেল: কোনো কোনো mixture বা negative-binomial-এ moment-দিয়ে dispersion estimate করলে MoM ঋণাত্মক variance বা probability \(>1\) দিতে পারে, যা অর্থহীন।
MLE-র তুলনায় দুর্বলতা: MLE সংজ্ঞা অনুযায়ী likelihood-কে প্যারামিটার-স্থানের ওপর সর্বোচ্চ করে, তাই তার সমাধান সবসময় বৈধ সীমানার মধ্যেই থাকে (যেমন Uniform-এ \(\hat\theta_{\text{MLE}}=X_{(n)}\ge\) সব observation)। MoM-এর এই "বাইরে বেরিয়ে যাওয়া"র ঝুঁকিই তার একটা পরিচিত দুর্বলতা — যদিও বড় \(n\)-এ consistency-র কারণে এটা ক্রমে কমে।
খ · গণনামূলক (computational)¶
সমাধান ৫ (★)¶
নমুনা \(\{0.8,1.2,2.0,0.5,1.5\}\), \(n=5\)।
(ক) \(\bar x=\dfrac{0.8+1.2+2.0+0.5+1.5}{5}=\dfrac{6.0}{5}=1.2\).
(খ) E2-তে \(\mathbb{E}[X]=1/\lambda\), MoM সমীকরণ \(\bar x=1/\hat\lambda\):
(গ) এই \(\hat\lambda\)-তে \(\mathbb{E}[X]=1/\hat\lambda=1/0.833=1.2=\bar x\) — হুবহু মেলে। (এটাই matching: প্রথম moment গঠনগতভাবে সমান।)
সমাধান ৬ (★)¶
নমুনা \(\{2.1,5.4,3.3,7.8,1.4\}\)।
(ক) \(\bar x=\dfrac{2.1+5.4+3.3+7.8+1.4}{5}=\dfrac{20.0}{5}=4.0\). E3-তে \(\mathbb{E}[X]=\theta/2\), তাই MoM সমীকরণ \(\bar x=\hat\theta/2\):
(খ) সবচেয়ে বড় observation \(X_{(n)}=7.8\), আর \(\hat\theta=8.0>7.8\) — এবার ঠিক আছে (estimate সব observation-কে ঢেকে রেখেছে)। কিন্তু এটা সবসময় হয় না: যদি নমুনায় একটা মান \(\theta\)-র খুব কাছে থাকে কিন্তু বাকিরা ছোট, তখন \(\bar x\) ছোট হয়ে \(2\bar x<X_{(n)}\) হতে পারে — যা অসম্ভব, কারণ Uniform\((0,\theta)\)-এ সব মান \(\le\theta\)। এটাই MoM-এর একটা দুর্বলতা (সমাধান ৪, ১১ দ্রষ্টব্য); MLE \(X_{(n)}\) এই সমস্যা এড়ায়।
সমাধান ৭ (★★)¶
E1 Normal\((\mu,\sigma^2)\) — দুই প্যারামিটার, তাই দুটো moment।
(ক) population moments:
(খ) MoM সমীকরণ ও সমাধান: sample moment-এর সমান ধরি —
প্রথমটা থেকে \(\hat\mu_{\text{MoM}}=\bar X\). দ্বিতীয়টায় \(\hat\mu=\bar X\) বসিয়ে:
(গ) পরিচিত পরিচয় ব্যবহার করে, variance-এর "shortcut" formula (সূত্র):
ডান পাশটা ঠিক (খ)-এ পাওয়া \(\hat\sigma^2_{\text{MoM}}\), তাই
লক্ষণীয়: হরে \(n\) (না \(n-1\)) — তাই এটা সামান্য biased (\(\mathbb{E}[\hat\sigma^2_{\text{MoM}}]=\frac{n-1}{n}\sigma^2<\sigma^2\)); বড় \(n\)-এ পার্থক্য নগণ্য (asymptotically unbiased)।
সমাধান ৮ (★★)¶
E4 Gamma\((\alpha,\beta)\) (shape–rate), \(\mathbb{E}[X]=\alpha/\beta\), \(\mathrm{Var}(X)=\alpha/\beta^2\); মাপা \(\bar x=3.75\), \(s^2=4.7\)।
দুটো moment-সমীকরণ: \(\bar x=\hat\alpha/\hat\beta\) এবং \(s^2=\hat\alpha/\hat\beta^2\)। ভাগ করে সরাসরি বিচ্ছিন্ন করা যায়:
সংখ্যায়:
(Figure 4-এর \((\hat\alpha,\hat\beta)\approx(2.93,\,0.75)\)-এর সাথে consistent — সেখানে data সামান্য ভিন্ন।)
গ · প্রমাণভিত্তিক (proof-based)¶
সমাধান ৯ (★★)¶
দাবি: \(\hat\theta_{\text{MoM}}=g^{-1}(\hat\mu_1')\xrightarrow{P}\theta\), যেখানে \(\mu_1'(\theta)=g(\theta)\), \(g\) continuous ও বিপরীতযোগ্য (তাই \(g^{-1}\)-ও continuous)।
প্রমাণ (দুই ধাপ): 1. LLN (3.3): যেহেতু \(\hat\mu_1'=\frac1n\sum_i X_i\) হলো i.i.d. পদের গড় এবং \(\mathbb{E}[X]=\mu_1'(\theta)=g(\theta)\) সসীম, দুর্বল-LLN অনুযায়ী
- Continuous Mapping Theorem: \(g^{-1}\) তার domain-এ (বিশেষত \(g(\theta)\) বিন্দুতে) continuous। CMT বলে: যদি \(Y_n\xrightarrow{P}c\) এবং \(h\) \(c\)-তে continuous, তবে \(h(Y_n)\xrightarrow{P}h(c)\). এখানে \(Y_n=\hat\mu_1'\), \(c=g(\theta)\), \(h=g^{-1}\):
একাধিক প্যারামিটার: vector \(\hat{\boldsymbol\mu}'=(\hat\mu_1',\dots,\hat\mu_p')\xrightarrow{P}\boldsymbol\mu'(\boldsymbol\theta)\) (componentwise LLN), আর moment→প্যারামিটার মানচিত্র \(\mathbf{g}^{-1}\) continuous হলে vector-CMT দিয়ে \(\hat{\boldsymbol\theta}_{\text{MoM}}\xrightarrow{P}\boldsymbol\theta\)। (Figure 2 এই consistency-রই অভিজ্ঞতামূলক প্রদর্শন।)
সমাধান ১০ (★★)¶
দাবি: E2-তে \(\mathbb{E}[\hat\lambda]>\lambda\) সসীম \(n\)-এ, যেখানে \(\hat\lambda=1/\bar X\)।
প্রমাণ (Jensen): \(\bar X=\frac1n\sum X_i\), প্রতিটি \(X_i\sim\text{Exp}(\lambda)\), তাই \(\mathbb{E}[\bar X]=\mathbb{E}[X]=1/\lambda\)। ফাংশন \(\phi(x)=1/x\) ধনাত্মক অঞ্চলে (\(x>0\)) কঠোরভাবে convex (কারণ \(\phi''(x)=2/x^3>0\))। যেহেতু \(\bar X>0\) প্রায় নিশ্চিতভাবে এবং অ-অবক্ষয়ী (non-degenerate, variance \(>0\) সসীম \(n\)-এ), Jensen-এর inequality (কঠোর রূপ) দেয়:
তাই \(\mathbb{E}[\hat\lambda]=\mathbb{E}[1/\bar X]>\lambda\) — estimator উপরের দিকে biased। \(\blacksquare\)
(Figure 3-এ ৮০০০-নমুনার গড় \(0.510>0.5\) ঠিক এটাই দেখায়। যেহেতু \(\mathrm{Var}(\bar X)=1/(n\lambda^2)\to0\), \(\bar X\) ক্রমে \(1/\lambda\)-তে কেন্দ্রীভূত হয় ও bias \(\to0\) — অর্থাৎ asymptotically unbiased এবং consistent। প্রকৃতপক্ষে Exponential-এ \(\mathbb{E}[1/\bar X]=\frac{n}{n-1}\lambda\), যা স্পষ্টভাবে \(>\lambda\) ও \(n\to\infty\)-এ \(\lambda\)-তে নামে।)
সমাধান ১১ (★★★)¶
E3 Uniform\((0,\theta)\), \(\hat\theta_{\text{MoM}}=2\bar X\)।
(ক) unbiased: \(X\sim\text{Unif}(0,\theta)\Rightarrow\mathbb{E}[X]=\theta/2\), তাই
(খ) variance: \(\mathrm{Var}(X)=\dfrac{\theta^2}{12}\) (Uniform-এর প্রমিত ফল)। independence-এ
(গ) MLE-র সাথে তুলনা: \(\hat\theta_{\text{MLE}}=X_{(n)}=\max_i X_i\), যার (দেওয়া) \(\mathrm{Var}(X_{(n)})=\dfrac{n\theta^2}{(n+1)^2(n+2)}\approx\dfrac{\theta^2}{n^2}\) (বড় \(n\)-এ)। অনুপাত:
অর্থাৎ MoM-এর variance MLE-র প্রায় \(n/3\) গুণ বড় — বড় \(n\)-এ MLE নাটকীয়ভাবে বেশি efficient (এর error \(\sim1/n\) হারে কমে, MoM-এর \(\sim1/\sqrt n\))। সিদ্ধান্ত: এখানে MoM সরল ও unbiased হলেও MLE অনেক বেশি নিখুঁত — 4.3-এর মূল বার্তার (efficiency) একটা নাটকীয় পূর্বাভাস। (টীকা: \(X_{(n)}\) সামান্য biased — \(\mathbb{E}[X_{(n)}]=\frac{n}{n+1}\theta\) — কিন্তু MSE-তে তবু MoM-কে অনেক ব্যবধানে হারায়; bias-corrected \(\frac{n+1}{n}X_{(n)}\) unbiased ও আরও ভালো।)
সমাধান ১২ (★★)¶
delta method: \(\hat\theta_{\text{MoM}}=g^{-1}(\hat\mu_1')\), \(g\) মসৃণ, \(g'(\theta)\ne0\)। CLT দেয় \(\sqrt n\,(\hat\mu_1'-\mu_1')\xrightarrow{d}\mathcal{N}(0,\mathrm{Var}(X))\), অর্থাৎ \(\mathrm{Var}(\hat\mu_1')\approx\mathrm{Var}(X)/n\)। \(h=g^{-1}\) ধরে, \(h'(\mu_1')=\dfrac{1}{g'(\theta)}\) (inverse-function নিয়ম)। delta method:
E2-তে verify: \(\mu_1'=g(\lambda)=1/\lambda\Rightarrow g'(\lambda)=-1/\lambda^2\); Exponential-এ \(\mathrm{Var}(X)=1/\lambda^2\)। বসিয়ে:
তাই \(\text{SE}(\hat\lambda)\approx\lambda/\sqrt n\) — ঠিক Figure 3-এর box-এর সূত্র (\(\lambda=0.5,n=50\Rightarrow0.5/\sqrt{50}=0.0707\), যা empirical SD \(0.073\)-এর সাথে মেলে)। \(\blacksquare\)
ঘ · কোডিং (coding)¶
সমাধান ১৩ (★)¶
import numpy as np
rng = np.random.default_rng(0)
theta_true = 10.0
x = rng.uniform(0, theta_true, size=200)
theta_mom = 2 * x.mean() # MoM: theta_hat = 2 * x-bar
x_max = x.max() # the largest observation (relates to the MLE)
print(f"theta_mom = {theta_mom:.3f}")
print(f"x_max = {x_max:.3f} (true theta = {theta_true})")
নমুনা আউটপুট (seed 0): theta_mom = 9.98x, x_max = 9.9x — দুটোই \(10\)-এর কাছাকাছি। সাধারণত \(X_{(n)}=\max\) সত্যি \(\theta\)-র একটু নিচে থাকে (কারণ সর্বোচ্চ observation \(\theta\)-তে পৌঁছায় না), আর \(2\bar X\) দুদিকেই হেলতে পারে কিন্তু গড়ে unbiased। (একই বীজ চালালে \(X_{(n)}\) প্রায়ই \(\theta\)-র সামান্য কাছে ও স্থিতিশীল — কারণ এর variance \(\sim\theta^2/n^2\), \(2\bar X\)-এর variance \(\sim\theta^2/n\)-এর চেয়ে অনেক ছোট; সমাধান ১১গ।)
সমাধান ১৪ (★★)¶
import numpy as np
rng = np.random.default_rng(1)
lam_true = 0.5
ns = [5, 10, 50, 100, 500, 2000]
reps = 300
print(f"{'n':>6} {'mean(lam_hat)':>14} {'sd(lam_hat)':>12} {'lam/sqrt(n)':>12}")
for n in ns:
# vectorised: reps samples of size n; lam_hat = 1 / row-mean
ests = 1.0 / rng.exponential(scale=1.0/lam_true, size=(reps, n)).mean(axis=1)
print(f"{n:>6} {ests.mean():>14.4f} {ests.std():>12.4f} {lam_true/np.sqrt(n):>12.4f}")
expected প্রবণতা (আনুমানিক):
| \(n\) | mean\((\hat\lambda)\) | sd\((\hat\lambda)\) | \(\lambda/\sqrt n\) |
|---|---|---|---|
| 5 | ~0.63 | ~0.34 | 0.224 |
| 10 | ~0.56 | ~0.20 | 0.158 |
| 50 | ~0.51 | ~0.073 | 0.071 |
| 100 | ~0.505 | ~0.051 | 0.050 |
| 500 | ~0.501 | ~0.022 | 0.022 |
| 2000 | ~0.500 | ~0.011 | 0.011 |
লক্ষণীয়: (১) mean → \(0.5\) এবং ছোট \(n\)-এ সামান্য উপরে (Jensen-bias, সমাধান ১০); (২) sd → \(0\) এবং বড় \(n\)-এ \(\lambda/\sqrt n\)-কে ভালোভাবে অনুসরণ করে (delta-method SE, সমাধান ১২) — এটাই consistency-র সংখ্যাগত নিশ্চিতকরণ (Figure 2)।
সমাধান ১৫ (★★★)¶
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng(42)
alpha_true, beta_true = 3.0, 0.8 # shape, rate
x = rng.gamma(shape=alpha_true, scale=1/beta_true, size=400)
# --- (খ) MoM fit: match mean & variance ---
m1, s2 = x.mean(), x.var() # 1st moment & central 2nd moment
beta_mom = m1 / s2
alpha_mom = m1 * beta_mom # = m1**2 / s2
print(f"MoM: alpha_hat={alpha_mom:.3f}, beta_hat={beta_mom:.3f}")
# --- (ঘ) MLE fit for comparison (scipy: shape a, scale=1/rate; fix loc=0) ---
a_ml, loc_ml, scale_ml = stats.gamma.fit(x, floc=0)
beta_ml = 1.0 / scale_ml
print(f"MLE: alpha_hat={a_ml:.3f}, beta_hat={beta_ml:.3f}")
# --- (গ) overlay on the data histogram ---
grid = np.linspace(0, x.max()*1.02, 400)
plt.hist(x, bins=32, density=True, alpha=0.5, color="#d9822b", label="data (n=400)")
plt.plot(grid, stats.gamma.pdf(grid, a=alpha_mom, scale=1/beta_mom),
color="#2e8b57", lw=2.5, label=f"MoM fit (a={alpha_mom:.2f}, b={beta_mom:.2f})")
plt.plot(grid, stats.gamma.pdf(grid, a=a_ml, scale=scale_ml),
color="#7b3fa0", lw=2.0, ls="--", label=f"MLE fit (a={a_ml:.2f}, b={beta_ml:.2f})")
plt.xlabel("x"); plt.ylabel("density"); plt.legend(); plt.tight_layout(); plt.show()
নমুনা আউটপুট (seed 42): MoM: alpha_hat≈2.93, beta_hat≈0.75; MLE: alpha_hat≈3.0x, beta_hat≈0.8x। পর্যবেক্ষণ: দুটো fit-ই data-র সাথে চমৎকার মেলে এবং একে অপরের কাছাকাছি (Figure 4-এর সবুজ-নীল curve প্রায় মিলেছিল)। সূক্ষ্ম পার্থক্য: MLE সাধারণত data-র full shape-কে (শুধু দুটো moment নয়) ব্যবহার করে বলে একটু কম "noisy" এবং ছোট \(n\)-এ পার্থক্য স্পষ্ট হয় — MLE তখন সামান্য ভালো fit ও কম variance দেয় (সমাধান ১১, এবং 4.3-এর efficiency-আলোচনার পূর্বাভাস)। MoM-এর বড় সুবিধা: এটা closed-form (কোনো optimization লাগে না), তাই দ্রুত — প্রায়ই MLE-র optimization-এর starting point হিসেবেও ব্যবহৃত হয়।