From 9142a28f5285e54d24cd2bc958d089bb01d5cc5f Mon Sep 17 00:00:00 2001
From: suda-yuga <yuga274274@gmail.com>
Date: Wed, 25 Jun 2025 18:20:24 +0900
Subject: [PATCH 1/3] Replace np.sum with @ for vector inner products, fix
 typos, and improve numerical stability using scipy.special.gammaln

---
 lectures/mle.md | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/lectures/mle.md b/lectures/mle.md
index ca9701b6d..343d984bb 100644
--- a/lectures/mle.md
+++ b/lectures/mle.md
@@ -334,7 +334,7 @@ $$
         = &
         \sum_{i=1}^{n} y_i \log{\mu_i} -
         \sum_{i=1}^{n} \mu_i -
-        \sum_{i=1}^{n} \log y!
+        \sum_{i=1}^{n} \log y_i!
 \end{split}
 $$
 
@@ -344,7 +344,7 @@ $$
 \underset{\beta}{\max} \Big(
 \sum_{i=1}^{n} y_i \log{\mu_i} -
 \sum_{i=1}^{n} \mu_i -
-\sum_{i=1}^{n} \log y! \Big)
+\sum_{i=1}^{n} \log y_i! \Big)
 $$
 
 However, no analytical solution exists to the above problem -- to find the MLE
@@ -458,7 +458,7 @@ class PoissonRegression:
     def logL(self):
         y = self.y
         μ = self.μ()
-        return np.sum(y * np.log(μ) - μ - np.log(factorial(y)))
+        return np.sum(y * np.log(μ) - μ - scipy.special.gammaln(y + 1))
 
     def G(self):
         y = self.y
@@ -869,13 +869,12 @@ class ProbitRegression:
 
     def logL(self):
         μ = self.μ()
-        return np.sum(y * np.log(μ) + (1 - y) * np.log(1 - μ))
+        return y @ np.log(μ) + (1 - y) @ np.log(1 - μ)
 
     def G(self):
         μ = self.μ()
         ϕ = self.ϕ()
-        return np.sum((X.T * y * ϕ / μ - X.T * (1 - y) * ϕ / (1 - μ)),
-                     axis=1)
+        return X.T @ (y * ϕ / μ - (1 - y) * ϕ / (1 - μ))
 
     def H(self):
         X = self.X

From a0c5032bf57b63b517219a978545bc23f3c58c00 Mon Sep 17 00:00:00 2001
From: suda-yuga <yuga274274@gmail.com>
Date: Wed, 25 Jun 2025 20:18:59 +0900
Subject: [PATCH 2/3] Add missing import of gammaln from scipy.special

---
 lectures/mle.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/lectures/mle.md b/lectures/mle.md
index 343d984bb..97191c28c 100644
--- a/lectures/mle.md
+++ b/lectures/mle.md
@@ -47,7 +47,7 @@ import matplotlib.pyplot as plt
 plt.rcParams["figure.figsize"] = (11, 5)  #set default figure size
 import numpy as np
 from numpy import exp
-from scipy.special import factorial
+from scipy.special import factorial, gammaln
 import pandas as pd
 from mpl_toolkits.mplot3d import Axes3D
 import statsmodels.api as sm
@@ -458,7 +458,7 @@ class PoissonRegression:
     def logL(self):
         y = self.y
         μ = self.μ()
-        return np.sum(y * np.log(μ) - μ - scipy.special.gammaln(y + 1))
+        return np.sum(y * np.log(μ) - μ - gammaln(y + 1))
 
     def G(self):
         y = self.y

From 01fbe3efd25045f37363644fb729cce91caffe87 Mon Sep 17 00:00:00 2001
From: suda-yuga <yuga274274@gmail.com>
Date: Mon, 30 Jun 2025 21:17:35 +0900
Subject: [PATCH 3/3] Add self. to instance variable references.md

---
 lectures/mle.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/lectures/mle.md b/lectures/mle.md
index 97191c28c..8031a9295 100644
--- a/lectures/mle.md
+++ b/lectures/mle.md
@@ -868,16 +868,20 @@ class ProbitRegression:
         return norm.pdf(self.X @ self.β.T)
 
     def logL(self):
+        y = self.y
         μ = self.μ()
         return y @ np.log(μ) + (1 - y) @ np.log(1 - μ)
 
     def G(self):
+        X = self.X
+        y = self.y  
         μ = self.μ()
         ϕ = self.ϕ()
         return X.T @ (y * ϕ / μ - (1 - y) * ϕ / (1 - μ))
 
     def H(self):
         X = self.X
+        y = self.y
         β = self.β
         μ = self.μ()
         ϕ = self.ϕ()