সমাধান — অধ্যায় ৫.৫ · GLM: Poisson Regression & Beyond¶

অধ্যায় ফাইল: part-5-modeling/05-05-glm-poisson-regression.md (§৭ অনুশীলনী)-এর সব প্রশ্নের পূর্ণ সমাধান। সব সংখ্যাগত ফল NumPy/SciPy দিয়ে যাচাই ও statsmodels-এর সাথে মিলিয়ে দেখা হয়েছে; seeded কোড reproducible। চলমান উদাহরণ — bike-rental count, seed np.random.default_rng(20260619), $n=250$: predictor temp (uniform $5$–$35$°C) ও weekend (Bernoulli $2/7$); true $\log\mu=1.5+0.06\cdot\text{temp}+0.4\cdot\text{weekend}$, count $\sim$ Gamma($6,1/6$)–Poisson mixture (overdispersed)।

স্মারক — canonical সংখ্যা (seed 20260619): count mean $19.56$, var $205$ (var/mean $\approx10.5$); Poisson $\hat\beta=[1.560,\,0.0597,\,0.301]$ (intercept, temp, weekend), $\mathrm{SE}=[0.0476,\,0.00178,\,0.0307]$; rate ratio temp $e^{\hat\beta_1}=1.0616$ (per $+1$°C), per $+5$°C $e^{5\hat\beta_1}=1.348$, weekend $e^{\hat\beta_2}=1.351$; residual deviance $1093.1$, Pearson $\chi^2=1096.3$, $df=247$, dispersion $\hat\phi=4.44$, AIC $2237.9$; NB $\hat\beta=[1.574,\,0.0592,\,0.300]$, $\mathrm{SE}=[0.0874,\,0.00374,\,0.0689]$, $\alpha=0.179$, AIC $1753.3$; prediction temp$=25$, weekend$=1\Rightarrow\hat\mu\approx28.6$। মূল সম্পর্ক: $\mu_i=e^{x_i^\top\beta}$, Poisson-এ $\operatorname{Var}=\mu$, NB-তে $\operatorname{Var}=\mu+\alpha\mu^2$।

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

সমাধান ১ (★)¶

count outcome-এ OLS (identity link) ব্যবহারের তিন মৌলিক সমস্যা:

(ক) ঋণাত্মক পূর্বাভাস — অর্থহীন count। OLS-এ $\hat\mu=x^\top\hat\beta$ যেকোনো বাস্তব সংখ্যা নিতে পারে। কম temp-এ (যেমন $5$°C, weekday) linear fit সহজেই $\hat\mu<0$ দিতে পারে — কিন্তু "$-3$টা bike" বলে কিছু নেই। log link এই সমস্যা মূলেই কাটায়: $\mu=e^{x^\top\beta}>0$ সব ইনপুটে, তাই পূর্বাভাস সর্বদা বৈধ (ধনাত্মক) count-rate।
(খ) heteroscedasticity — 'E' (equal variance) ভাঙে। count data Poisson-সদৃশ, যেখানে $\operatorname{Var}(y\mid x)=\mu$ — অর্থাৎ যত বেশি প্রত্যাশিত count, তত বেশি spread। OLS ধরে নেয় error-variance সব $x$-এ স্থির ($\sigma^2$, LINE-এর 'E')। এখানে variance mean-এর সাথে বাড়ে, তাই OLS-এর SE ও inference ভুল হয়।
(গ) Normal নয় — 'N' ভাঙে। ছোট count-এ ($\mu\approx2$–$3$) outcome discrete ($0,1,2,\dots$) ও right-skewed (নিচে $0$-তে আটকানো, উপরে লম্বা লেজ)। OLS ধরে নেয় error Normal (অবিচ্ছিন্ন, প্রতিসম)। small-count regime-এ এই অনুমান স্পষ্টভাবে মিথ্যা, ফলে CI/p-value বিশ্বাসযোগ্য নয়।

(সংক্ষেপে: OLS-এর তিন ভিত্তি — ধনাত্মক support, equal variance, Normality — তিনটাই count outcome-এ ভাঙে; GLM ঠিক এই তিনটা সারাতে link + distribution বদলায়।)

সমাধান ২ (★)¶

(ক) positivity। exponential function $e^{(\cdot)}$ যেকোনো বাস্তব ইনপুট $\eta\in(-\infty,\infty)$-কে $(0,\infty)$-তে map করে। তাই $\mu=e^{x^\top\beta}$ সর্বদা কঠোরভাবে ধনাত্মক — coefficient বা predictor যাই হোক, $\mu\le0$ হওয়া অসম্ভব। এ-ই log link-এর মূল আকর্ষণ: positivity constraint স্বয়ংক্রিয়, আলাদা করে চাপাতে হয় না।

(খ) multiplicative effect। $\mu=e^{\beta_0+\beta_1 x_1+\beta_2 x_2}=e^{\beta_0}\cdot e^{\beta_1 x_1}\cdot e^{\beta_2 x_2}$। অর্থাৎ linear predictor-এ predictor-গুলো যোগ হলেও, mean count-এ তারা গুণ হয়ে কাজ করে — প্রতিটা predictor একটা গুণক $e^{\beta_j}$ অবদান রাখে। তাই "$x_1$ এক একক বাড়লে $\mu$ একটা নির্দিষ্ট গুণে ($e^{\beta_1}$ গুণ) বদলায়, একটা নির্দিষ্ট পরিমাণ যোগ হয় না" — এ-ই multiplicative interpretation।

(গ) canonical link ও সুবিধা। Poisson-এর canonical link হলো log ($g(\mu)=\log\mu$)। canonical link-এ score equation সরলতম রূপ নেয় — $X^\top(y-\mu)=0$ (সমাধান ১৩) — এবং observed Fisher information = expected Fisher information ($X^\top WX$, $W=\operatorname{diag}(\mu_i)$), ফলে Newton–Raphson ও Fisher scoring একই হয়ে IRLS সরল ও সংখ্যাগতভাবে স্থিতিশীল হয়।

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

(ক) equidispersion ধর্ম। Poisson-এর সংজ্ঞাগত ধর্ম: $\operatorname{Var}(y\mid x)=\mathbb E[y\mid x]=\mu$ — mean আর variance সমান। চলমান উদাহরণে temp$=25$, weekend$=1$-এ $\hat\mu=e^{1.560+0.0597\cdot25+0.301}=e^{3.3535}\approx28.6$। সুতরাং Poisson model দাবি করে সেখানে $\operatorname{Var}(y)\approx28.6$ (SD $\approx\sqrt{28.6}\approx5.35$) — অর্থাৎ ওই পরিস্থিতির count-গুলো মোটামুটি $28.6\pm5.3$-এ ঘোরাঘুরি করবে।

(খ) data-র সাথে সংঘর্ষ। কিন্তু পুরো data-তে count mean $19.56$, var $205$ — অনুপাত $205/19.56\approx10.5$। Poisson বলে এই অনুপাত $\approx1$ হওয়ার কথা; এখানে variance mean-এর $\sim10$ গুণ। মানে: data-তে যত spread আছে, Poisson model তার একটা ভগ্নাংশ মাত্র আশা করে — এ-ই overdispersion। (এটা শুধু overall দেখা; conditional dispersion §৮-এ $\hat\phi=4.44$ দিয়ে মাপা হয়েছে।)

(গ) কারণ — unobserved heterogeneity। এমন overdispersion সাধারণত আসে অদৃশ্য variability থেকে: প্রতিটি দিনের আসল rate শুধু temp আর weekend-এ ঠিক হয় না — আবহাওয়ার সূক্ষ্ম তারতম্য, স্থানীয় event, ছুটির ধরন ইত্যাদি rate-কে দিনে-দিনে ওঠানামা করায়। গাণিতিকভাবে এটা Gamma-mixed Poisson (rate নিজেই random) — যা ঠিক আমাদের data-generating process (rate = rng.gamma(6,1/6))। এই extra layer-ই variance বাড়িয়ে দেয়।

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

(ক) SE underestimate, ফলাফল। Poisson model variance $=\mu$ ধরে নেয়, কিন্তু overdispersion-এ আসল variance $\phi\mu$ ($\phi>1$)। ফলে Poisson-হিসাবের SE আসলের চেয়ে ছোট (underestimated, প্রায় $\sqrt\phi$ গুণ ছোট)। এর শৃঙ্খল-প্রতিক্রিয়া: $$z=\frac{\hat\beta}{\mathrm{SE}}\ \text{(স্ফীত)}\ \Rightarrow\ p\text{-value (কৃত্রিমভাবে ছোট)}\ \Rightarrow\ \text{CI (অতি-সংকীর্ণ)}.$$

(খ) বিপদ — over-confident / false positive। SE ছোট হওয়ায় inference anti-conservative: যে predictor আসলে significant নয়, তাকেও "significant" দেখানোর ঝুঁকি বাড়ে — অর্থাৎ বেশি false positive (Type I error)। বাস্তব গবেষণায় এটা ভুল আবিষ্কারের দিকে নিয়ে যায়।

(গ) দুই প্রতিকার। একই log-link mean model রেখে শুধু variance ঠিক করা: (i) quasi-Poisson — SE-কে $\sqrt{\hat\phi}$ দিয়ে স্ফীত করা ($\operatorname{Var}=\phi\mu$); (ii) negative binomial (NB) — variance-এ explicit extra-dispersion parameter $\alpha$ যোগ করা ($\operatorname{Var}=\mu+\alpha\mu^2$)।

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

(ক) variance functions। | model | variance | বৈশিষ্ট্য | |---|---|---| | Poisson | $\operatorname{Var}=\mu$ | equidispersion (single param) | | quasi-Poisson | $\operatorname{Var}=\phi\mu$ | mean-এর সাথে linear (রৈখিক) scaling ($\phi$ constant) | | negative binomial | $\operatorname{Var}=\mu+\alpha\mu^2$ | quadratic — variance mean-এর চেয়ে দ্রুত বাড়ে |

(খ) likelihood বনাম quasi-likelihood। NB একটা সত্যিকারের probability distribution (Gamma-mixed Poisson), তাই তার পূর্ণ likelihood আছে ⇒ log-likelihood, AIC, BIC সবই সংজ্ঞায়িত ও model-তুলনায় ব্যবহারযোগ্য (NB AIC $=1753.3$)। quasi-Poisson কোনো পূর্ণ distribution ধরে না — শুধু mean-variance সম্পর্ক ($\operatorname{Var}=\phi\mu$) নির্দিষ্ট করে estimating equation সমাধান করে (quasi-likelihood)। তাই এর কোনো প্রকৃত likelihood নেই ⇒ সাধারণ AIC প্রযোজ্য নয়; এটা দেয় কেবল সংশোধিত SE ও একই $\hat\beta$।

(গ) $\alpha\to0$। NB variance $\mu+\alpha\mu^2\xrightarrow{\alpha\to0}\mu$ — অর্থাৎ NB তখন Poisson-এ ফিরে যায়। তাই Poisson হলো NB-র একটা nested special case ($\alpha=0$); NB হলো Poisson-এর Gamma-mixed সাধারণীকরণ যা overdispersion ধরতে পারে।

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

সমাধান ৬ (★)¶

(ক) linear predictor। $\hat\beta=[1.560,0.0597,0.301]$, temp$=25$, weekend$=1$: $$\eta=1.560+0.0597(25)+0.301(1)=1.560+1.4925+0.301=3.3535.$$

(খ) mean count। $$\mu=e^{\eta}=e^{3.3535}\approx28.60.$$

(গ) যাচাই ও variance। canonical $\hat\mu\approx28.6$ ✓ মিলেছে। Poisson equidispersion ধর্ম অনুযায়ী ওই দিনে $\operatorname{Var}(y)\approx\mu\approx28.6$ (SD $\approx5.35$)। লক্ষণীয় — log scale-এ যোগ, count scale-এ গুণ: $\mu=e^{1.560}\cdot e^{0.0597\cdot25}\cdot e^{0.301}=4.758\times4.444\times1.351\approx28.6$।

সমাধান ৭ (★)¶

(ক) per-$1$°C rate ratio। $e^{\hat\beta_1}=e^{0.0597}\approx1.0616$। → "তাপমাত্রা $1$°C বাড়লে (weekend স্থির) প্রত্যাশিত bike-count $\sim6.16\%$ বাড়ে।"

(খ) $+5$°C rate ratio। $e^{5\hat\beta_1}=e^{5\times0.0597}=e^{0.2985}\approx1.348$। → "$5$°C বাড়লে count $\sim34.8\%$ বাড়ে।" (যাচাই: rate ratio multiplicative, তাই $1.0616^5=1.348$ — একই উত্তর।)

(গ) weekend rate ratio। $e^{\hat\beta_2}=e^{0.301}\approx1.351$। → "weekend-এ (weekday-র তুলনায়, temp স্থির) প্রত্যাশিত count $\sim1.35$ গুণ — অর্থাৎ প্রায় $35\%$ বেশি।"

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

(ক) dispersion statistic। $$\hat\phi=\frac{\text{Pearson }\chi^2}{df}=\frac{1096.3}{247}=4.438\approx4.44\ \checkmark.$$ deviance দিয়েও: $1093.1/247=4.426$ — কাছাকাছি, দুটোই $\gg1$ (দুই পরিমাপ সামঞ্জস্যপূর্ণ)।

(খ) ব্যাখ্যা। $\hat\phi\approx1$ হলে equidispersion (Poisson যথেষ্ট)। এখানে $\hat\phi=4.44$ মানে আসল conditional variance Poisson-প্রত্যাশিত variance-এর $\sim4.4$ গুণ — শক্তিশালী overdispersion (rule of thumb: $\hat\phi>1.5$-ই উদ্বেগজনক, $4.44$ স্পষ্ট সমস্যা)।

(গ) বিকৃতি। Poisson SE তখন আসলের চেয়ে $\sim\sqrt{4.44}\approx2.107$ গুণ ছোট। তাই temp-এর Wald $z=\hat\beta_1/\mathrm{SE}_1$ প্রায় $2.1$ গুণ স্ফীত ⇒ over-confident inference: $p$-value কৃত্রিমভাবে ছোট, CI অতি-সংকীর্ণ, false-positive ঝুঁকি বাড়ে। (এখানে temp এত প্রবলভাবে significant যে সংশোধনের পরও থাকবে — কিন্তু দুর্বল predictor হলে এই স্ফীতি ভুল সিদ্ধান্ত দিতে পারত।)

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

(ক) quasi-Poisson SE। $\sqrt{\hat\phi}=\sqrt{4.44}\approx2.107$। Poisson SE $=[0.0476,0.00178,0.0307]$: $$\mathrm{SE}_{\text{quasi}}=2.107\times[0.0476,\,0.00178,\,0.0307]=[0.1003,\,0.00375,\,0.0647].$$

(খ) temp-এর $z$ তুলনা। $$z_{\text{Poisson}}=\frac{0.0597}{0.00178}\approx33.5,\qquad z_{\text{quasi}}=\frac{0.0597}{0.00375}\approx15.9.$$ $z$ প্রায় অর্ধেক হয়ে গেল ($33.5\to15.9$, ঠিক $\sqrt{4.44}$ ভাগ) — তবু $\gg1.96$, তাই temp এখনো অত্যন্ত significant। (এটাই দেখায়: overdispersion-সংশোধন significance "মুছে" দেয় না, কেবল over-confidence ছাঁটে।)

(গ) NB-র সাথে তুলনা। quasi SE $[0.100,\,0.0038,\,0.065]$ আর NB SE $[0.0874,\,0.00374,\,0.0689]$ — খুবই কাছাকাছি (বিশেষত temp ও weekend-এ প্রায় অভিন্ন)। দুই স্বাধীন পদ্ধতি (একটা quasi-likelihood, অন্যটা full NB likelihood) একই "SE অনেক বড় করতে হবে" সিদ্ধান্তে পৌঁছাচ্ছে — overdispersion-এর শক্ত নিশ্চয়তা ✓।

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

(ক) NB-implied variance। $\mu\approx28.6$, $\alpha=0.179$: $$\operatorname{Var}_{\text{NB}}=\mu+\alpha\mu^2=28.6+0.179\,(28.6)^2=28.6+0.179(817.96)=28.6+146.4\approx175.0.$$ Poisson-implied variance ছিল মাত্র $\mu=28.6$। অর্থাৎ NB ওই বিন্দুতে প্রায় $6$ গুণ বেশি variance অনুমোদন করে — যা data-র observed spread (var/mean $\approx10.5$)-এর সাথে অনেক বেশি সঙ্গতিপূর্ণ।

(খ) AIC তুলনা। $$\Delta\text{AIC}=\text{AIC}_{\text{Poisson}}-\text{AIC}_{\text{NB}}=2237.9-1753.3=484.6.$$ ছোট AIC ভালো → NB জেতে। $\Delta>10$-ই সাধারণত "strong evidence"; এখানে $\Delta\approx485$ — অপ্রতিরোধ্যভাবে decisive। AIC ($=-2\ell+2k$) NB-কে শাস্তি দেয় মাত্র $2\times1$ (একটা extra parameter $\alpha$), কিন্তু log-likelihood-উন্নতি তার তুলনায় বিশাল।

(গ) nested তাৎপর্য। NB ও Poisson nested ($\alpha=0$-তে NB $=$ Poisson)। মাত্র এক ডিগ্রি স্বাধীনতা খরচে ($\alpha$) AIC এত নামা মানে — overdispersion সত্যিকারের ও তীব্র, এবং NB-র quadratic variance তা ঠিকভাবে ধরছে। (চাইলে formal: likelihood-ratio test $H_0:\alpha=0$-ও এখানে অতি-significant হবে, কারণ $2\Delta\ell\approx486\gg\chi^2_{1}$।)

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

(ক) offset কেন coefficient $1$। rate-model করতে চাইলে শুরু করি $\log(\mu_i/t_i)=x_i^\top\beta$ থেকে। বাঁপাশ ভাঙলে: $$\log\mu_i-\log t_i=x_i^\top\beta\ \Longrightarrow\ \log\mu_i=\underbrace{\log t_i}_{\text{offset}}+x_i^\top\beta.$$ এখানে $\log t_i$ একটা জানা (data থেকেই পাওয়া, estimate-নয়) পদ যার coefficient জোর করে $1$-এ আটকানো — এ-ই offset। ফলে $x^\top\beta$ সরাসরি rate ($\mu_i/t_i$)-এর log মডেল করে, raw count নয়।

(খ) exposure দ্বিগুণ। একই $x_i$ কিন্তু $t=24$ বনাম $t=12$: যেহেতু rate স্থির, $\mu\propto t$, তাই $$\frac{\mu(t=24)}{\mu(t=12)}=\frac{24}{12}=2.$$ exposure দ্বিগুণ হলে প্রত্যাশিত count ঠিক দ্বিগুণ — স্বজ্ঞাতই সঠিক।

(গ) offset না দিলে। $t_i$-কে সাধারণ predictor হিসেবে রাখলে তার coefficient $1$-এর বদলে কিছু-একটা estimate হবে — তখন $\mu\propto t^{\hat\gamma}$ ($\hat\gamma\ne1$), অর্থাৎ exposure-proportionality নষ্ট ও rate-interpretation চলে যায়। তাছাড়া exposure-এর প্রভাব আর true predictor-গুলোর প্রভাব মিশে inference বিভ্রান্ত করে। offset ঠিক এই proportionality-কে জোর করে রাখে।

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

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

$Y\sim\text{Poisson}(\mu)$, $P(Y=k)=\dfrac{e^{-\mu}\mu^k}{k!}$, $k=0,1,2,\dots$।

(ক) $\mathbb E[Y]=\mu$। $$ \mathbb E[Y]=\sum_{k=0}^\infty k\,\frac{e^{-\mu}\mu^k}{k!} =\sum_{k=1}^\infty \frac{e^{-\mu}\mu^k}{(k-1)!} \overset{j=k-1}{=}\mu\sum_{j=0}^\infty\frac{e^{-\mu}\mu^{j}}{j!} =\mu\cdot 1=\mu, $$ কারণ ভেতরের যোগফল আবার সম্পূর্ণ Poisson pmf ($=1$)।

(খ) factorial moment ও variance। $$ \mathbb E[Y(Y-1)]=\sum_{k=2}^\infty k(k-1)\frac{e^{-\mu}\mu^k}{k!} =\sum_{k=2}^\infty\frac{e^{-\mu}\mu^k}{(k-2)!} \overset{j=k-2}{=}\mu^2\sum_{j=0}^\infty\frac{e^{-\mu}\mu^{j}}{j!}=\mu^2. $$ সুতরাং $\mathbb E[Y^2]=\mathbb E[Y(Y-1)]+\mathbb E[Y]=\mu^2+\mu$, এবং $$ \operatorname{Var}(Y)=\mathbb E[Y^2]-(\mathbb E[Y])^2=(\mu^2+\mu)-\mu^2=\mu.\qquad\blacksquare $$ অর্থাৎ $\mathbb E[Y]=\operatorname{Var}(Y)=\mu$ — equidispersion।

(গ) overdispersion-এর মূল। যেহেতু একটিমাত্র parameter $\mu$ একসাথে center ও spread ঠিক করে, variance-কে আলাদাভাবে বড়/ছোট করার কোনো স্বাধীনতা Poisson-এ নেই। তাই বাস্তব data-তে spread বেশি হলে (overdispersion) Poisson তা ধরতে অক্ষম — extra freedom আনতে quasi-Poisson-এর $\phi$ বা NB-র $\alpha$ লাগে।

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

log link: $\mu_i=e^{\eta_i}$, $\eta_i=x_i^\top\beta$। Poisson log-likelihood: $$\ell(\beta)=\sum_{i=1}^n\bigl[y_i\,x_i^\top\beta-e^{x_i^\top\beta}-\log(y_i!)\bigr].$$

পদ-পদ derivative ($\beta$-vector-এর সাপেক্ষে): - $\dfrac{\partial}{\partial\beta}\bigl[y_i\,x_i^\top\beta\bigr]=y_i\,x_i$ (রৈখিক পদ)। - $\dfrac{\partial}{\partial\beta}\bigl[e^{x_i^\top\beta}\bigr]=e^{x_i^\top\beta}\,x_i=\mu_i\,x_i$ (chain rule: $\frac{d}{d\eta}e^\eta=e^\eta=\mu$, তারপর $\frac{\partial\eta_i}{\partial\beta}=x_i$)। - $\dfrac{\partial}{\partial\beta}\bigl[\log(y_i!)\bigr]=0$ ($\beta$-মুক্ত ধ্রুবক)।

যোগ করলে: $$ \frac{\partial\ell}{\partial\beta}=\sum_{i=1}^n\bigl[y_i\,x_i-\mu_i\,x_i\bigr]=\sum_{i=1}^n(y_i-\mu_i)\,x_i=X^\top(y-\mu). $$ MLE-তে gradient শূন্য, তাই $\;X^\top(y-\mu)=\mathbf 0\;\blacksquare$। (এটি logistic-এর $X^\top(y-p)=0$-র হুবহু সমান্তরাল — canonical link-এর সাধারণ score-রূপ।)

conservation property। intercept column সব-$1$ ⇒ $X^\top(y-\mu)$-এর প্রথম উপাদান $\sum_i(y_i-\mu_i)=0$, অর্থাৎ $$\sum_{i=1}^n\hat\mu_i=\sum_{i=1}^n y_i.$$ fitted count-এর যোগফল ঠিক observed count-এর যোগফলের সমান — তাই sample mean সংরক্ষিত: $\frac1n\sum\hat\mu_i=\frac1n\sum y_i=19.56$। (logistic-এ এটাই ছিল $\sum\hat p_i=\sum y_i$।)

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

সমাধান ১৩-র score $\nabla\ell=\sum_i(y_i-\mu_i)x_i$। Hessian পেতে আবার $\beta^\top$-এর সাপেক্ষে derivative নিই। শুধু $\mu_i$ পদ $\beta$-নির্ভর, এবং $\dfrac{\partial\mu_i}{\partial\beta}=\mu_i x_i$ (যেহেতু $\mu_i=e^{x_i^\top\beta}$)। তাই: $$ \frac{\partial^2\ell}{\partial\beta\,\partial\beta^\top} =\frac{\partial}{\partial\beta^\top}\sum_i(y_i-\mu_i)x_i =-\sum_i\Bigl(\frac{\partial\mu_i}{\partial\beta^\top}\Bigr)x_i =-\sum_i\mu_i\,x_i x_i^\top=-X^\top W X, $$ যেখানে $W=\operatorname{diag}(\mu_1,\dots,\mu_n)$। $\blacksquare$

concavity। যেকোনো nonzero vector $v$-র জন্য $$ v^\top(X^\top W X)v=\sum_{i=1}^n\mu_i\,(x_i^\top v)^2\ge0, $$ কারণ প্রতিটি $\mu_i=e^{x_i^\top\beta}>0$ এবং $(x_i^\top v)^2\ge0$। তাই $X^\top WX\succeq0$ ⇒ Hessian $-X^\top WX\preceq0$ (negative semi-definite) ⇒ $\ell(\beta)$ concave। design $X$ full-rank ও কোনো degenerate পরিস্থিতি না-থাকলে $X^\top WX\succ0$ (strictly positive definite), তাই $\ell$ strictly concave ⇒ stationary point অনন্য ⇒ একক global MLE। concave function-এ স্থানীয়-সর্বোচ্চই বৈশ্বিক-সর্বোচ্চ, তাই Newton–Raphson/IRLS অবধারিতভাবে সেই global maximum-এ পৌঁছায় (logistic-এর মতোই local-minima-মুক্ত)।

IRLS সম্পর্ক। প্রতিটি Newton ধাপ গঠনগতভাবে একটা weighted least-squares solve, weight $w_i=\mu_i$ (logistic-এ ছিল $p_i(1-p_i)$); canonical (log) link হওয়ায় observed Hessian = expected Fisher information, তাই Newton ও Fisher scoring অভিন্ন হয়ে IRLS সরল হয়।

ঘ · কোডিং (Python)¶

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

নিচের একক runnable script ১–৮ সব ধাপ করে এবং প্রতিটি ফল canonical-এর সাথে মেলায়।

import numpy as np
import statsmodels.api as sm
from statsmodels.discrete.discrete_model import NegativeBinomial

# ---------- data (seed 20260619) ----------
rng = np.random.default_rng(20260619)
n = 250
temp    = rng.uniform(5, 35, n)
weekend = rng.binomial(1, 2/7, n)
mu_true = np.exp(1.5 + 0.06*temp + 0.4*weekend)
rate    = rng.gamma(6, 1/6, n)        # mean 1, overdispersion আনতে
count   = rng.poisson(mu_true * rate) # Gamma-Poisson mixture

X = sm.add_constant(np.column_stack([temp, weekend]))   # design matrix

# ---------- 1. EDA: var/mean ----------
print("count mean =", round(count.mean(), 2),          # ~19.56
      " var =", round(count.var(), 1),                  # ~205
      " var/mean =", round(count.var()/count.mean(), 2))# ~10.5  -> overdispersion সংকেত

# ---------- 2. Poisson fit ----------
pois = sm.GLM(count, X, family=sm.families.Poisson()).fit()
print("\nPoisson beta =", np.round(pois.params, 4))     # [1.560, 0.0597, 0.301]
print("Poisson SE   =", np.round(pois.bse, 5))          # [0.0476, 0.00178, 0.0307]
print("Poisson z    =", np.round(pois.tvalues, 2))
print("Poisson AIC  =", round(pois.aic, 1))             # 2237.9

# ---------- 3. rate ratios + 95% CI ----------
rr   = np.exp(pois.params)
ci   = np.exp(pois.conf_int())
print("\nrate ratios =", np.round(rr, 4))               # temp 1.0616, weekend 1.351
print("RR 95% CI:\n", np.round(ci, 4))
print("temp per +5C =", round(np.exp(5*pois.params[1]), 3))  # 1.348

# ---------- 4. overdispersion ----------
pearson = ((count - pois.mu)**2 / pois.mu).sum()
print("\nPearson chi2 =", round(pearson, 1),            # 1096.3
      " deviance =", round(pois.deviance, 1),            # 1093.1
      " df =", int(pois.df_resid))                       # 247
phi = pearson / pois.df_resid
print("dispersion phi =", round(phi, 2))                # 4.44  (>> 1)

# ---------- 5. quasi-Poisson SE = sqrt(phi) * Poisson SE ----------
quasi = sm.GLM(count, X, family=sm.families.Poisson()).fit(scale="X2")
print("\nquasi-Poisson SE =", np.round(quasi.bse, 5))   # ~[0.100, 0.0038, 0.065]
print("check sqrt(phi)*SE=", np.round(np.sqrt(phi)*pois.bse, 5))  # একই
print("temp z: Poisson", round(pois.tvalues[1], 1),     # ~33.5
      " quasi", round(quasi.tvalues[1], 1))             # ~15.9

# ---------- 6. negative binomial (alpha estimate) ----------
nb = NegativeBinomial(count, X).fit(disp=0)
alpha = nb.params[-1]            # শেষ param = alpha (lnalpha নয়, এই API-তে alpha)
print("\nNB beta  =", np.round(nb.params[:3], 4))       # [1.574, 0.0592, 0.300]
print("NB SE    =", np.round(nb.bse[:3], 5))            # [0.0874, 0.00374, 0.0689]
print("NB alpha =", round(alpha, 3))                    # 0.179
print("NB AIC   =", round(nb.aic, 1))                   # 1753.3

# ---------- 7. AIC তুলনা ----------
print("\nAIC: Poisson", round(pois.aic, 1),
      " NB", round(nb.aic, 1),
      " delta", round(pois.aic - nb.aic, 1))            # ~485 -> NB জেতে

# NB-implied variance at mu=28.6
mu_pt = 28.6
print("NB Var at mu=28.6 :", round(mu_pt + alpha*mu_pt**2, 1),  # ~175
      " vs Poisson Var", mu_pt)

# ---------- 8. prediction temp=25, weekend=1 ----------
x_new = np.array([1, 25, 1.0])
eta   = x_new @ pois.params
print("\nby-hand: eta =", round(eta, 4), " mu =", round(np.exp(eta), 2))  # 3.3535, 28.60
print("model.predict :", round(float(pois.predict(x_new[None, :])), 2))   # 28.60

# ---------- (ঐচ্ছিক) score-check: X^T (y - mu) ≈ 0 ----------
print("\nscore X^T(y-mu) =", np.round(X.T @ (count - pois.mu), 6))  # ≈ [0,0,0]

প্রত্যাশিত আউটপুট (canonical-এর সাথে মিলবে):

পরিমাণ	মান
count mean / var / ratio	$19.56$ / $\approx205$ / $\approx10.5$
Poisson $\hat\beta$	$[1.560,\,0.0597,\,0.301]$
Poisson SE	$[0.0476,\,0.00178,\,0.0307]$
rate ratio: temp / +5°C / weekend	$1.0616$ / $1.348$ / $1.351$
Pearson $\chi^2$ / deviance / $df$	$1096.3$ / $1093.1$ / $247$
dispersion $\hat\phi$	$4.44$
quasi-Poisson SE	$\approx[0.100,\,0.0038,\,0.065]$
NB $\hat\beta$ / SE / $\alpha$	$[1.574,0.0592,0.300]$ / $[0.0874,0.00374,0.0689]$ / $0.179$
AIC: Poisson / NB / $\Delta$	$2237.9$ / $1753.3$ / $\approx485$
NB Var @ $\mu=28.6$	$\approx175$ (vs Poisson $28.6$)
prediction (temp$=25$,weekend$=1$)	$\eta=3.3535,\ \hat\mu\approx28.6$
score $X^\top(y-\mu)$	$\approx\mathbf 0$ (§১৩ যাচাই)

মূল পাঠ: EDA-তেই var/mean $\approx10.5$ overdispersion-এর প্রথম ইঙ্গিত দেয়; Poisson fit চমৎকার $\hat\beta$ দিলেও dispersion $\hat\phi=4.44$ স্পষ্ট সমস্যা দেখায়; quasi-Poisson SE-কে $\sqrt{4.44}$ গুণ বাড়িয়ে ও NB তা explicit $\alpha=0.179$-এ ধরে — দুই পদ্ধতিই একই দিকে যায়, আর AIC ($\Delta\approx485$) decisively NB-কে বেছে নেয়। (statsmodels-এর সংস্করণভেদে NB API-তে $\alpha$ বনাম $\ln\alpha$-এর parametrization আলাদা হতে পারে; nb.params[-1] সরাসরি $\alpha\approx0.179$ দিলে ঠিক আ

পরিমাণ	মান
count mean / var / ratio	\(19.56\) / \(\approx205\) / \(\approx10.5\)
Poisson \(\hat\beta\)	\([1.560,\,0.0597,\,0.301]\)
Poisson SE	\([0.0476,\,0.00178,\,0.0307]\)
rate ratio: temp / +5°C / weekend	\(1.0616\) / \(1.348\) / \(1.351\)
Pearson \(\chi^2\) / deviance / \(df\)	\(1096.3\) / \(1093.1\) / \(247\)
dispersion \(\hat\phi\)	\(4.44\)
quasi-Poisson SE	\(\approx[0.100,\,0.0038,\,0.065]\)
NB \(\hat\beta\) / SE / \(\alpha\)	\([1.574,0.0592,0.300]\) / \([0.0874,0.00374,0.0689]\) / \(0.179\)
AIC: Poisson / NB / \(\Delta\)	\(2237.9\) / \(1753.3\) / \(\approx485\)
NB Var @ \(\mu=28.6\)	\(\approx175\) (vs Poisson \(28.6\))
prediction (temp\(=25\),weekend\(=1\))	\(\eta=3.3535,\ \hat\mu\approx28.6\)
score \(X^\top(y-\mu)\)	\(\approx\mathbf 0\) (§১৩ যাচাই)