সমাধান — অধ্যায় ৫.৪ · GLM: Logistic Regression¶

অধ্যায় ফাইল: part-5-modeling/05-04-glm-logistic-regression.md (§৭ অনুশীলনী)-এর সব প্রশ্নের পূর্ণ সমাধান। সব সংখ্যাগত ফল NumPy/SciPy দিয়ে যাচাই ও statsmodels/scikit-learn-এর সাথে মিলিয়ে দেখা হয়েছে; seeded কোড reproducible। চলমান উদাহরণ — পরীক্ষায় পাস-prediction, seed np.random.default_rng(20260619), $n=200$: predictor hours (uniform $0$–$10$) ও attend (uniform $40$–$100$); true $\eta=-6+0.9\cdot\text{hours}+0.04\cdot\text{attend}$, $p=\sigma(\eta)$, passed $\sim\text{Bernoulli}(p)$।

স্মারক — canonical সংখ্যা (seed 20260619): pass rate $0.62$; $\hat\beta=[-6.320,\,0.9224,\,0.04047]$ (intercept, hours, attend), $\mathrm{SE}=[1.173,\,0.1239,\,0.01278]$, Wald $z=[-5.39,\,7.44,\,3.17]$; odds ratio hours $e^{\hat\beta_1}=2.515$ (per $+1$ ঘণ্টা), attend per $+10$ point $e^{10\hat\beta_2}=1.499$; $\ell=-67.92$, null $\ell_0=-132.81$, McFadden pseudo-$R^2=0.489$, deviance $D=135.83$, LLR $p=6.5\times10^{-29}$; threshold $0.5$-এ accuracy $0.855$, confusion $\begin{bmatrix}\text{TN }61&\text{FP }15\\\text{FN }14&\text{TP }110\end{bmatrix}$, precision $0.88$, recall $0.887$, AUC $0.924$; example hours$=5$, attend$=75\Rightarrow\eta=1.327,\ p=0.790$। মূল সম্পর্ক: $p=\sigma(\eta)=\frac1{1+e^{-\eta}}$, $\operatorname{logit}(p)=\eta=x^\top\beta$, odds $=\frac{p}{1-p}=e^{\eta}$।

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

সমাধান ১ (★)¶

binary $y$-তে OLS (linear probability model) চালানোর তিনটি মূল সমস্যা:

(ক) পূর্বাভাস $[0,1]$-এর বাইরে যায়। OLS fitted মান $\hat y=x^\top\hat\beta$ যেকোনো বাস্তব সংখ্যা — তাই অনায়াসে $\hat y<0$ বা $\hat y>1$ পাওয়া যায়, যা "probability" হিসেবে অর্থহীন। চলমান data-তে অনেক বেশি hours/attend-এ OLS $\hat y>1$ দেবে। logistic regression-এর $\sigma(\cdot)$ পূর্বাভাসকে $(0,1)$-তে আটকে রাখে।

(খ) error-variance equal নয় (LINE-এর 'E' ভাঙে)। Bernoulli-তে $\operatorname{Var}(y\mid x)=p(1-p)$ — অর্থাৎ variance $p$-নির্ভর, $p=0.5$-এ সর্বোচ্চ, প্রান্তে কম। এটি heteroscedasticity: OLS-এর equal-variance (homoscedasticity) অনুমান ভাঙে, তাই SE/test invalid। উপরন্তু error দুটো মাত্র মান নেয় ($y-\hat y$), তাই Normality-ও ভাঙে।

(গ) রৈখিক প্রভাব saturate করে না। OLS-এ predictor বাড়লে প্রভাব অসীমভাবে রৈখিক বাড়ে; কিন্তু বাস্তবে hours $\to\infty$ গেলে পাসের সম্ভাবনা $1$-তে saturate হওয়া উচিত (আর $\to-\infty$-এ $0$-তে)। logistic curve ঠিক এই asymptotic আচরণ দেয়; সরলরেখা দেয় না।

সারমর্ম: binary outcome-এর জন্য চাই এমন মডেল যার পূর্বাভাস বৈধ probability, যার variance mean-নির্ভরতা স্বীকার করে, এবং যা প্রান্তে saturate করে — logistic regression তিনটাই দেয়।

সমাধান ২ (★)¶

(ক) logit ↔ sigmoid inverse। শুরু $\eta=\operatorname{logit}(p)=\log\frac{p}{1-p}$: $$e^{\eta}=\frac{p}{1-p}\;\Rightarrow\;e^{\eta}(1-p)=p\;\Rightarrow\;e^{\eta}=p(1+e^{\eta})\;\Rightarrow\;p=\frac{e^{\eta}}{1+e^{\eta}}=\frac{1}{1+e^{-\eta}}=\sigma(\eta).$$ অর্থাৎ $\sigma=\operatorname{logit}^{-1}$।

(খ) মান। $\eta=0\Rightarrow p=\frac1{1+1}=0.5$; $\eta\to+\infty\Rightarrow e^{-\eta}\to0\Rightarrow p\to1$; $\eta\to-\infty\Rightarrow e^{-\eta}\to\infty\Rightarrow p\to0$।

(গ) দুই স্কেল। মডেলটা log-odds স্কেলে linear (রৈখিক) — কারণ $\operatorname{logit}(p)=x^\top\beta$ হলো $x$-এ সরলরৈখিক; কিন্তু সেই একই সম্পর্ক probability স্কেলে S-আকৃতির, কারণ $\sigma(\cdot)$ log-odds-কে $(0,1)$-তে nonlinear (অরৈখিক)-ভাবে চেপে আনে। তাই coefficient log-odds-এ সরল additive, probability-তে নয়।

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

(ক) odds। $\text{odds}=\frac{p}{1-p}$ — event ঘটার সম্ভাবনা বনাম না-ঘটার সম্ভাবনার অনুপাত। উদাহরণ-পূর্বাভাসে $p=0.79\Rightarrow\text{odds}=\frac{0.79}{0.21}\approx3.76$, অর্থাৎ "পাসের পক্ষে প্রায় $3.76$ বনাম $1$"। (probability $0.79$ আর odds $3.76$ একই তথ্য, ভিন্ন স্কেল।)

(খ) odds ratio $e^{\hat\beta_1}=2.515$। logistic regression-এ predictor $x_j$ এক একক বাড়লে log-odds বাড়ে $\hat\beta_j$, তাই odds গুণ হয় $e^{\hat\beta_j}$। এখানে: hours $1$ ঘণ্টা বাড়লে (attend স্থির রেখে) পাসের odds $2.515$ গুণ হয় — অর্থাৎ $\sim151.5\%$ বৃদ্ধি। এটা probability-র ratio নয়, odds-এর ratio।

(গ) odds ratio-র চিহ্ন-পাঠ। $e^{\hat\beta}>1$ ⇒ predictor বাড়লে event-এর odds বাড়ে (positive প্রভাব); $=1$ ($\hat\beta=0$) ⇒ কোনো প্রভাব নেই; $<1$ ⇒ odds কমে (negative)। attend-এর $\hat\beta_2=0.04047$ ছোট দেখালেও significant, কারণ (i) এর SE আরও ছোট ($0.01278$), তাই $z=3.17>1.96$; এবং (ii) attend-এর scale বড় ($40$–$100$), তাই বাস্তব-মাপের পরিবর্তন (যেমন $+10$ point) odds $e^{10\times0.04047}=1.499$ গুণ করে — মোটেও তুচ্ছ নয়।

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

(ক) log-likelihood। $n$টা স্বাধীন Bernoulli observation, $p_i=\sigma(x_i^\top\beta)$: $$L(\beta)=\prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i},\qquad \ell(\beta)=\sum_{i=1}^n\bigl[y_i\log p_i+(1-y_i)\log(1-p_i)\bigr].$$ $\log p_i=\eta_i-\log(1+e^{\eta_i})$ আর $\log(1-p_i)=-\log(1+e^{\eta_i})$ বসিয়ে সরল রূপ: $$\ell(\beta)=\sum_{i=1}^n\bigl[y_i\eta_i-\log(1+e^{\eta_i})\bigr],\qquad \eta_i=x_i^\top\beta.$$

(খ) closed form নেই কেন। maximize করতে $\nabla\ell=X^\top(y-p)=0$ সমাধান করতে হয় (সমাধান ১৩), কিন্তু $p_i=\sigma(x_i^\top\beta)$ হলো $\beta$-র nonlinear (অরৈখিক) function — তাই $\beta$-কে বীজগণিতে আলাদা করে আনা যায় না (OLS-এ $X^\top(y-X\beta)=0$ রৈখিক বলে $\hat\beta=(X^\top X)^{-1}X^\top y$ পাওয়া যেত; এখানে $\sigma$-র কারণে তা সম্ভব নয়)।

(গ) iterative সমাধান ও WLS-সম্পর্ক। ব্যবহৃত হয় Newton–Raphson: $$\beta^{(t+1)}=\beta^{(t)}+(X^\top W^{(t)}X)^{-1}X^\top(y-p^{(t)}),\qquad W=\operatorname{diag}(p_i(1-p_i)).$$ প্রতিটি ধাপ পুনর্লিখন করলে এটি একটা weighted least squares solve হয়ে ওঠে (working response-এর উপর), তাই নাম IRLS — iteratively reweighted least squares। weight $w_i=p_i(1-p_i)$ ঠিক Bernoulli variance, যা প্রতি iteration-এ $\hat p_i$-র সাথে আপডেট হয়। (concavity-র কারণে — সমাধান ১৪ — এটি সবসময় global maximum-এ পৌঁছায়।)

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

logistic regression-কে GLM-এর তিন উপাদানে ভাঙা:

উপাদান	logistic regression	ordinary linear regression (৫.১)
(ক) random component (response distribution)	$y_i\sim\text{Bernoulli}(p_i)$ — exponential family-র সদস্য	$y_i\sim\mathcal N(\mu_i,\sigma^2)$
(খ) systematic component (linear predictor)	$\eta_i=x_i^\top\beta$	$\eta_i=x_i^\top\beta$
(গ) link function $g(\text{mean})=\eta$	$g(p)=\operatorname{logit}(p)=\log\frac{p}{1-p}$ (canonical link)	$g(\mu)=\mu$ (identity link)

তিনটি ক্ষেত্রেই systematic component এক ($\eta=x^\top\beta$); শুধু distribution ও link বদলায়। ordinary linear regression হলো GLM-এর সেই special case যেখানে distribution Normal এবং link identity ($g(\mu)=\mu\Rightarrow\mathbb E[y]=x^\top\beta$ সরাসরি)। logistic regression Bernoulli + logit বেছে নেয় বলে পূর্বাভাস বৈধ probability হয় এবং mean-নির্ভর variance-কে স্বাভাবিকভাবে ধরে। (৫.৫-এ একই ছাঁচে Poisson + log link বসবে।)

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

সমাধান ৬ (★)¶

$\hat\beta=[-6.320,\,0.9224,\,0.04047]$, hours$=5$, attend$=75$।

(ক) linear predictor: $$\eta=-6.320+0.9224(5)+0.04047(75)=-6.320+4.612+3.035=\mathbf{1.327}.$$

(খ) probability: $$p=\sigma(1.327)=\frac{1}{1+e^{-1.327}}=\frac{1}{1+0.2653}=\frac{1}{1.2653}\approx\mathbf{0.790}.$$

(গ) odds, এবং $e^{\eta}$-এর সাথে মিল: $$\text{odds}=\frac{p}{1-p}=\frac{0.790}{0.210}\approx3.76,\qquad e^{\eta}=e^{1.327}\approx3.77.$$ দুটো (rounding বাদে) সমান — কারণ সংজ্ঞা-অনুসারে $\frac{p}{1-p}=e^{\eta}$ যখন $p=\sigma(\eta)$ ✓। canonical $\eta=1.327,\ p=0.790$ হুবহু মিলেছে।

সমাধান ৭ (★)¶

(ক) hours-এর odds ratio: $e^{\hat\beta_1}=e^{0.9224}\approx\mathbf{2.515}$। ব্যাখ্যা: "attend স্থির রেখে অধ্যয়ন $1$ ঘণ্টা বাড়লে পাসের odds $2.515$ গুণ হয়" ($\sim151\%$ বৃদ্ধি)। ✓ canonical।

(খ) attend-এর $10$-point odds ratio: coefficient per-unit, তাই $10$ point-এ exponent $10\hat\beta_2$: $$e^{10\times0.04047}=e^{0.4047}\approx\mathbf{1.499}.$$ "উপস্থিতি $10$ percentage-point বাড়লে পাসের odds $\sim1.5$ গুণ" ✓।

(গ) $2$ ঘণ্টা বেশি: log-odds additive, তাই odds ratio multiplicative — $$e^{2\hat\beta_1}=\bigl(e^{\hat\beta_1}\bigr)^2=2.515^2\approx\mathbf{6.33}\ \text{গুণ}.$$ ($k$ একক বাড়লে odds $e^{k\hat\beta}$ গুণ — এই log-linearity-ই odds-ratio interpretation-কে এত পরিষ্কার করে।)

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

$\hat\beta=[-6.320,\,0.9224,\,0.04047]$, $\mathrm{SE}=[1.173,\,0.1239,\,0.01278]$।

(ক) Wald $z=\hat\beta_j/\mathrm{SE}_j$: $$z_0=\frac{-6.320}{1.173}\approx-5.39,\quad z_1=\frac{0.9224}{0.1239}\approx7.44,\quad z_2=\frac{0.04047}{0.01278}\approx3.17.$$ canonical $[-5.39,\,7.44,\,3.17]$ মিলেছে ✓।

(খ) significance ($\lvert z\rvert>1.96$): তিনটাই $\lvert z\rvert>1.96$ — intercept, hours, attend সবই $\alpha=0.05$-এ significant।

(গ) hours-এর CI ও odds-ratio CI: $$\hat\beta_1\pm1.96\,\mathrm{SE}_1=0.9224\pm1.96(0.1239)=0.9224\pm0.2429=[0.6795,\ 1.1653].$$ odds-ratio স্কেলে exponentiate: $$\bigl[e^{0.6795},\ e^{1.1653}\bigr]=[1.97,\ 3.21].$$ পুরো interval $>1$ (অর্থাৎ $1$-কে অন্তর্ভুক্ত করে না), তাই hours-এর প্রভাব $\alpha=0.05$-এ significant ও positive। কেন $1$ গুরুত্বপূর্ণ: odds ratio $=1\Leftrightarrow\hat\beta=0\Leftrightarrow$ "কোনো প্রভাব নেই"; CI যদি $1$-কে ঘিরে ফেলত, তবে প্রভাব null থেকে আলাদা বলা যেত না (Wald test-এর সাথে সম্পূর্ণ সঙ্গতিপূর্ণ)।

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

$\ell=-67.92$, $\ell_0=-132.81$।

(ক) deviance। residual deviance $D=-2\ell=-2(-67.92)=\mathbf{135.83}$ ✓ (canonical)। null deviance $D_0=-2\ell_0=-2(-132.81)=265.63$।

(খ) likelihood-ratio test। $$G^2=2(\ell-\ell_0)=2\bigl(-67.92-(-132.81)\bigr)=2(64.89)=\mathbf{129.78}=D_0-D.$$ দুই predictor (hours, attend), তাই $G^2\sim\chi^2_2$ under null। $\chi^2_{0.95,2}=5.99$, আর $129.78\gg5.99$ — null (দুই coefficient-ই শূন্য) জোরালোভাবে বাতিল; মডেল intercept-only-এর চেয়ে অত্যন্ত ভালো (canonical LLR $p=6.5\times10^{-29}$)।

(গ) McFadden pseudo-$R^2$। $$R^2_{\text{McF}}=1-\frac{\ell}{\ell_0}=1-\frac{-67.92}{-132.81}=1-0.5114=\mathbf{0.489}\ \checkmark.$$ OLS-$R^2$ থেকে পার্থক্য: এটি ব্যাখ্যাকৃত variance-এর অংশ নয়; বরং null মডেলের তুলনায় log-likelihood কতটা উন্নত হলো তার আপেক্ষিক পরিমাপ। তাই $0.489$-কে "৪৯% variance ব্যাখ্যা" বলা ভুল — McFadden pseudo-$R^2$ সাধারণত OLS-$R^2$-এর চেয়ে ছোট দেখায়, আর $0.2$–$0.4$-ই প্রায়ই "ভালো fit" ধরা হয়।

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

confusion (threshold $0.5$, $n=200$): TN $=61$, FP $=15$, FN $=14$, TP $=110$।

(ক) accuracy: $$\text{accuracy}=\frac{\mathrm{TP}+\mathrm{TN}}{n}=\frac{110+61}{200}=\frac{171}{200}=\mathbf{0.855}\ \checkmark.$$

(খ) precision ও recall (positive = "pass"): $$\text{precision}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}=\frac{110}{110+15}=\frac{110}{125}=\mathbf{0.880},$$ $$\text{recall}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}=\frac{110}{110+14}=\frac{110}{124}\approx\mathbf{0.887}.$$ অর্থাৎ "pass" বলা case-এর $88\%$ সত্যিই pass (precision), আর সত্যিকারের pass-দের $88.7\%$ ধরা পড়েছে (recall)।

(গ) FP বেশি ক্ষতিকর হলে threshold। FP = আসলে-fail-কে "pass" বলা। এটি কমাতে threshold বাড়াতে হবে (যেমন $0.5\to0.7$) — তখন কেবল উঁচু-$\hat p$ case-কেই "pass" বলা হবে, FP কমবে ⇒ precision বাড়বে। কিন্তু এতে কিছু সত্যিকারের pass-ও বাদ পড়ে FN বাড়ে ⇒ recall কমবে। এটাই precision–recall trade-off; সঠিক threshold দুই ভুলের আপেক্ষিক খরচের উপর নির্ভর করে (ROC/PR curve দেখে বাছাই করা হয়)।

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

AUC $=0.924$।

(ক) ROC-এর অক্ষ। threshold $0\to1$ ঘোরালে — $y$-অক্ষে TPR (true positive rate $=$ recall $=$ sensitivity $=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$) বনাম $x$-অক্ষে FPR (false positive rate $=1-$specificity $=\frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}$)।

(খ) AUC-র অর্থ। AUC $=0.924$ মানে: "একটা random সত্যিকারের-positive case-কে মডেল একটা random সত্যিকারের-negative case-এর চেয়ে বেশি $\hat p$ score দেওয়ার সম্ভাবনা $92.4\%$"। $0.5$ = random guessing, $1.0$ = perfect ranking।

(গ) threshold-নিরপেক্ষতা ও imbalance। accuracy একটা নির্দিষ্ট cutoff-এ গোনা হয় বলে threshold বদলালেই বদলায় (threshold-নির্ভর); AUC সব threshold জুড়ে TPR-FPR curve-এর নিচের পুরো ক্ষেত্রফল মাপে, তাই কোনো একক cutoff-এর উপর নির্ভর করে না। imbalance: majority class বড় হলে (যেমন $95\%$ negative) শুধু সবাইকে "negative" বলেই accuracy $95\%$ পাওয়া যায় — অথচ মডেল আসলে কিছুই শেখেনি; AUC তখনও ranking-গুণ মাপে বলে এই ফাঁদে পড়ে না, তাই imbalanced সমস্যায় বেশি নির্ভরযোগ্য। ($0.924$ এখানে শক্তিশালী discrimination নির্দেশ করে।)

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

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

দাবি: $\sigma'(z)=\sigma(z)\bigl(1-\sigma(z)\bigr)$, যেখানে $\sigma(z)=\dfrac{1}{1+e^{-z}}$।

প্রমাণ। $\sigma(z)=(1+e^{-z})^{-1}$। chain rule: $$\sigma'(z)=-(1+e^{-z})^{-2}\cdot\frac{d}{dz}(1+e^{-z})=-(1+e^{-z})^{-2}\cdot(-e^{-z})=\frac{e^{-z}}{(1+e^{-z})^2}.$$ এবার ভাঙি: $$\frac{e^{-z}}{(1+e^{-z})^2}=\underbrace{\frac{1}{1+e^{-z}}}_{=\,\sigma(z)}\cdot\underbrace{\frac{e^{-z}}{1+e^{-z}}}_{=\,1-\sigma(z)}.$$ শেষ ধাপে $\dfrac{e^{-z}}{1+e^{-z}}=\dfrac{(1+e^{-z})-1}{1+e^{-z}}=1-\sigma(z)$। তাই $$\boxed{\sigma'(z)=\sigma(z)(1-\sigma(z))}.\qquad\blacksquare$$

কেন চাবিকাঠি। $p=\sigma(\eta)$ হওয়ায় $\dfrac{\partial p}{\partial\eta}=p(1-p)$ — ঠিক Bernoulli variance $\operatorname{Var}(y\mid x)=p(1-p)$। এটিই (i) score-derivation-এ $\log(1+e^{\eta})$-এর derivative $=p$ হওয়ার মূল (সমাধান ১৩), এবং (ii) Hessian/IRLS-এর weight $w_i=p_i(1-p_i)$ (সমাধান ১৪)। অর্থাৎ একটা সরল identity logistic regression-এর গোটা MLE-যন্ত্রকে চালায়।

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

দাবি: $\dfrac{\partial\ell}{\partial\beta}=\displaystyle\sum_{i=1}^n(y_i-p_i)x_i=X^\top(y-p)$, এবং MLE-তে এটি $=\mathbf 0$।

প্রমাণ। $\ell(\beta)=\sum_i\bigl[y_i\eta_i-\log(1+e^{\eta_i})\bigr]$, $\eta_i=x_i^\top\beta$, তাই $\dfrac{\partial\eta_i}{\partial\beta}=x_i$।

পদ ১: $\dfrac{\partial}{\partial\beta}\bigl[y_i\eta_i\bigr]=y_i x_i$।

পদ ২: $\dfrac{\partial}{\partial\beta}\log(1+e^{\eta_i})=\dfrac{1}{1+e^{\eta_i}}\cdot e^{\eta_i}\cdot x_i=\dfrac{e^{\eta_i}}{1+e^{\eta_i}}\,x_i=\sigma(\eta_i)\,x_i=p_i x_i$।

(এখানে $\dfrac{e^{\eta}}{1+e^{\eta}}=\sigma(\eta)$ — সমাধান ১২-র identity-র সরাসরি অনুসিদ্ধান্ত।) যোগ করে: $$\frac{\partial\ell}{\partial\beta}=\sum_{i=1}^n\bigl(y_i-p_i\bigr)x_i.$$ matrix-আকারে, $X$-এর $i$-তম সারি $x_i^\top$, $y=(y_i)$, $p=(p_i)$ হলে $$\frac{\partial\ell}{\partial\beta}=X^\top(y-p).$$ MLE $\hat\beta$ stationary point, তাই $X^\top(y-\hat p)=\mathbf 0$। $\blacksquare$

ফলিত অর্থ (conservation)। intercept থাকায় $X$-এর একটা column সব-$1$; সেই সারির সমীকরণ: $$\sum_{i=1}^n(y_i-\hat p_i)=0\;\Rightarrow\;\sum_i\hat p_i=\sum_i y_i.$$ অর্থাৎ fitted probability-র যোগফল ঠিক observed positive-সংখ্যার সমান। চলমান data-তে $\sum_i y_i=0.62\times200=124$, তাই $\sum_i\hat p_i=124$ — মডেল মোট পাস-সংখ্যা হুবহু পুনরুৎপাদন করে। (OLS-এ residual-যোগফল $0$ হওয়ার সরাসরি logistic-সমতুল্য।)

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

দাবি: Hessian $H=\dfrac{\partial^2\ell}{\partial\beta\,\partial\beta^\top}=-X^\top W X$, $W=\operatorname{diag}(p_i(1-p_i))$, এবং $\ell$ concave।

প্রমাণ। সমাধান ১৩ থেকে $\nabla\ell=\sum_i(y_i-p_i)x_i$। $\beta^\top$-র সাপেক্ষে আবার derivative; শুধু $p_i$ $\beta$-নির্ভর, আর সমাধান ১২ দিল $\dfrac{\partial p_i}{\partial\beta^\top}=p_i(1-p_i)\,x_i^\top$: $$H=\frac{\partial}{\partial\beta^\top}\sum_i(y_i-p_i)x_i=-\sum_i x_i\,\frac{\partial p_i}{\partial\beta^\top}=-\sum_i p_i(1-p_i)\,x_i x_i^\top=-X^\top W X.$$

negative semi-definite। যেকোনো $v\in\mathbb R^{d+1}$-এর জন্য $$v^\top(X^\top W X)v=\sum_i p_i(1-p_i)\,(x_i^\top v)^2\ge0,$$ কারণ প্রতিটি $p_i\in(0,1)\Rightarrow p_i(1-p_i)>0$ এবং বর্গ $\ge0$। তাই $H=-X^\top WX\preceq0$ ⇒ $\ell$ concave। যদি $X$ full column rank এবং কোনো $x_i^\top v$ সব $i$-তে $0$ না হয় (no perfect separation), তবে $v^\top X^\top WX\,v>0$ for $v\ne0$ ⇒ strictly concave।

ফল। concave function-এর যেকোনো stationary point-ই global maximum; strictly concave হলে তা একক। তাই logistic-regression-এ MLE অনন্য (OLS-এর মতোই — local-minima সমস্যা নেই), এবং Newton–Raphson/IRLS অবধারিতভাবে সেই global maximum-এ converge করে। $\blacksquare$ (সতর্কতা: classes perfectly separable হলে likelihood $1$-এর দিকে monotone বাড়ে, MLE অসীমে চলে যায় — তখন regularization/penalized fit লাগে।)

ঘ · কোডিং (Python)¶

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

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

import numpy as np
import statsmodels.api as sm
from sklearn.metrics import confusion_matrix, roc_auc_score, roc_curve

# ---------- ডেটা (canonical seed) ----------
rng = np.random.default_rng(20260619)
n = 200
hours  = rng.uniform(0, 10, n)
attend = rng.uniform(40, 100, n)
eta_true = -6 + 0.9*hours + 0.04*attend
p_true   = 1/(1 + np.exp(-eta_true))
passed   = rng.binomial(1, p_true)
print(f"pass rate = {passed.mean():.4f}")          # 0.62

X = sm.add_constant(np.column_stack([hours, attend]))   # [1, hours, attend]
names = ["const", "hours", "attend"]

# ---------- 1) fit ----------
model = sm.Logit(passed, X).fit(disp=0)
b, se, z = model.params, model.bse, model.tvalues
print("\n--- coefficient summary ---")
for j, nm in enumerate(names):
    print(f"{nm:7s}  beta={b[j]:+.4f}  SE={se[j]:.4f}  z={z[j]:+.2f}  p={model.pvalues[j]:.2e}")
# beta = [-6.3204, 0.9224, 0.0405], z = [-5.39, 7.44, 3.17]

# ---------- 2) odds ratios ----------
print("\n--- odds ratios (exp beta) + 95% CI ---")
ci = model.conf_int()
for j, nm in enumerate(names):
    print(f"{nm:7s}  OR={np.exp(b[j]):.4f}  CI=[{np.exp(ci[j,0]):.3f}, {np.exp(ci[j,1]):.3f}]")
print(f"hours per +1 hr : OR = {np.exp(b[1]):.4f}")          # 2.515
print(f"attend per +10  : OR = {np.exp(10*b[2]):.4f}")        # 1.499

# ---------- 3) goodness-of-fit ----------
D  = -2*model.llf
print("\n--- goodness of fit ---")
print(f"loglik={model.llf:.2f}  null_loglik={model.llnull:.2f}")   # -67.92, -132.81
print(f"deviance={D:.2f}  pseudoR2(McFadden)={model.prsquared:.4f}")  # 135.83, 0.489
print(f"LLR stat={model.llr:.2f}  LLR p={model.llr_pvalue:.2e}")     # 129.78, 6.5e-29

# ---------- 4) classification + ROC ----------
phat = model.predict(X)
pred = (phat >= 0.5).astype(int)
cm = confusion_matrix(passed, pred)          # [[TN, FP], [FN, TP]]
(TN, FP), (FN, TP) = cm
acc  = (TP + TN)/n
prec = TP/(TP + FP)
rec  = TP/(TP + FN)
auc  = roc_auc_score(passed, phat)
print("\n--- classification @ threshold 0.5 ---")
print("confusion [[TN,FP],[FN,TP]] =", cm.tolist())   # [[61,15],[14,110]]
print(f"accuracy={acc:.4f}  precision={prec:.4f}  recall={rec:.4f}  AUC={auc:.4f}")
# 0.855, 0.880, 0.887, 0.924
fpr, tpr, thr = roc_curve(passed, phat)      # ROC curve পয়েন্ট (প্লট করা যায়)

# ---------- 5) by-hand prediction check ----------
eta5 = b[0] + b[1]*5 + b[2]*75
p5   = 1/(1 + np.exp(-eta5))
print("\n--- by-hand: hours=5, attend=75 ---")
print(f"eta={eta5:.4f}  p(hand)={p5:.4f}  p(model)={model.predict([[1,5,75]])[0]:.4f}")
# eta=1.327, p=0.790

# ঐচ্ছিক: score-equation যাচাই (সমাধান ১৩)
print("\nX^T (y - phat) =", np.round(X.T @ (passed - phat), 8))   # ~ [0,0,0]
print("sum(phat) =", round(phat.sum(), 4), " sum(y) =", passed.sum())  # 124, 124

প্রত্যাশিত output (canonical-এর সাথে মিল):

pass rate = 0.6200

--- coefficient summary ---
const    beta=-6.3204  SE=1.1732  z=-5.39  p=7.1e-08
hours    beta=+0.9224  SE=0.1239  z=+7.44  p=1.0e-13
attend   beta=+0.0405  SE=0.0128  z=+3.17  p=1.6e-03

--- odds ratios (exp beta) + 95% CI ---
hours    OR=2.5154  CI=[1.973, 3.207]
attend   OR=1.0413  CI=[1.016, 1.068]
hours per +1 hr : OR = 2.5154
attend per +10  : OR = 1.4988

--- goodness of fit ---
loglik=-67.92  null_loglik=-132.81
deviance=135.83  pseudoR2(McFadden)=0.4886
LLR stat=129.78  LLR p=6.54e-29

--- classification @ threshold 0.5 ---
confusion [[TN,FP],[FN,TP]] = [[61, 15], [14, 110]]
accuracy=0.8550  precision=0.8800  recall=0.8871  AUC=0.9236

--- by-hand: hours=5, attend=75 ---
eta=1.3266  p(model)=0.7903  ...

X^T (y - phat) = [ 0.  0.  0.]
sum(phat) = 124.0  sum(y) = 124

প্রতিটি মান উপরের canonical তালিকার সাথে মিলেছে: $\hat\beta=[-6.320,0.9224,0.04047]$, $z=[-5.39,7.44,3.17]$, OR hours $2.515$, attend per-$10$ $1.499$, $\ell=-67.92$, $D=135.83$, pseudo-$R^2=0.489$, LLR $p=6.5\times10^{-29}$, confusion $[[61,15],[14,110]]$, accuracy $0.855$, AUC $0.924$, by-hand $\eta=1.327,\ p=0.790$। score-যাচাইয়ে $X^\top(y-\hat p)=\mathbf 0$ ও $\sum\hat p=\sum y=124$ — সমাধান ১৩-র তাত্ত্বিক ফল কো

উপাদান	logistic regression	ordinary linear regression (৫.১)
(ক) random component (response distribution)	\(y_i\sim\text{Bernoulli}(p_i)\) — exponential family-র সদস্য	\(y_i\sim\mathcal N(\mu_i,\sigma^2)\)
(খ) systematic component (linear predictor)	\(\eta_i=x_i^\top\beta\)	\(\eta_i=x_i^\top\beta\)
(গ) link function \(g(\text{mean})=\eta\)	\(g(p)=\operatorname{logit}(p)=\log\frac{p}{1-p}\) (canonical link)	\(g(\mu)=\mu\) (identity link)