Skip to content

সমাধান — অধ্যায় ৮.২ · Simulation Study

অধ্যায় ফাইল: part-8-capstone/08-02-simulation-study.md (§৭ অনুশীলনী)। সব সিমুলেশন seed np.random.default_rng(20260619)-এ চালানো; সংখ্যাগত উত্তর reproducible। প্রয়োজনীয় import: import numpy as np, from scipy import stats। মনে রাখবেন সংখ্যা-টানার ক্রম ও size অভিন্ন না-হলে canonical মান হুবহু নাও মিলতে পারে (draw-order dependence, ← §৪.১)।

মূল সংজ্ঞা ও নীতি। Monte-Carlo simulation — একটি পরিচিত তাত্ত্বিক ফল পুনরুৎপাদনের জন্য নকশা-করা দৈব পরীক্ষা, যা একই সঙ্গে তত্ত্ব ও কোড উভয় যাচাই করে (← §৪.৬)। Monte-Carlo standard error (MC SE): গড়-অনুমানে \(\approx s/\sqrt R\), অনুপাত-অনুমানে \(\approx\sqrt{\hat p(1-\hat p)/R}\); অর্ধেক করতে \(R\) চারগুণ। KS-দূরত্ব সর্বত্র \(D_n=\sup_x\lvert F_n(x)-\Phi(x)\rvert\)common random numbers (CRN): তুলনীয় শর্তে একই দৈব ইনপুট ব্যবহার করে পার্থক্যের ভেদ কমানো (variance reduction)। canonical মান। E1 — CLT log-log ঢাল \(-0.4998\) (Berry–Esseen পূর্বানুমান \(-0.5\)), \(C=0.1326\), KS \(n=5\to0.0590\), \(n=40\to0.0225\), \(n=160\to0.0108\), \(n=320\to0.0071\); E2 — bootstrap coverage \(0.9130\) (MC SE \(0.0063\)), গড় CI-প্রস্থ \(0.5869\), নামমাত্র \(0.95\); E3 — optimal degree \(5\), ন্যূনতম total test error \(0.5784\), irreducible \(\sigma^2=0.4900\); E4 — consistency RMSE log-log ঢাল \(-0.5570\) (পূর্বানুমান \(-0.5\)), asymptotic sd \(=\lambda=2.0000\), standardized-KS \(0.0180\) (\(n=200\), \(R=40000\)), bias ধনাত্মক কিন্তু \(\to0\)


ক · ধারণাগত (conceptual)

সমাধান ১ (★)

(ক) কেন পুনরুৎপাদন একই সঙ্গে দুটি জিনিস প্রমাণ করে। একটি পরিচিত তাত্ত্বিক ফল পুনরুৎপাদন করলে একই কাজে দুটো ভিন্ন দাবি একসঙ্গে যাচাই হয়। (i) তত্ত্বটি সংখ্যায় সত্য: যেমন E1-এ পর্যবেক্ষিত log-log ঢাল \(-0.4998\) Berry–Esseen-ঘোষিত তাত্ত্বিক \(-0.5\)-এর সঙ্গে মেলা নিশ্চিত করে CLT-অভিসরণ প্রকৃতই \(1/\sqrt n\) হারে ঘটে। (ii) সিমুলেশন-pipeline ঠিকঠাক: যেহেতু তত্ত্বটি আগে থেকেই প্রমাণিত (← Part VII), মিল পাওয়া মানে কোড, নমুনায়ন ও পরিসংখ্যান-গণনার পুরো শৃঙ্খল নির্ভুলভাবে কাজ করছে — নয়তো জানা উত্তরটা বেরোত না (← §৪.৬)।

(খ) ঢাল \(-0.25\) বেরোলে সন্দেহ কোথায়। সন্দেহ প্রথমে কোডে, তত্ত্বে (Berry–Esseen) নয় — কারণ Berry–Esseen সুপ্রতিষ্ঠিত ও প্রমাণিত, তাই অসঙ্গতির সবচেয়ে সম্ভাব্য উৎস কোডের একটি বাগ (যেমন মানক-করণে \(\sqrt n\) ভুল জায়গায়, বা KS-গণনায় ভুল রেফারেন্স-বণ্টন, বা draw-order-জনিত সংখ্যা-দূষণ)।

(গ) কেন পরিচিত ফল বাছা হয়। একটি জানা উত্তরই একমাত্র যা তত্ত্ব ও কোড উভয়কে একসঙ্গে যাচাই করার সুযোগ দেয় — এই অন্তর্নির্মিত যাচাই পেরিয়ে তবেই সেই সরঞ্জাম কোনো নতুন, অজানা প্রশ্নে বিশ্বাসযোগ্যভাবে প্রয়োগ করা যায়।

এক বাক্যে: একটা জানা ফল পুনরুৎপাদন একই সঙ্গে তত্ত্বকে সংখ্যায় নিশ্চিত করে ও pipeline যাচাই করে — মিল না-এলে সন্দেহটা প্রথমে কোডের বাগে, প্রমাণিত তত্ত্বে নয়।


সমাধান ২ (★★)

(ক) Berry–Esseen উপপাদ্য। iid \(X_1,\dots,X_n\), mean (গড়) \(\mu\), variance (ভেদ) \(\sigma^2\in(0,\infty)\), সসীম তৃতীয় absolute moment \(\rho=\mathbb E\lvert X-\mu\rvert^3<\infty\) হলে, মানক যোগফলের empirical CDF \(F_n\)\(N(0,1)\)-CDF \(\Phi\)-এর মধ্যে Kolmogorov দূরত্ব একটি সর্বজনীন ধ্রুবক \(C\) দিয়ে আবদ্ধ: $$ \sup_x\lvert F_n(x)-\Phi(x)\rvert\ \le\ \frac{C\,\rho}{\sigma^3\sqrt n}. $$ অর্থাৎ CLT-এর convergence (অভিসরণ) কেবল ঘটে না, তার হার-ও \(1/\sqrt n\)-এ বাঁধা।

(খ) log-log ঢাল \(-\tfrac12\) সীমার ডান পাশে একমাত্র \(n\)-নির্ভর অংশ \(1/\sqrt n\); বাকি সব (\(C,\rho,\sigma\)) \(n\)-নিরপেক্ষ ধ্রুবক। তাই বাস্তবে KS একটি ধ্রুবক \(C'=C\rho/\sigma^3\) গুণ \(n^{-1/2}\) হারে চলে: $$ \mathrm{KS}\ \approx\ C'\,n^{-1/2}\ \Longrightarrow\ \log(\mathrm{KS})\ \approx\ \log C'-\tfrac12\log n. $$ এটি \(\log(\mathrm{KS})\)-বনাম-\(\log n\)-এ একটি সরলরেখা, যার ঢাল ঠিক \(-\tfrac12\) — যার সঙ্গে canonical \(-0.4998\) মেলে।

(গ) skewness কেবল ধ্রুবক বদলায়, ঢাল নয়। Exp(1)-এর skewness \(2\) ধনাত্মক ও বড়, তাই \(\rho/\sigma^3\) বড়, ফলে \(C'\) (অর্থাৎ intercept, canonical \(C=0.1326\)) বাড়ে; কিন্তু \(n\)-এর ঘাত \(-\tfrac12\) অপরিবর্তিত থাকে — হার বণ্টন-নিরপেক্ষ, কেবল রেখার উচ্চতা বণ্টন-নির্ভর।

এক বাক্যে: Berry–Esseen-এর একমাত্র \(n\)-নির্ভরতা \(1/\sqrt n\), তাই \(\log(\mathrm{KS})\)-বনাম-\(\log n\) ঢাল \(-\tfrac12\) (canonical \(-0.4998\)); base-বণ্টনের skewness কেবল intercept \(C\) বদলায়, ঢাল নয়।


সমাধান ৩ (★★)

(ক) কেন right-skewed ডেটায় percentile bootstrap under-cover করে। base বণ্টন Exp(1) দৃঢ়ভাবে right-skewed, তাই \(n=40\)-এর মতো মাঝারি নমুনায় নমুনা-গড়ের বণ্টন এখনো পুরোপুরি প্রতিসম নয় — এতে সামান্য ডান-লেজ থাকে। percentile bootstrap সরাসরি resample-বণ্টনের \(2.5\)\(97.5\) শতাংশক কেটে CI বানায়, কিন্তু সেই পদ্ধতি এই অবশিষ্ট অপ্রতিসমতা ঠিকভাবে ধরতে পারে না — দুই প্রান্তকে সমান-ভাবে ঢাকে, ফলে প্রকৃত সত্য-মান কিছুটা বেশি বার বাইরে পড়ে ⇒ empirical coverage \(0.9130\), নামমাত্র \(0.95\)-এর সামান্য নিচে (← 4.9)।

(খ) \(n\) বাড়ালে। \(n\to\infty\)-এ CLT-এর ফলে নমুনা-গড় ক্রমে normal (তাই প্রতিসম) হয়ে যায়, অপ্রতিসমতা মুছে যায়, তাই coverage \(0.95\)-এর দিকে উঠে আসার কথা।

(গ) কেন under-coverage প্রকৃত প্রভাব। পার্থক্য \(0.95-0.9130=0.0370\) হলো MC SE \(0.0063\)-এর প্রায় \(0.0370/0.0063\approx5.9\) গুণ — এত বড় বিচ্যুতি নিছক সিমুলেশন-কোলাহলে (দৈব ওঠানামায়) ব্যাখ্যাযোগ্য নয়, তাই under-coverage-টি বাস্তব প্রভাব (← §৪.৩)।

এক বাক্যে: Exp(1)-এর skewness \(n=40\)-এ নমুনা-গড়ের বণ্টনকে অপ্রতিসম রাখে বলে percentile bootstrap সামান্য under-cover করে (\(0.9130\)); \(n\) বাড়লে coverage \(0.95\)-এ উঠবে, আর \(0.0370\approx5.9\,\)MC SE বলে এই ঘাটতি প্রকৃত, কোলাহল নয়।


সমাধান ৪ (★★)

(ক) bias–variance বিয়োজন। একটি স্থির test-বিন্দু \(x\)-এ প্রত্যাশিত prediction error তিন অংশে ভাঙে: $$ \mathbb E\bigl[(y-\hat f(x))^2\bigr]\ =\ \underbrace{\bigl(\mathbb E[\hat f(x)]-f(x)\bigr)^2}{\operatorname{Bias}^2}\ +\ \underbrace{\operatorname{Var}!\bigl(\hat f(x)\bigr)}, $$ যেখানে }}\ +\ \underbrace{\sigma^2}_{\text{irreducible}\(\sigma^2\) হলো noise-জনিত irreducible ত্রুটি — যেকোনো মডেল দিয়েও কমানো অসম্ভব (এখানে noise sd \(0.7\), তাই \(\sigma^2=0.4900\))।

(খ) কেন U-আকৃতি। degree বাড়লে মডেল নমনীয়তর হয়, তাই signal-এর বাঁক ভালো ধরে ⇒ \(\operatorname{Bias}^2\) কমে; কিন্তু একই নমনীয়তা noise-ও ধরে ফেলে (over-fit) ⇒ \(\operatorname{Var}\) বাড়ে। সরল মডেল signal ধরতে পারে না (under-fit, বড় bias), জটিল মডেল noise ধরে ফেলে (over-fit, বড় variance); দুইয়ের যোগফল তাই একটি optimal degree-এ ন্যূনতম — canonical সারণিতে degree \(5\)-এ total test error \(0.5784\) (← 6.1)। Variance monotone বাড়ে, Bias² monotone কমে/সমতল হয়, তাই Total U-আকৃতির।

(গ) কেন \(\operatorname{Bias}^2\) হঠাৎ \(2\to3\)-এ পড়ে। signal \(f(x)=\sin(1.5x)+0.5x\)-এর \(\sin\)-অংশের Taylor-প্রসারণে কেবল বিজোড় (odd) ঘাত থাকে (\(x,x^3,x^5,\dots\)); তাই একটি জোড়-সর্বোচ্চ-ঘাতের বহুপদী নতুন কিছু ধরে না, কিন্তু বিজোড় ঘাত যোগ হওয়ামাত্র (degree \(2\to3\)) প্রথম বড় bias-হ্রাস আসে (\(0.4812\to0.0741\)), এবং ফের \(4\to5\)-এ (\(0.0742\to0.0020\))।

এক বাক্যে: \(\mathbb E[(y-\hat f)^2]=\operatorname{Bias}^2+\operatorname{Var}+\sigma^2\); degree বাড়লে Bias² কমে কিন্তু Var বাড়ে, তাই Total U-আকৃতির (min degree \(5\), \(0.5784\)), আর \(\sin\)-এর বিজোড়-ঘাত-প্রকৃতির জন্য Bias² \(2\to3\)-এ হঠাৎ পড়ে।


খ · গণনামূলক (computational)

সমাধান ৫ (★)

(ক) MC SE গণনা। coverage একটি অনুপাত, তাই প্রতিটি ডেটাসেট একটি Bernoulli ফল ("CI সত্য মান ঢাকল কি না"); অনুপাত-অনুমানের MC SE সূত্র প্রযোজ্য। \(\hat p=0.9130\), \(D=2000\): $$ \operatorname{MC\,SE}(\hat p)=\sqrt{\frac{\hat p(1-\hat p)}{D}}=\sqrt{\frac{0.9130\times0.0870}{2000}}=\sqrt{\frac{0.079431}{2000}}=\sqrt{3.9716\times10^{-5}}=0.0063, $$ যা canonical \(0.0063\)-এর সঙ্গে হুবহু মেলে।

(খ) ৯৫% Monte-Carlo interval। আনুমানিক interval \(\hat p\pm1.96\,\operatorname{MC\,SE}\): $$ 0.9130\pm1.96\times0.0063=0.9130\pm0.0123=[0.9007,\ 0.9253]. $$ নামমাত্র \(0.95\) এই interval-এর বাইরে (\(0.9253<0.95\)) — অর্থাৎ পর্যবেক্ষিত under-coverage সিমুলেশন-কোলাহলে ব্যাখ্যাযোগ্য নয়, তাৎপর্যপূর্ণ।

(গ) MC SE অর্ধেক করতে। যেহেতু \(\operatorname{MC\,SE}\propto1/\sqrt D\), MC SE অর্ধেক (\(0.00315\)) করতে \(D\)-কে চারগুণ (\(8000\)) করতে হবে, কারণ \(\sqrt4=2\) — নির্ভুলতা ব্যয়বহুল (← §৪.২)।

এক বাক্যে: \(\sqrt{0.9130\times0.0870/2000}=0.0063\) (canonical); \(95\%\) MC interval \([0.9007,0.9253]\) নামমাত্র \(0.95\)-কে বাদ দেয় ⇒ তাৎপর্যপূর্ণ under-coverage; MC SE অর্ধেক করতে \(D\) চারগুণ।


সমাধান ৬ (★★)

(ক) asymptotic variance \(=\lambda^2\) MLE-এর সাধারণ asymptotic-normality ফল বলে, regularity-শর্ত পূরণে \(\sqrt n(\hat\lambda-\lambda)\xrightarrow{d}N\bigl(0,1/I(\lambda)\bigr)\), যেখানে \(I(\lambda)\) একক-পর্যবেক্ষণের Fisher information (← 4.3)। Exp(rate \(\lambda\))-এ \(I(\lambda)=1/\lambda^2\), তাই $$ \text{asymptotic variance}=\frac{1}{I(\lambda)}=\frac{1}{1/\lambda^2}=\lambda^2\ \Longrightarrow\ \text{asymptotic sd}=\sqrt{\lambda^2}=\lambda. $$

(খ) সংখ্যাগত মান। \(\lambda=2\) বসিয়ে asymptotic sd \(=\lambda=2.0000\), যা canonical \(2.0000\)-এর সঙ্গে মেলে।

(গ) কেন \(S\)-এর sd \(\approx1\) তত্ত্ব সমর্থন করে। \(S=\sqrt n(\hat\lambda-\lambda)/\lambda\)-কে \(\lambda\) দিয়ে ভাগ করায় তার তাত্ত্বিক asymptotic sd হয় \(\lambda/\lambda=1\); \(n=200\)-এ পর্যবেক্ষিত নমুনা-sd \(1.0091\) প্রায় \(1\), তাই \(N(0,1)\)-সীমা ঠিকঠাক (সঙ্গে KS \(0.0180\)-ও ছোট)।

এক বাক্যে: Exp-rate MLE-এর asymptotic variance \(=1/I(\lambda)=\lambda^2\), তাই asymptotic sd \(=\lambda=2.0000\); standardized \(S\)-এর sd \(1.0091\approx1\) সেই \(N(0,1)\)-সীমা নিশ্চিত করে।


সমাধান ৭ (★)

(ক) \(n=40\to160\) স্কেলিং। Berry–Esseen বলে \(\mathrm{KS}\approx C'\,n^{-1/2}\)\(n\) চারগুণ হলে (\(40\to160\)) \(n^{-1/2}\) ঠিক অর্ধেক হয় (\(\sqrt{40/160}=\sqrt{1/4}=1/2\)), তাই $$ \mathrm{KS}(160)\ \approx\ \tfrac12\times\mathrm{KS}(40)=\tfrac12\times0.0225=0.01125. $$

(খ) canonical-এর সঙ্গে তুলনা। পূর্বানুমান \(0.01125\) বনাম canonical প্রকৃত \(0.0108\) — খুব কাছাকাছি। হুবহু সমান না-হওয়ার কারণ, Berry–Esseen কেবল প্রধান পদ \(n^{-1/2}\) দেয়, কিন্তু সসীম \(n\)-এ উচ্চতর-ক্রম সংশোধন (\(o(n^{-1/2})\)-পদ) থাকে, যা মাঝারি \(n\)-এ ক্ষীণ কিন্তু শূন্য নয়; তাই সামান্য ব্যবধান।

(গ) কেন সাত-বিন্দু ধারা প্রায়-সরলরেখা। \(n=5\to0.0590\) থেকে \(n=320\to0.0071\) — এই পুরো ধারায় প্রধান পদ \(n^{-1/2}\)-ই আধিপত্য করে, তাই \(\log(\mathrm{KS})\)\(\log n\)-এর সম্পর্ক প্রায়-রৈখিক (ঢাল \(-0.4998\))।

এক বাক্যে: \(\mathrm{KS}\propto n^{-1/2}\)-এ \(n:40\to160\) চারগুণ হলে KS অর্ধেক (\(0.0225\to0.01125\)), যা canonical \(0.0108\)-এর কাছাকাছি — সামান্য ব্যবধান \(o(n^{-1/2})\)-সংশোধনের জন্য, আর প্রধান পদের আধিপত্যে সাত-বিন্দু log-log ধারা প্রায়-সরলরেখা।


গ · প্রমাণ/ডিজাইন-ভিত্তিক (proof / design-based)

সমাধান ৮ (★★)

দাবি ও নকশা। তত্ত্ব: iid \(X_i\sim N(\mu,\sigma^2)\) হলে \(\dfrac{(n-1)s^2}{\sigma^2}\sim\chi^2_{n-1}\), যেখানে \(s^2\) নমুনা-ভেদ (← 4.1)। এই দাবি যাচাইয়ের Monte-Carlo পরীক্ষা:

(ক) নকশা। - স্থির factor: \(\mu=0\), \(\sigma=1\), \(n=10\) (একটিমাত্র স্থির নমুনা-আকার; এখানে কোনো factor বদলানো হয় না, কেবল একটি বণ্টন-দাবি যাচাই)। - replication-সংখ্যা: \(R=50000\) বড় রাখা হয়, যাতে empirical বণ্টনের MC-কোলাহল (\(\sim1/\sqrt R\)) ছোট থাকে ও KS-অনুমান নির্ভরযোগ্য হয়। - প্রতি replication-এ গণনা: একটি নমুনা \(X_1,\dots,X_n\) টেনে \(s_r^2\) (ddof \(=1\)) বের করে \(T_r=(n-1)s_r^2/\sigma^2\)। - মিলিয়ে-দেখা: সংগৃহীত \(\{T_r\}\)-এর empirical বণ্টনকে \(\chi^2_{n-1}\)-এর সঙ্গে দুইভাবে মেলানো যায় — (i) KS-দূরত্ব stats.kstest(T, 'chi2', args=(n-1,)).statistic ছোট কি না; এবং (ii) moment-মিল: \(\mathbb E[\chi^2_{n-1}]=n-1\)\(\operatorname{Var}=2(n-1)\), তাই \(\bar T\approx n-1=9\) ও নমুনা-ভেদ \(\approx2(n-1)=18\) কি না।

import numpy as np
from scipy import stats
rng = np.random.default_rng(20260619)
mu, sigma, n, R = 0.0, 1.0, 10, 50_000
X  = rng.normal(mu, sigma, size=(R, n))            # iid N(0,1) নমুনা
s2 = X.var(axis=1, ddof=1)                          # নমুনা-ভেদ
T  = (n - 1) * s2 / sigma**2                        # তত্ত্বে ~ chi2_{n-1}
print("mean(T) =", round(T.mean(), 4), "(তত্ত্ব n-1 =", n-1, ")")
print("var(T)  =", round(T.var(ddof=0), 4), "(তত্ত্ব 2(n-1) =", 2*(n-1), ")")
print("KS(T, chi2_{n-1}) =", round(stats.kstest(T, 'chi2', args=(n-1,)).statistic, 4))

(খ) কেন বীজ স্থির জরুরি। বীজ default_rng(20260619) স্থির না-রাখলে ফল পুনরুৎপাদনযোগ্য নয় — অন্য কেউ (বা ভিন্ন সময়ে একই কোড) হুবহু একই সংখ্যা পাবে না, তাই canonical-মান মেলানো বা স্বাধীন যাচাই অসম্ভব হয়ে পড়ে (← §৪.১)।

(গ) কেন এটিও pipeline-যাচাইয়ের দৃষ্টান্ত। \(\chi^2_{n-1}\)-ফল একটি প্রমাণিত তত্ত্ব; empirical বণ্টন এর সঙ্গে মিললে একই কাজে তত্ত্বটি সংখ্যায় নিশ্চিত হয় এবং নমুনা-ভেদ-গণনা কোড (ddof, স্কেলিং) নির্ভুল বলে যাচাই হয় — ঠিক "জানা ফল পুনরুৎপাদন করে pipeline যাচাই"-এর নীতি (← §৪.৬)।

এক বাক্যে: \(\mu,\sigma,n\) স্থির রেখে বড় \(R\)-এ \(T_r=(n-1)s_r^2/\sigma^2\) সংগ্রহ করে \(\chi^2_{n-1}\)-এর সঙ্গে KS বা moment-মিলে যাচাই করলে তত্ত্ব ও নমুনা-ভেদ-কোড উভয়ই একসঙ্গে পরীক্ষিত হয় — বীজ স্থির থাকায় ফল পুনরুৎপাদনযোগ্য।


সমাধান ৯ (★★)

(ক) কেন এটি variance reduction (common random numbers)। ধরা যাক দুই degree-এর total error \(A\)\(B\)। আগ্রহ তাদের পার্থক্য \(A-B\)-এ (কোন degree ভালো)। ভেদের নিয়ম: $$ \operatorname{Var}(A-B)=\operatorname{Var}A+\operatorname{Var}B-2\operatorname{Cov}(A,B). $$ দুই degree-এ একই noise realisation ব্যবহার করলে \(A\)\(B\) একই দৈব ইনপুট ভাগ করে, তাই তারা ধনাত্মকভাবে সহ-সম্পর্কিত (\(\operatorname{Cov}(A,B)>0\)); ফলে \(\operatorname{Var}(A-B)\) কমে — অর্থাৎ degree-তুলনা তীক্ষ্ণতর হয়। এটিই common random numbers, একটি variance-reduction কৌশল (← §৪.৪)।

(খ) স্বাধীন noise হলে কী ক্ষতি। প্রতিটি degree-এ স্বাধীন noise ব্যবহার করলে \(\operatorname{Cov}(A,B)=0\), তাই পার্থক্যের ভেদ বাড়ে; U-আকৃতি বক্ররেখা অমসৃণ হয়ে যায়, আর দুই কাছাকাছি degree-এর প্রকৃত পার্থক্য অতিরিক্ত কোলাহলে ঢেকে যেতে পারে ⇒ optimal-degree সিদ্ধান্ত অস্থির।

(গ) কেন এটি factor-বিচ্ছিন্নকরণের সরাসরি প্রয়োগ। noise স্থির রাখায় দুই degree-এর total-error-এর পরিবর্তন কেবল degree-এর ফল, দৈব-ভিন্নতার নয় — অর্থাৎ একমাত্র পরিবর্তনশীল factor থাকে degree, যা এক-factor (one-factor-at-a-time) নকশার শর্ত পূরণ করে।

এক বাক্যে: একই noise ব্যবহারে \(\operatorname{Cov}(A,B)>0\) হওয়ায় \(\operatorname{Var}(A-B)\) কমে (common random numbers), তাই degree-তুলনা তীক্ষ্ণ ও optimal-degree স্থির; স্বাধীন noise হলে বক্ররেখা কোলাহলে দুলত — এটি degree-কে একমাত্র factor রাখার সরাসরি প্রয়োগ।


সমাধান ১০ (★★★)

দাবি। Exp-rate MLE \(\hat\lambda=1/\bar X\)-এর asymptotic ভেদ \(\lambda^2\), অর্থাৎ \(\sqrt n(\hat\lambda-\lambda)\xrightarrow{d}N(0,\lambda^2)\)

(ক) \(\bar X\)-এর CLT। \(X_i\sim\text{Exp}(\lambda)\)-এ \(\mathbb E[X]=1/\lambda\)\(\operatorname{Var}(X)=1/\lambda^2\)। iid গড়ের CLT সরাসরি দেয় (← 3.4): $$ \sqrt n\Bigl(\bar X-\tfrac1\lambda\Bigr)\ \xrightarrow{d}\ N!\Bigl(0,\ \tfrac1{\lambda^2}\Bigr). $$

(খ) delta method। MLE \(\hat\lambda=g(\bar X)\) যেখানে \(g(t)=1/t\), তাই \(g'(t)=-1/t^2\); মূল্যায়ন-বিন্দু \(t=\mathbb E[X]=1/\lambda\)-তে \(g'(1/\lambda)=-1/(1/\lambda)^2=-\lambda^2\)। delta method বলে একটি মসৃণ রূপান্তরের asymptotic ভেদ \(=(g')^2\times(\text{মূল ভেদ})\): $$ \sqrt n(\hat\lambda-\lambda)\ \xrightarrow{d}\ N!\Bigl(0,\ [g'(1/\lambda)]^2\cdot\tfrac1{\lambda^2}\Bigr)=N!\Bigl(0,\ \lambda^4\cdot\tfrac1{\lambda^2}\Bigr)=N(0,\lambda^2). $$

(গ) কেন \(\lambda^2=1/I(\lambda)\) ও ঢাল \(-\tfrac12\) Exp(rate \(\lambda\))-এর Fisher information \(I(\lambda)=1/\lambda^2\), তাই \(1/I(\lambda)=\lambda^2\) — delta-method-জাত ভেদ ঠিক Cramér–Rao নিম্নসীমার সমান, অর্থাৎ MLE asymptotically efficient। যেহেতু asymptotic sd \(=\lambda\), RMSE \(\approx\lambda/\sqrt n\propto n^{-1/2}\), তাই \(\log(\mathrm{RMSE})\)-বনাম-\(\log n\)-এ ঢাল \(-\tfrac12\) (consistency, canonical \(-0.5570\))। \(\blacksquare\)

এক বাক্যে: \(\bar X\)-এর CLT-তে \(g(t)=1/t\)-এর delta method বসালে \((g')^2\cdot\tfrac1{\lambda^2}=\lambda^4\cdot\tfrac1{\lambda^2}=\lambda^2\) পাওয়া যায়, যা \(1/I(\lambda)\)-এর সমান (efficient), আর RMSE \(\propto n^{-1/2}\) বলে log-log ঢাল \(-\tfrac12\)


ঘ · কোডিং (coding)

সব স্নিপেট seed np.random.default_rng(20260619)-এ চালানো; সংখ্যাগত উত্তর reproducible। import: import numpy as np, from scipy import stats। সংখ্যা-টানার ক্রম ও size অভিন্ন না-হলে canonical মান হুবহু নাও মিলতে পারে (draw-order dependence, ← §৪.১)।

সমাধান ১১ (★★)

লক্ষ্য। iid \(\text{Exp}(1)\) (mean \(1\), sd \(1\), skewness \(2\))-এর মানক নমুনা-গড় \(Z_n=\sqrt n(\bar X_n-1)\)-এর \(N(0,1)\)-এর সাপেক্ষে KS-দূরত্ব \(n\in\{5,10,20,40,80,160,320\}\)-এ (প্রতিটিতে \(R=60000\)) মেপে log-log ঢাল বের করা — Berry–Esseen-পূর্বানুমিত \(-0.5\) পুনরুৎপাদন।

import numpy as np
from scipy import stats
rng = np.random.default_rng(20260619)
ns = np.array([5, 10, 20, 40, 80, 160, 320])
R1 = 60_000
ks = np.empty(len(ns))
for j, n in enumerate(ns):
    x = rng.standard_exponential(size=(R1, n))       # iid Exp(1)
    Z = np.sqrt(n) * (x.mean(axis=1) - 1.0)          # মানক নমুনা-গড়
    ks[j] = stats.kstest(Z, "norm").statistic        # sup|F_R - Phi|
slope, intercept = np.polyfit(np.log(ns.astype(float)), np.log(ks), 1)
for n, k in zip(ns, ks):
    print(f"  n={n:>3} | KS={k:.4f}")
print("log-log slope =", round(slope, 4), " C = e^intercept =", round(np.exp(intercept), 4))

ফল (canonical)।

\(n\) KS দূরত্ব \(C/\sqrt n\) ফিট
\(5\) \(0.0590\) \(0.0593\)
\(10\) \(0.0417\) \(0.0419\)
\(20\) \(0.0285\) \(0.0297\)
\(40\) \(0.0225\) \(0.0210\)
\(80\) \(0.0148\) \(0.0148\)
\(160\) \(0.0108\) \(0.0105\)
\(320\) \(0.0071\) \(0.0074\)

log-log ঢাল \(=-0.4998\) (Berry–Esseen পূর্বানুমান \(-0.5\)), \(C=e^{\text{intercept}}=0.1326\)

(খ) মিল। KS তালিকা \(n=5\to0.0590\), \(n=40\to0.0225\), \(n=160\to0.0108\), \(n=320\to0.0071\) — canonical লক্ষ্যের সঙ্গে হুবহু; ঢাল \(-0.4998\) তাত্ত্বিক \(-0.5\)-এর অত্যন্ত কাছে।

(গ) কেন এটি তত্ত্ব ও pipeline দুটোই যাচাই করে। ঢাল \(-0.5\)-এর এত কাছে আসা একই সঙ্গে Berry–Esseen-এর \(1/\sqrt n\)-হার সংখ্যায় নিশ্চিত করে এবং মানক-করণ, KS-গণনা ও নমুনায়ন-কোডের নির্ভুলতা প্রমাণ করে — মিল না-এলে সন্দেহ যেত প্রথমে কোডে (← §৪.৬)। (\(R=60000\) বড় রাখা হয় যাতে KS-অনুমানের কোলাহল \(1/\sqrt R\) বৃহত্তম \(n\)-এর সংকেত \(1/\sqrt n\)-এর নিচে থাকে।)


সমাধান ১২ (★★)

লক্ষ্য। Exp(1) থেকে সত্য গড় \(\theta=1\); \(n=40\); \(D=2000\) ডেটাসেট; প্রতিটিতে \(B=1000\) resample-এর percentile bootstrap CI নমুনা-গড়ের ওপর — empirical coverage \(0.9130\) (MC SE \(0.0063\)) ও গড় CI-প্রস্থ \(0.5869\) পুনরুৎপাদন।

import numpy as np
rng = np.random.default_rng(20260619)
n, D, B, theta = 40, 2000, 1000, 1.0
covered = np.zeros(D, dtype=bool)
widths  = np.empty(D)
for d in range(D):
    data = rng.standard_exponential(size=n)          # একটি Exp(1) ডেটাসেট
    idx  = rng.integers(0, n, size=(B, n))           # B টি bootstrap resample
    boot = data[idx].mean(axis=1)                    # গড়ের bootstrap বণ্টন
    lo, hi = np.percentile(boot, [2.5, 97.5])        # ৯৫% percentile CI
    covered[d] = (lo <= theta <= hi)
    widths[d]  = hi - lo
coverage = covered.mean()
mc_se = np.sqrt(coverage * (1 - coverage) / D)
print("coverage   =", round(coverage, 4))            # 0.9130
print("MC SE      =", round(mc_se, 4))               # 0.0063
print("mean width =", round(widths.mean(), 4))       # 0.5869

ফল (canonical)।

রাশি সিমুলেশন তত্ত্ব/নামমাত্র
empirical coverage \(0.9130\) \(0.95\)
Monte-Carlo SE \(0.0063\) \(\sqrt{0.9130\times0.0870/2000}\)
গড় CI-প্রস্থ \(0.5869\)

(গ) কেন coverage \(0.95\)-এর নিচে ও তা তাৎপর্যপূর্ণ। Exp(1) right-skewed, তাই \(n=40\)-এ নমুনা-গড়ের বণ্টন এখনো সামান্য অপ্রতিসম, যা percentile bootstrap ঠিকভাবে ধরতে না-পারায় সামান্য under-cover (\(0.9130\))। পার্থক্য \(0.95-0.9130=0.0370\) MC SE \(0.0063\)-এর প্রায় \(5.9\) গুণ ⇒ এই under-coverage প্রকৃত প্রভাব, সিমুলেশন-কোলাহল নয় (← §৪.৩)।


সমাধান ১৩ (★★)

লক্ষ্য। signal \(f(x)=\sin(1.5x)+0.5x\) on \([-3,3]\), noise sd \(\sigma=0.7\) (\(\sigma^2=0.4900\)); polynomial regression degree \(d=1,\dots,12\); \(n_{\text{train}}=30\) (স্থির design); \(R=500\) প্রশিক্ষণ-সেট, প্রতিটি degree-এ একই noise (common random numbers)। প্রতি degree-এ \(\operatorname{Bias}^2\), \(\operatorname{Var}\), Total বের করে optimal degree \(5\) (Total \(0.5784\)) ও U-আকৃতি পুনরুৎপাদন।

import numpy as np
SEED = 20260619
sigma, n_train, R3 = 0.7, 30, 500
degrees = np.arange(1, 13)
x_train = np.linspace(-3, 3, n_train)                # স্থির design
x_test  = np.linspace(-3, 3, 200)                    # ঘন টেস্ট-গ্রিড
signal  = lambda x: np.sin(1.5 * x) + 0.5 * x
f_train, f_test = signal(x_train), signal(x_test)
bias2, variance, total = (np.empty(len(degrees)) for _ in range(3))
for k, d in enumerate(degrees):
    rng_d = np.random.default_rng(SEED)              # প্রতি degree-এ একই noise (CRN)
    Vd, Vt = np.vander(x_train, d + 1), np.vander(x_test, d + 1)
    preds = np.empty((R3, len(x_test)))
    for r in range(R3):
        y = f_train + sigma * rng_d.standard_normal(n_train)
        coef, *_ = np.linalg.lstsq(Vd, y, rcond=None)
        preds[r] = Vt @ coef
    mean_pred   = preds.mean(axis=0)
    bias2[k]    = np.mean((mean_pred - f_test) ** 2)
    variance[k] = np.mean(preds.var(axis=0))
    total[k]    = bias2[k] + variance[k] + sigma**2
best = int(degrees[np.argmin(total)])
for d, b, v, t in zip(degrees, bias2, variance, total):
    print(f"  {d:>2} | Bias^2={b:.4f} | Var={v:.4f} | Total={t:.4f}")
print("optimal degree =", best, " min total =", round(total.min(), 4))

ফল (canonical)।

deg Bias² Var Total
\(1\) \(0.4812\) \(0.0305\) \(1.0016\)
\(2\) \(0.4812\) \(0.0450\) \(1.0162\)
\(3\) \(0.0741\) \(0.0598\) \(0.6239\)
\(4\) \(0.0742\) \(0.0737\) \(0.6378\)
\(5\) \(0.0020\) \(0.0865\) \(0.5784\) (min)
\(6\) \(0.0019\) \(0.1001\) \(0.5920\)
\(7\) \(0.0002\) \(0.1131\) \(0.6033\)
\(8\) \(0.0004\) \(0.1288\) \(0.6192\)
\(9\) \(0.0004\) \(0.1455\) \(0.6359\)
\(10\) \(0.0004\) \(0.1641\) \(0.6545\)
\(11\) \(0.0004\) \(0.1826\) \(0.6731\)
\(12\) \(0.0006\) \(0.2041\) \(0.6947\)

irreducible \(\sigma^2=0.4900\); optimal degree \(=5\), ন্যূনতম total test error \(=0.5784\)। Bias² monotone কমে (বিশেষত \(2\to3\)\(4\to5\)-এ তীব্র পতন), Var monotone বাড়ে, তাই Total U-আকৃতির — degree \(5\)-এ ন্যূনতম।

(গ) কেন একই noise optimal-degree সিদ্ধান্ত তীক্ষ্ণ করে। সব degree-এ একই \(R=500\)টি noise realisation পুনর্ব্যবহার করায় (common random numbers) পাশাপাশি degree-এর Total-এর পার্থক্য কম ভেদ পায় (\(\operatorname{Cov}>0\)), তাই U-বক্ররেখা মসৃণ ও ন্যূনতম-বিন্দু (degree \(5\)) স্থিরভাবে চিহ্নিত হয় — স্বাধীন noise হলে কাছাকাছি degree-গুলোর তুলনা কোলাহলে দুলত (← §৪.৪)।


সমাধান ১৪ (★★)

লক্ষ্য। মডেল Exp(rate \(\lambda=2\)), MLE \(\hat\lambda=1/\bar X\), \(I(\lambda)=1/\lambda^2\) ⇒ asymptotic sd \(=\lambda=2\)। (a) Consistency: \(n\in\{10,40,160,640\}\), \(R=4000\)-এ bias ও RMSE, এবং RMSE log-log ঢাল (\(-0.5570\))। (b) Asymptotic normality: \(n=200\), \(R=40000\)-এ standardized \(S=\sqrt n(\hat\lambda-\lambda)/\lambda\)-এর \(N(0,1)\)-KS (\(0.0180\))।

import numpy as np
from scipy import stats
rng = np.random.default_rng(20260619)
lam = 2.0
# (a) consistency
ns = np.array([10, 40, 160, 640]); R4 = 4000
bias = np.empty(len(ns)); rmse = np.empty(len(ns))
for j, n in enumerate(ns):
    x = rng.exponential(scale=1.0 / lam, size=(R4, n))   # Exp(rate 2)
    lam_hat = 1.0 / x.mean(axis=1)                       # MLE = 1/Xbar
    bias[j] = lam_hat.mean() - lam
    rmse[j] = np.sqrt(np.mean((lam_hat - lam) ** 2))
sl, ic = np.polyfit(np.log(ns.astype(float)), np.log(rmse), 1)
for n, b, r in zip(ns, bias, rmse):
    print(f"  n={n:>3} | bias={b:+.4f} | RMSE={r:.4f}")
print("RMSE log-log slope =", round(sl, 4))              # -0.5570
# (b) asymptotic normality
x = rng.exponential(scale=1.0 / lam, size=(40_000, 200))
lam_hat = 1.0 / x.mean(axis=1)
S = np.sqrt(200) * (lam_hat - lam) / lam                 # -> N(0,1)
print("asymptotic sd = lambda =", round(lam, 4))         # 2.0000
print("KS(S, N(0,1)) =", round(stats.kstest(S, "norm").statistic, 4))  # 0.0180

ফল (canonical)।

\(n\) bias RMSE \(C/\sqrt n\) ফিট
\(10\) \(+0.2182\) \(0.7965\)
\(40\) \(+0.0484\) \(0.3404\)
\(160\) \(+0.0114\) \(0.1639\)
\(640\) \(+0.0028\) \(0.0775\)

RMSE log-log ঢাল \(=-0.5570\) (consistency পূর্বানুমান \(-0.5\)); asymptotic sd \(=\lambda=2.0000\); standardized-KS \(=0.0180\) (ছোট ⇒ Normal-সীমা ঠিক)। bias ধনাত্মক কিন্তু \(\to0\) (\(n=10\to+0.2182\), \(n=640\to+0.0028\))।

(গ) কেন bias ধনাত্মক অথচ consistent। \(g(t)=1/t\) একটি convex (উত্তল) ফাংশন, তাই Jensen-এর অসমতায় \(\mathbb E[1/\bar X]>1/\mathbb E[\bar X]=\lambda\) — সসীম \(n\)-এ \(\hat\lambda=1/\bar X\)-এর bias ধনাত্মক। কিন্তু \(\bar X\xrightarrow{p}1/\lambda\) হওয়ায় \(\hat\lambda\xrightarrow{p}\lambda\), তাই MLE consistent; bias \(n\to\infty\)-এ \(0\)-এ নামে (RMSE-ও \(\propto n^{-1/2}\), ঢাল \(-0.5570\approx-0.5\))।

এক বাক্যে: RMSE log-log ঢাল \(-0.5570\) (\(1/\sqrt n\)-হার) ও standardized-KS \(0.0180\) (asymptotic sd \(=\lambda=2.0000\)) একসঙ্গে MLE-এর consistency ও asymptotic normality নিশ্চিত করে; \(1/\bar X\)-এর ধনাত্মক bias Jensen-জাত, কিন্তু \(\bar X\to1/\lambda\)