সমাধান — অধ্যায় ৪.২ · Point Estimation: The Method of Moments¶

অধ্যায় ফাইল: part-4-inference/04-02-method-of-moments.md (§৭ অনুশীলনী)। সংখ্যাগত উত্তর numpy/scipy দিয়ে যাচাইযোগ্য (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$টা স্বাধীন সমীকরণ লাগে। কম মেলালে (যেমন $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$ স্বাধীন ও সমবণ্টিত হলে $$ \hat\mu_k'=\frac1n\sum_{i=1}^n X_i^k \xrightarrow{P} \mathbb{E}[X^k]=\mu_k' \quad (n\to\infty). $$ অর্থাৎ বড় নমুনায় 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 বা সম্ভাবনা $>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_{\text{MoM}}=\frac{1}{\bar x}=\frac{1}{1.2}\approx 0.833. $$

(গ) এই $\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$: $$ \hat\theta_{\text{MoM}}=2\bar x=2\times4.0=8.0. $$

(খ) সবচেয়ে বড় 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: $$ \mu_1'=\mathbb{E}[X]=\mu,\qquad \mu_2'=\mathbb{E}[X^2]=\mathrm{Var}(X)+(\mathbb{E}[X])^2=\sigma^2+\mu^2. $$

(খ) MoM সমীকরণ ও সমাধান: sample moment-এর সমান ধরি — $$ \hat\mu=\hat\mu_1'=\bar X,\qquad \hat\sigma^2+\hat\mu^2=\hat\mu_2'=\frac1n\sum_i X_i^2. $$ প্রথমটা থেকে $\hat\mu_{\text{MoM}}=\bar X$. দ্বিতীয়টায় $\hat\mu=\bar X$ বসিয়ে: $$ \hat\sigma^2_{\text{MoM}}=\hat\mu_2'-\bar X^2=\frac1n\sum_i X_i^2-\bar X^2. $$

(গ) পরিচিত পরিচয় ব্যবহার করে, variance-এর "shortcut" formula (সূত্র): $$ \frac1n\sum_{i=1}^n (X_i-\bar X)^2 =\frac1n\sum_i X_i^2-2\bar X\cdot\underbrace{\frac1n\sum_i X_i}{=\,\bar X}+\bar X^2 =\frac1n\sum_i X_i^2-2\bar X^2+\bar X^2 =\frac1n\sum_i X_i^2-\bar X^2. $$ ডান পাশটা ঠিক (খ)-এ পাওয়া $\hat\sigma^2_{\text{MoM}}$, তাই $$ \boxed{\hat\sigma^2 $$ }}=\frac1n\sum_{i=1}^n (X_i-\bar X)^2.লক্ষণীয়: হরে $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$। ভাগ করে সরাসরি বিচ্ছিন্ন করা যায়: $$ \frac{(\mathbb{E}[X])^2}{\mathrm{Var}(X)}=\frac{(\alpha/\beta)^2}{\alpha/\beta^2}=\alpha \quad\Rightarrow\quad \hat\alpha=\frac{\bar x^2}{s^2}, \qquad \frac{\mathbb{E}[X]}{\mathrm{Var}(X)}=\frac{\alpha/\beta}{\alpha/\beta^2}=\beta \quad\Rightarrow\quad \hat\beta=\frac{\bar x}{s^2}. $$ সংখ্যায়: $$ \hat\beta=\frac{3.75}{4.7}\approx 0.798,\qquad \hat\alpha=\bar x\,\hat\beta=3.75\times0.798\approx 2.99 \;\;\left(\text{বা }\hat\alpha=\frac{3.75^2}{4.7}\approx2.99\right). $$ (Figure 4-এর $(\hat\alpha,\hat\beta)\approx(2.93,\,0.75)$-এর সাথে সঙ্গতিপূর্ণ — সেখানে data সামান্য ভিন্ন।)

গ · প্রমাণভিত্তিক (proof-based)¶

সমাধান ৯ (★★)¶

দাবি: $\hat\theta_{\text{MoM}}=g^{-1}(\hat\mu_1')\xrightarrow{P}\theta$, যেখানে $\mu_1'(\theta)=g(\theta)$, $g$ সন্তত ও বিপরীতযোগ্য (তাই $g^{-1}$-ও সন্তত)।

প্রমাণ (দুই ধাপ): 1. LLN (3.3): যেহেতু $\hat\mu_1'=\frac1n\sum_i X_i$ হলো i.i.d. পদের গড় এবং $\mathbb{E}[X]=\mu_1'(\theta)=g(\theta)$ সসীম, দুর্বল-LLN অনুযায়ী $$ \hat\mu_1'\xrightarrow{P} g(\theta)\quad(n\to\infty). $$ 2. Continuous Mapping Theorem: $g^{-1}$ তার domain-এ (বিশেষত $g(\theta)$ বিন্দুতে) সন্তত। CMT বলে: যদি $Y_n\xrightarrow{P}c$ এবং $h$ $c$-তে সন্তত, তবে $h(Y_n)\xrightarrow{P}h(c)$. এখানে $Y_n=\hat\mu_1'$, $c=g(\theta)$, $h=g^{-1}$: $$ \hat\theta_{\text{MoM}}=g^{-1}(\hat\mu_1')\xrightarrow{P} g^{-1}(g(\theta))=\theta. \qquad\blacksquare $$

একাধিক প্যারামিটার: vector $\hat{\boldsymbol\mu}'=(\hat\mu_1',\dots,\hat\mu_p')\xrightarrow{P}\boldsymbol\mu'(\boldsymbol\theta)$ (componentwise LLN), আর moment→প্যারামিটার মানচিত্র $\mathbf{g}^{-1}$ সন্তত হলে 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$) কঠোরভাবে উত্তল (কারণ $\phi''(x)=2/x^3>0$)। যেহেতু $\bar X>0$ প্রায় নিশ্চিতভাবে এবং অ-অবক্ষয়ী (non-degenerate, variance $>0$ সসীম $n$-এ), Jensen-এর অসমতা (কঠোর রূপ) দেয়: $$ \mathbb{E}[\phi(\bar X)]>\phi(\mathbb{E}[\bar X]) \;\Longrightarrow\; \mathbb{E}!\left[\frac{1}{\bar X}\right]>\frac{1}{\mathbb{E}[\bar X]}=\frac{1}{1/\lambda}=\lambda. $$ তাই $\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$, তাই $$ \mathbb{E}[2\bar X]=2\,\mathbb{E}[\bar X]=2\cdot\frac1n\sum_i\mathbb{E}[X_i]=2\cdot\frac{\theta}{2}=\theta. \quad\Rightarrow\text{unbiased.} $$

(খ) variance: $\mathrm{Var}(X)=\dfrac{\theta^2}{12}$ (Uniform-এর প্রমিত ফল)। স্বাধীনতায় $$ \mathrm{Var}(\bar X)=\frac{1}{n}\,\mathrm{Var}(X)=\frac{\theta^2}{12n}, \qquad \mathrm{Var}(2\bar X)=4\,\mathrm{Var}(\bar X)=\frac{4\theta^2}{12n}=\boxed{\frac{\theta^2}{3n}}. $$

(গ) 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$-এ)। অনুপাত: $$ \frac{\mathrm{Var}(\hat\theta_{\text{MoM}})}{\mathrm{Var}(\hat\theta_{\text{MLE}})} \approx\frac{\theta^2/(3n)}{\theta^2/n^2}=\frac{n^2}{3n}=\frac{n}{3}. $$ অর্থাৎ 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: $$ \mathrm{Var}\big(g^{-1}(\hat\mu_1')\big)\approx\big[h'(\mu_1')\big]^2\,\mathrm{Var}(\hat\mu_1') =\frac{1}{[g'(\theta)]^2}\cdot\frac{\mathrm{Var}(X)}{n} =\boxed{\frac{\mathrm{Var}(X)}{n\,[g'(\theta)]^2}}. $$

E2-তে যাচাই: $\mu_1'=g(\lambda)=1/\lambda\Rightarrow g'(\lambda)=-1/\lambda^2$; Exponential-এ $\mathrm{Var}(X)=1/\lambda^2$। বসিয়ে: $$ \mathrm{Var}(\hat\lambda)\approx\frac{1/\lambda^2}{n\,(1/\lambda^2)^2}=\frac{1/\lambda^2}{n/\lambda^4}=\frac{\lambda^4}{n\lambda^2}=\frac{\lambda^2}{n}. $$ তাই $\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}")

প্রত্যাশিত প্রবণতা (আনুমানিক):

$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 হিসেবেও ব্যবহৃত হ