সমাধান — অধ্যায় ৪.৩ · Maximum Likelihood Estimation (MLE)¶
অধ্যায় ফাইল:
part-4-inference/04-03-maximum-likelihood-estimation.md(§৭ অনুশীলনী)। সংখ্যাগত উত্তরnumpy/scipyদিয়ে যাচাইযোগ্য (seed উল্লেখ থাকলে reproducible)। মূল বস্তু — likelihood \(L(\theta)=\prod_i f(X_i;\theta)\), log-likelihood \(\ell(\theta)=\sum_i\log f(X_i;\theta)\), score \(\ell'(\theta)\), এবং MLE-নিয়ম "\(\ell(\theta)\) সর্বোচ্চ করো" (মসৃণ অভ্যন্তরীণ চূড়ায় \(\ell'(\theta)=0\) সমাধান)। চলমান উদাহরণ: E1 Bernoulli\((p)\) → \(\hat p=\bar X\); E2 Normal\((\mu,\sigma^2)\) → \((\bar X,\ \frac1n\sum(X_i-\bar X)^2)\); E3 Exponential\((\lambda)\) → \(\hat\lambda=1/\bar X\); E4 Uniform\((0,\theta)\) → \(\hat\theta=\max_i X_i\)।
ক · ধারণাগত (conceptual)¶
সমাধান ১ (★)¶
MLE-র ধাপগুলো: 1. likelihood লেখো: ধরে নেওয়া মডেল থেকে \(L(\theta)=\prod_{i=1}^n f(X_i;\theta)\) — observed data-কে \(\theta\)-এর ফাংশন হিসেবে দেখার সম্ভাব্যতা। 2. log নাও: \(\ell(\theta)=\log L(\theta)=\sum_{i=1}^n\log f(X_i;\theta)\)। 3. derivative নিয়ে শূন্য বসাও (score equation): \(\ell'(\theta)=0\) (একাধিক প্যারামিটারে প্রতিটার partial = 0)। 4. সমাধান করো ও যাচাই করো: এই সমীকরণ থেকে \(\hat\theta_{\text{MLE}}\) পাও; \(\ell''(\theta)<0\) (বা Hessian negative-definite) দেখে নিশ্চিত করো এটা সর্বোচ্চ, সর্বনিম্ন নয়।
কেন \(\ell=\log L\), \(L\) নয় — দুটো কারণ: - অন্তরকলন সহজ: \(\log\) গুণফলকে যোগফলে বদলায় — \(\log\prod_i f_i=\sum_i\log f_i\)। যোগফলের derivative পদ-পদ নেওয়া যায়, কিন্তু \(n\)টা পদের গুণফলের derivative (product rule বারবার) দুঃস্বপ্ন। - সংখ্যাগত স্থিতিশীলতা: অনেক ছোট সম্ভাবনার (প্রতিটা \(<1\)) গুণফল দ্রুত শূন্যের দিকে যায় (floating-point underflow), কিন্তু তাদের \(\log\)-এর যোগফল একটা স্বাভাবিক মাঝারি সংখ্যা থাকে।
চূড়া একই কেন: \(\log\) একটা কঠোরভাবে বর্ধমান (strictly increasing/monotonic) ফাংশন। তাই \(L(\theta_1)>L(\theta_2)\iff\ell(\theta_1)>\ell(\theta_2)\) — মানের ক্রম অপরিবর্তিত, কাজেই argmax একই (\(\arg\max_\theta L=\arg\max_\theta\ell\))। Figure 1-এ এটাই দৃশ্যমান: সবুজ রেখা (MLE) উপরের (\(L\)) ও নিচের (\(\ell\)) প্যানেলে হুবহু একই \(\lambda\)-তে।
সমাধান ২ (★)¶
বাক্যটির অর্থ: একই সূত্র \(f(x;\theta)\) দুই ভূমিকায় ব্যবহৃত হয়, কোন চলকটা স্থির তার ওপর নির্ভর করে: - density হিসেবে: \(\theta\) স্থির, \(x\) চলক — তখন \(f(x;\theta)\) data-র probability density (বা mass)। এর বৈশিষ্ট্য: \(x\)-জুড়ে যোগ/integral = 1 (\(\int f(x;\theta)\,dx=1\) বা \(\sum_x f=1\))। - likelihood হিসেবে: \(x\) (observed data) স্থির, \(\theta\) চলক — তখন একই সংখ্যাকে আমরা \(\theta\)-এর ফাংশন \(L(\theta)=f(x;\theta)\) হিসেবে পড়ি।
কেন \(\theta\)-জুড়ে integral ১ নয়: normalization (\(\int=1\)) প্রযোজ্য কেবল \(x\)-চলক-এর density-র জন্য — কারণ ওই অক্ষে probability "ভাগ হয়"। কিন্তু \(L(\theta)\) একটা probability নয়, এটা শুধু "এই data কতটা মানানসই বিভিন্ন \(\theta\)-তে" তার একটা পরিমাপ; \(\theta\)-অক্ষ বরাবর তার কোনো normalization-শর্ত নেই, integral যেকোনো কিছু হতে পারে। তাই MLE শুধু \(L(\theta)\)-এর সর্বোচ্চ বিন্দু (argmax) খোঁজে — area-র সাথে কোনো সম্পর্ক নেই (Figure 1-এ আমরা ইচ্ছেমতো \(L\)-কে scale করেছি peak = 1, কারণ scaling argmax বদলায় না)।
সমাধান ৩ (★★)¶
কেন derivative খাটে না: Uniform\((0,\theta)\)-এ likelihood $$ L(\theta)=\theta^{-n}\quad(\theta\ge\max_i x_i),\qquad L(\theta)=0\quad(\theta<\max_i x_i). $$ এর আকৃতি (Figure 4): \(\max_i x_i\)-এর বাঁয়ে শূন্য, ঠিক \(\max_i x_i\)-তে হঠাৎ লাফিয়ে সর্বোচ্চে, তারপর \(\theta\) বাড়লে \(\theta^{-n}\) ধরে একটানা হ্রাসমান। কোথাও ঢাল শূন্য হয় না — \(\theta>\max\) অঞ্চলে \(\ell'(\theta)=-n/\theta<0\) (কখনো \(0\) নয়); আর \(\max\) বিন্দুতে ফাংশনটা অন্তরকলনযোগ্যই নয় (একটা লাফ/কোণ আছে)। সর্বোচ্চ likelihood তাই একটা প্রান্ত-বিন্দুতে (boundary), মসৃণ অভ্যন্তরে নয়।
সাধারণ সাবধান-বার্তা: "\(\ell'(\theta)=0\) সমাধান করো" নিয়মটা MLE দেয় কেবল তখনই যখন (ক) সর্বোচ্চ parameter-space-এর অভ্যন্তরে পড়ে, এবং (খ) likelihood সেখানে মসৃণ (অন্তরকলনযোগ্য)। support নিজেই \(\theta\)-র ওপর নির্ভর করলে (যেমন Uniform-এ উপরের সীমা \(\theta\)) এই শর্ত ভাঙে, আর তখন MLE বের করতে হয় likelihood-এর আকৃতি সরাসরি বিচার করে — যান্ত্রিকভাবে derivative নিয়ে নয়। (এটা পরীক্ষায় একটা প্রিয় ফাঁদ।)
সমাধান ৪ (★★)¶
Invariance property: যদি \(\hat\theta_{\text{MLE}}\) হয় \(\theta\)-র MLE এবং \(g\) একটা ফাংশন, তবে \(g(\theta)\)-র MLE হলো $$ \widehat{g(\theta)}{\text{MLE}}=g\big(\hat\theta\big). $$ অর্থাৎ নতুন প্যারামিটারের MLE পেতে শুধু পুরোনো MLE-টা }\(g\)-এর ভেতরে বসিয়ে দাও — নতুন করে likelihood লিখে আবার maximize করার দরকার নেই।
প্রয়োগ (E3): Exponential-এ rate-এর MLE \(\hat\lambda_{\text{MLE}}=1/\bar X\)। mean lifetime \(\tau=g(\lambda)=1/\lambda\)। invariance দিয়ে সরাসরি: $$ \hat\tau_{\text{MLE}}=g(\hat\lambda_{\text{MLE}})=\frac{1}{\hat\lambda_{\text{MLE}}}=\frac{1}{1/\bar X}=\bar X. $$ অর্থপূর্ণ: গড় আয়ুর MLE হলো নিছক নমুনা-গড় আয়ু \(\bar X\) — কোনো নতুন গণনা ছাড়াই।
খ · গণনামূলক (computational)¶
সমাধান ৫ (★)¶
Bernoulli, \(n=10\), \(\sum_i x_i=7\)।
(ক) \(\ell(p)=\Big(\sum_i x_i\Big)\log p+\Big(n-\sum_i x_i\Big)\log(1-p)=7\log p+3\log(1-p)\).
(খ) \(\displaystyle \ell'(p)=\frac{7}{p}-\frac{3}{1-p}=0\ \Rightarrow\ 7(1-p)=3p\ \Rightarrow\ 7=10p\ \Rightarrow\ \hat p=0.7\).
(গ) হ্যাঁ — \(\hat p=\dfrac{\sum_i x_i}{n}=\dfrac{7}{10}=\bar x=0.7\). (Bernoulli-তে MLE সর্বদা নমুনা-অনুপাত।)
সমাধান ৬ (★)¶
Exponential, নমুনা \(\{0.5,1.5,1.0,2.0,0.5\}\), \(n=5\)।
(ক) \(\bar x=\dfrac{0.5+1.5+1.0+2.0+0.5}{5}=\dfrac{5.5}{5}=1.1\).
(খ) \(\hat\lambda_{\text{MLE}}=\dfrac{1}{\bar x}=\dfrac{1}{1.1}\approx0.909\).
(গ) invariance: \(\tau=1/\lambda\Rightarrow\hat\tau_{\text{MLE}}=1/\hat\lambda_{\text{MLE}}=\bar x=1.1\).
সমাধান ৭ (★★)¶
Normal, নমুনা \(\{4,6,5,7,3\}\), \(n=5\)।
(ক) MLE formula (সূত্র): \(\hat\mu_{\text{MLE}}=\bar X\), \(\hat\sigma^2_{\text{MLE}}=\dfrac1n\sum_i(X_i-\bar X)^2\) (§৪-এ উৎপাদিত)।
(খ) \(\bar x=\dfrac{4+6+5+7+3}{5}=\dfrac{25}{5}=5\). বিচ্যুতি \((x_i-\bar x)\): \(\{-1,1,0,2,-2\}\); বর্গ: \(\{1,1,0,4,4\}\), যোগ \(=10\). তাই $$ \hat\mu=5,\qquad \hat\sigma^2_{\text{MLE}}=\frac{10}{5}=2. $$
(গ) হরে \(n\) (এখানে \(5\)), \(n-1\) নয় — তাই MLE-র \(\hat\sigma^2\) পরিচিত নিরপেক্ষ (unbiased) নমুনা-variance \(s^2=\frac{1}{n-1}\sum(X_i-\bar X)^2=10/4=2.5\)-এর চেয়ে ছোট। কাজেই \(\hat\sigma^2_{\text{MLE}}\) একটু নিচের দিকে biased: \(\mathbb{E}[\hat\sigma^2_{\text{MLE}}]=\frac{n-1}{n}\sigma^2<\sigma^2\) (বিস্তারিত 4.4-এ; bias \(\to0\) যখন \(n\to\infty\))।
সমাধান ৮ (★★)¶
Uniform\((0,\theta)\), নমুনা \(\{2.3,5.1,4.4,1.8,3.9\}\)।
(ক) \(\hat\theta_{\text{MLE}}=\max_i x_i=5.1\).
(খ) \(\bar x=\dfrac{2.3+5.1+4.4+1.8+3.9}{5}=\dfrac{17.5}{5}=3.5\); \(\hat\theta_{\text{MoM}}=2\bar x=7.0\).
(গ) এখানে MoM \((7.0)\) > MLE \((5.1)\)। MLE কখনো \(\max_i x_i\)-এর চেয়ে ছোট হতে পারে না — কারণ সংজ্ঞা অনুযায়ী এটা ঠিক \(\max_i x_i\), আর model-শর্ত \(\theta\ge\max_i x_i\) স্বয়ংক্রিয়ভাবে মেটে। বিপরীতে MoM-এর \(2\bar x\) অন্য নমুনায় \(\max_i x_i\)-এর চেয়ে ছোট হয়ে অবৈধ হতে পারত (4.2 §৭ Q4/Q6 দ্রষ্টব্য) — MLE-তে সেই সমস্যা নেই।
গ · প্রমাণভিত্তিক (proof-based)¶
সমাধান ৯ (★★)¶
Bernoulli, \(X_i\overset{iid}{\sim}\text{Bernoulli}(p)\), \(f(x;p)=p^x(1-p)^{1-x}\), \(x\in\{0,1\}\)।
(ক) Likelihood: $$ L(p)=\prod_{i=1}^n p^{X_i}(1-p)^{1-X_i}=p^{\sum_i X_i}(1-p)^{\,n-\sum_i X_i}. $$ log-likelihood (ধরি \(T=\sum_i X_i\)): $$ \ell(p)=T\log p+(n-T)\log(1-p). $$
(খ) score: $$ \ell'(p)=\frac{T}{p}-\frac{n-T}{1-p}=0 \ \Rightarrow\ T(1-p)=(n-T)p \ \Rightarrow\ T-Tp=np-Tp \ \Rightarrow\ T=np, $$ অর্থাৎ \(\boxed{\hat p_{\text{MLE}}=\dfrac{T}{n}=\dfrac1n\sum_i X_i=\bar X}\)।
(গ) দ্বিতীয় derivative: $$ \ell''(p)=-\frac{T}{p^2}-\frac{n-T}{(1-p)^2}<0\quad\text{(সব }p\in(0,1)\text{-এ)}, $$ কারণ \(T\ge0\) ও \(n-T\ge0\) এবং হরগুলো ধনাত্মক। \(\ell\) সর্বত্র কঠোরভাবে অবতল (concave), তাই critical point-টি অনন্য সর্বোচ্চ। (প্রান্ত-কেস: \(T=0\) হলে \(\hat p=0\), \(T=n\) হলে \(\hat p=1\) — সীমানায়।)
সমাধান ১০ (★★)¶
Exponential, \(X_i\overset{iid}{\sim}\text{Exp}(\lambda)\), \(f(x;\lambda)=\lambda e^{-\lambda x}\), \(x>0\)।
(ক) \(\displaystyle \ell(\lambda)=\sum_{i=1}^n\log\!\big(\lambda e^{-\lambda X_i}\big)=\sum_{i=1}^n(\log\lambda-\lambda X_i)=n\log\lambda-\lambda\sum_{i=1}^n X_i.\)
(খ) score: $$ \ell'(\lambda)=\frac{n}{\lambda}-\sum_{i=1}^n X_i=0 \ \Rightarrow\ \lambda=\frac{n}{\sum_i X_i}=\frac{1}{\bar X}, $$ অর্থাৎ \(\boxed{\hat\lambda_{\text{MLE}}=1/\bar X}\)।
(গ) \(\displaystyle \ell''(\lambda)=-\frac{n}{\lambda^2}<0\) — সর্বত্র অবতল, তাই critical point অনন্য সর্বোচ্চ (Figure 1-এর নিচের প্যানেলের concave আকৃতি ঠিক এটাই)।
সমাধান ১১ (★★★)¶
Uniform\((0,\theta)\), \(X_i\overset{iid}{\sim}\text{Uniform}(0,\theta)\), \(f(x;\theta)=\frac1\theta\) for \(0\le x\le\theta\), else \(0\)।
(ক) Likelihood: $$ L(\theta)=\prod_{i=1}^n f(X_i;\theta)=\prod_{i=1}^n\frac1\theta\,\mathbf{1}{0\le X_i\le\theta} =\theta^{-n}\prod_{i=1}^n\mathbf{1}{X_i\le\theta}. $$ এখন \(\prod_i\mathbf{1}\{X_i\le\theta\}=1\) কেবল তখনই যখন প্রতিটি \(X_i\le\theta\), অর্থাৎ যখন \(\max_i X_i\le\theta\); নয়তো অন্তত একটা indicator \(0\), পুরো গুণফল \(0\)। কাজেই $$ L(\theta)= \begin{cases} \theta^{-n}, & \theta\ge \max_i X_i,\[2pt] 0, & \theta<\max_i X_i. \end{cases} $$
(খ) \(\theta\ge\max_i X_i\) অঞ্চলে \(L(\theta)=\theta^{-n}\) একটা কঠোরভাবে হ্রাসমান ফাংশন (\(\theta\) বাড়লে \(\theta^{-n}\) কমে)। তাই এই অঞ্চলে \(L\) যত বেশি, \(\theta\) তত ছোট চাই; কিন্তু \(\theta\) অন্তত \(\max_i X_i\) হতেই হবে (নয়তো \(L=0\))। দুই শর্ত মিলিয়ে সর্বোচ্চ \(L\) ঠিক সীমানায়: $$ \boxed{\hat\theta_{\text{MLE}}=\max_i X_i}. $$ (Figure 4-এ এটাই — likelihood \(\max_i x_i\)-তে চূড়ায়, তার ডানে ক্ষয়ে যায়, বাঁয়ে শূন্য।)
(গ) \(\ell'(\theta)=0\) পদ্ধতি ব্যর্থ কারণ: (i) \(\theta>\max\) অঞ্চলে \(\ell(\theta)=-n\log\theta\), \(\ell'(\theta)=-n/\theta\ne0\) — কোথাও শূন্য নয়; (ii) চূড়াটা \(\theta=\max_i X_i\) বিন্দুতে, যেখানে \(L\) অন্তরকলনযোগ্য নয় (লাফ আছে); (iii) সর্বোচ্চ অভ্যন্তরে নয়, প্রান্তে (boundary)। derivative-শূন্য পদ্ধতির পূর্বশর্ত (মসৃণ, অভ্যন্তরীণ চূড়া) এখানে মেটে না — তাই আকৃতি সরাসরি বিচার করতে হলো।
সমাধান ১২ (★★★)¶
দাবি: \(g\) এক-এক (bijective) হলে \(\widehat{g(\theta)}_{\text{MLE}}=g(\hat\theta_{\text{MLE}})\)।
প্রমাণ: ধরি \(\psi=g(\theta)\), তাই \(\theta=g^{-1}(\psi)\) (\(g\) এক-এক বলে \(g^{-1}\) সুসংজ্ঞায়িত)। নতুন প্যারামিটারে likelihood: $$ \tilde L(\psi)=L\big(g^{-1}(\psi)\big). $$ আমরা চাই \(\tilde L(\psi)\)-এর argmax। যেহেতু \(g\) (ও \(g^{-1}\)) একটা এক-এক সংগতি (bijection), \(\psi\)-জুড়ে maximize করা আর \(\theta=g^{-1}(\psi)\)-জুড়ে maximize করা হুবহু একই কাজ — শুধু চলকের নাম বদলেছে, মানগুলোর সেট অপরিবর্তিত। তাই $$ \hat\psi=\arg\max_\psi \tilde L(\psi)=\arg\max_\psi L(g^{-1}(\psi)). $$ \(L(g^{-1}(\psi))\) সর্বোচ্চ হয় ঠিক যখন তার ভেতরের আর্গুমেন্ট \(L\)-কে সর্বোচ্চ করে, অর্থাৎ যখন \(g^{-1}(\psi)=\hat\theta_{\text{MLE}}\), অর্থাৎ \(\psi=g(\hat\theta_{\text{MLE}})\)। কাজেই $$ \hat\psi=g(\hat\theta_{\text{MLE}}),\qquad\text{অর্থাৎ}\quad \widehat{g(\theta)}{\text{MLE}}=g(\hat\theta).\qquad\blacksquare $$}
মন্তব্য: \(g\) এক-এক না হলেও ফলটা টেকে, তবে তখন "induced/profile likelihood" \(\tilde L(\psi)=\sup_{\theta:\,g(\theta)=\psi}L(\theta)\) দিয়ে সংজ্ঞায়িত করতে হয় (Wasserman §9.5; Casella–Berger Thm 7.2.10)। এই invariance-ই MLE-কে ব্যবহারিকভাবে এত সুবিধাজনক করে — যেকোনো রূপান্তরিত পরিমাণের MLE বিনা-পরিশ্রমে পাওয়া যায়।
ঘ · কোডিং (coding)¶
সমাধান ১৩ (★)¶
import numpy as np
rng = np.random.default_rng(0)
true_lam = 0.5
x = rng.exponential(scale=1/true_lam, size=200) # scale = 1/lambda
lam_hat = 1 / x.mean() # MLE of rate
tau_hat = x.mean() # MLE of mean lifetime = 1/lam_hat (invariance)
print(f"lambda_hat_MLE = {lam_hat:.4f} (true 0.5)")
print(f"tau_hat_MLE = {tau_hat:.4f} (true 1/0.5 = 2.0)")
নমুনা আউটপুট:
দুটোই সত্যিকারের মানের কাছাকাছি (\(n=200\) যথেষ্ট বড়)। লক্ষ করুন \(\hat\tau=\bar X=1/\hat\lambda\) — invariance-এর সরাসরি যাচাই।সমাধান ১৪ (★★)¶
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(7)
true_lam, n = 0.7, 12
x = rng.exponential(scale=1/true_lam, size=n)
sx, xbar = x.sum(), x.mean()
lam_hat = 1 / xbar # closed-form MLE
g = np.linspace(0.15, 2.2, 600)
loglik = n*np.log(g) - g*sx # l(lambda) = n log(lam) - lam * sum(x)
plt.plot(g, loglik, color="#1b6ca8", lw=2)
plt.axvline(lam_hat, color="#2e8b57", lw=2, label=f"1/x-bar = {lam_hat:.3f}")
plt.axvline(g[loglik.argmax()], color="#c0392b", ls="--",
label=f"grid argmax = {g[loglik.argmax()]:.3f}")
plt.xlabel("lambda"); plt.ylabel("log-likelihood"); plt.legend(); plt.show()
print("closed-form 1/x-bar :", round(lam_hat, 4))
print("grid argmax :", round(g[loglik.argmax()], 4))
যাচাই: grid-এর argmax আর বিশ্লেষণী \(1/\bar X\) কার্যত মেলে (পার্থক্য কেবল grid-এর সূক্ষ্মতার সীমায়, এখানে \(\approx0.003\)); grid যত ঘন (linspace-এ বেশি বিন্দু), তত নিখুঁত মেলে। বক্ররেখাটা concave, একটিমাত্র চূড়া — Figure 1-এর নিচের প্যানেলের পুনরুৎপাদন।
সমাধান ১৫ (★★★)¶
import numpy as np
from scipy.optimize import minimize_scalar
rng = np.random.default_rng(11)
true_theta, n, reps = 10.0, 20, 5000
S = rng.uniform(0, true_theta, size=(reps, n))
mle = S.max(axis=1) # MLE = max_i X_i
mom = 2 * S.mean(axis=1) # MoM = 2 * X-bar
mse = lambda est, t: np.mean((est - t)**2)
# (খ) summary
for name, e in [("MLE (max)", mle), ("MoM (2*xbar)", mom)]:
print(f"{name:14s} mean={e.mean():.3f} SD={e.std():.3f} "
f"bias={e.mean()-true_theta:+.3f} MSE={mse(e, true_theta):.3f}")
নমুনা আউটপুট:
MLE (max) mean=9.522 SD=0.452 bias=-0.478 MSE=0.434
MoM (2*xbar) mean=9.998 SD=1.304 bias=-0.002 MSE=1.700
(গ) কোনটার MSE ছোট, কেন: MLE-র MSE অনেক ছোট (\(\approx0.43\) বনাম \(\approx1.70\), প্রায় ৪ গুণ)। MoM প্রায় unbiased (bias \(\approx0\)) কিন্তু variance বিশাল (SD \(\approx1.30\)); MLE সামান্য নিচে-biased (\(\max_i X_i\le\theta\) সর্বদা, তাই গড়ে একটু কম — bias \(\approx-0.48\)) কিন্তু variance খুবই ছোট (SD \(\approx0.45\))। যেহেতু MSE = variance + bias², আর MLE-র variance-সুবিধা তার ছোট bias²-কে অনেক ছাড়িয়ে যায়, MLE সামগ্রিকভাবে জেতে — Figure 3-এর "MLE far tighter"-এর সংখ্যাগত নিশ্চিতকরণ। (এখানে MLE variance \(\sim\theta^2/n^2\), MoM variance \(\sim\theta^2/(3n)\) — বড় \(n\)-এ ব্যবধান আরও বাড়ে।)
(ঘ) বোনাস — সংখ্যাগত MLE = বিশ্লেষণী \(1/\bar X\) (Exponential):
x = rng.exponential(scale=1/0.5, size=300)
nll = lambda L: -(len(x)*np.log(L) - L*x.sum()) # negative log-likelihood
res = minimize_scalar(nll, bounds=(1e-3, 5), method="bounded")
print("numerical MLE :", round(res.x, 4))
print("1 / x-bar :", round(1/x.mean(), 4))
# numerical MLE : 0.5097
# 1 / x-bar : 0.5097