최대가능도 추정¶

import pandas as pd
from sklearn.model_selection import train_test_split

def load_data(seed=1):
    """보험 가입 여부 자료를 읽어 훈련·시험으로 나눈다.

    설명변수는 나이 하나, 반응은 가입 여부(0/1)다. 자료가 27건뿐이라
    절반씩 나눈다.
    """
    url = ('https://raw.githubusercontent.com/codebasics/py/'
           '3ee4bde332ae7a10499c2679900094fa22ae191f/ML/7_logistic_reg/insurance_data.csv')
    df = pd.read_csv(url)
    x = df[['age']].values.reshape((-1, 1))
    y = df.bought_insurance.values.reshape((-1,))
    x_train, x_test, y_train, y_test = train_test_split(
        x, y, train_size=0.5, random_state=seed)
    return x_train, x_test, y_train, y_test


import numpy as np
import statsmodels.api as sm

x_tr, x_te, y_tr, y_te = load_data()
print(f"훈련 {len(y_tr)}건(가입 {y_tr.sum()}, 비율 {y_tr.mean():.4f}),  "
      f"시험 {len(y_te)}건(가입 {y_te.sum()}, 비율 {y_te.mean():.4f})")
order = np.argsort(x_tr.ravel())
print("나이순:", list(zip(x_tr.ravel()[order].tolist(), y_tr[order].tolist())))

# (25, 1) 하나를 빼면 완전분리가 된다.
keep = ~((x_tr.ravel() == 25) & (y_tr == 1))
x2, y2 = x_tr[keep], y_tr[keep]
A2 = sm.add_constant(x2)
for mi in (50, 500, 5000):
    res = sm.Logit(y2, A2).fit(disp=0, maxiter=mi)
    print(f"maxiter {mi:5d}  beta = {np.round(res.params, 4)}  "
          f"loglik = {res.llf:.3e}  수렴 {res.mle_retvals['converged']}")

훈련 13건(가입 8, 비율 0.6154),  시험 14건(가입 6, 비율 0.4286)
나이순: [(22, 0), (23, 0), (25, 0), (25, 1), (27, 0), (28, 0), (50, 1), (54, 1), (55, 1), (56, 1), (60, 1), (61, 1), (62, 1)]
maxiter    50  beta = [-183.0834    5.076 ]  loglik = 0.000e+00  수렴 False
maxiter   500  beta = [-183.0834    5.076 ]  loglik = 0.000e+00  수렴 False
maxiter  5000  beta = [-183.0834    5.0762]  loglik = 0.000e+00  수렴 False

import numpy as np

class LogisticRegression:
    """경사하강법으로 로지스틱 회귀를 직접 구현한다.

    최대가능도 추정이 무엇을 하는 일인지 보이려고 풀어 쓴 것이다. 선형회귀와
    달리 닫힌 해가 없어 수치적으로 찾아야 한다.
    """

    def __init__(self, x, y, lr=2e-4, epochs=100_000, theta=None):
        self.x = x
        self.y = y
        self.lr = lr
        self.epochs = epochs
        self.theta = (theta if theta is not None
                      else np.random.normal(size=(x.shape[1] + 1, 1)))

    @staticmethod
    def design_matrix(x):
        """1 로 채운 열을 앞에 붙여 절편을 만든다."""
        ones = np.ones((x.shape[0], 1))
        return np.concatenate((ones, x), axis=1)

    @staticmethod
    def sigmoid(z):
        """실수 전체를 (0, 1) 로 눌러 담는다. 확률로 읽을 수 있게 하는 장치다."""
        return 1 / (1 + np.exp(-z))

    def predict_proba(self, x):
        A = self.design_matrix(x)
        z = A @ self.theta
        return self.sigmoid(z).reshape((-1,))

    def predict(self, x):
        p = self.predict_proba(x)
        return (p > 0.5).astype(float)

    def loss(self):
        """음의 로그가능도(교차엔트로피). 이것을 최소화하는 것이 곧 MLE 다.

        eps 를 더하는 것은 p 가 정확히 0 이나 1 이 될 때 로그가 발산하는 것을
        막기 위함이다.
        """
        p = self.predict_proba(self.x)
        eps = 1e-6
        return -np.mean(
            self.y * np.log(p + eps) + (1 - self.y) * np.log(1 - p + eps))

    def gradient(self):
        """기울기. 놀랍게도 선형회귀와 똑같은 꼴인 X'(p - y) 로 나온다.

        시그모이드의 미분과 로그가능도의 미분이 서로 약분되면서 이렇게 된다.
        일반화선형모형 전체에서 되풀이되는 구조다.
        """
        A = self.design_matrix(self.x)
        p = self.predict_proba(self.x).reshape((-1, 1))
        y = self.y.reshape((-1, 1))
        return A.T @ (p - y)

    def train(self):
        for _ in range(self.epochs):
            self.theta -= self.lr * self.gradient()


# 헤세행렬 A^T W A 를 손으로 만들어 고유값을 본다.
def hessian(theta, x=x_tr):
    A = LogisticRegression.design_matrix(x)
    p = LogisticRegression.sigmoid(A @ theta).ravel()
    return A.T @ np.diag(p * (1 - p)) @ A

for name, th in [("초기값 (0, 0)", np.zeros((2, 1))),
                 ("최대가능도해", np.array([[-7.580328], [0.237100]]))]:
    ev = np.linalg.eigvalsh(hessian(th))
    print(f"{name:14s} 고유값 {ev[0]:.6f}, {ev[1]:.4f}  "
          f"조건수 {ev[1] / ev[0]:.1f}")

# loss() 와 gradient() 가 같은 함수의 값·기울기인가?
m = LogisticRegression(x_tr, y_tr, theta=np.array([[-7.0], [0.22]]))
eps = 1e-6
m.theta = np.array([[-7.0 + eps], [0.22]])
f1 = m.loss()
m.theta = np.array([[-7.0 - eps], [0.22]])
f0 = m.loss()
m.theta = np.array([[-7.0], [0.22]])
print(f"\nloss() 의 수치미분(절편) = {(f1 - f0) / (2 * eps):.6f}")
print(f"gradient() 의 첫 성분    = {m.gradient()[0, 0]:.6f}")
print(f"비 = {m.gradient()[0, 0] / ((f1 - f0) / (2 * eps)):.4f}   "
      f"(훈련 자료 수 n = {len(y_tr)})")

초기값 (0, 0)     고유값 0.418688, 6632.3313  조건수 15840.8
최대가능도해         고유값 0.034358, 631.4089  조건수 18377.3

loss() 의 수치미분(절편) = 0.008564
gradient() 의 첫 성분    = 0.111332
비 = 12.9996   (훈련 자료 수 n = 13)

# 직접 구현한 모형으로 적합한다. 학습률이 작고 반복이 10만 번이라 시간이 걸린다.
x_train, x_test, y_train, y_test = load_data()

# __init__ 이 theta 를 np.random.normal 로 뽑으므로 씨앗을 고정해야
# 결과를 되풀이할 수 있다.
np.random.seed(0)
model = LogisticRegression(x_train, y_train)
model.train()

y_pred = model.predict(x_test)
y_prob = model.predict_proba(x_test)

import statsmodels.api as sm

A = sm.add_constant(x_train)
mle = sm.Logit(y_train, A).fit(disp=0).params
nll = lambda t: -np.sum(
    y_train * np.log(1 / (1 + np.exp(-(A @ t))))
    + (1 - y_train) * np.log(1 - 1 / (1 + np.exp(-(A @ t)))))

print(f"10만 번 뒤 theta = {np.round(model.theta.ravel(), 6)}")
print(f"         NLL  = {nll(model.theta.ravel()):.6f}")
print(f"         기울기 = {np.round(model.gradient().ravel(), 6)}  "
      f"(0 이어야 한다)")

# 100만 번으로 늘려 본다.
np.random.seed(0)
long_model = LogisticRegression(x_train, y_train, epochs=1_000_000)
long_model.train()
print(f"\n100만 번 뒤 theta = {np.round(long_model.theta.ravel(), 6)}")
print(f"          NLL   = {nll(long_model.theta.ravel()):.6f}")

print(f"\n최대가능도해        = {np.round(mle, 6)}")
print(f"          NLL   = {nll(mle):.6f}")

# 씨앗을 바꾸면 10만 번 결과가 달라진다 — 미수렴의 증거다.
for s in (1, 2, 7):
    np.random.seed(s)
    m = LogisticRegression(x_train, y_train)
    m.train()
    print(f"씨앗 {s}: theta = {np.round(m.theta.ravel(), 6)}, "
          f"NLL = {nll(m.theta.ravel()):.6f}")

10만 번 뒤 theta = [-6.31561   0.193328]
         NLL  = 2.871060
         기울기 = [ 0.069215 -0.002281]  (0 이어야 한다)

100만 번 뒤 theta = [-7.578611  0.237038]
          NLL   = 2.833261

최대가능도해        = [-7.580328  0.2371  ]
          NLL   = 2.833261
씨앗 1: theta = [-6.319317  0.19345 ], NLL = 2.870803
씨앗 2: theta = [-6.393905  0.195917], NLL = 2.865875
씨앗 7: theta = [-6.317316  0.193384], NLL = 2.870942

# 같은 일을 sklearn 으로 하면 세 줄이면 된다. lbfgs 는 기울기뿐 아니라
# 곡률까지 어림해 쓰므로 경사하강법보다 훨씬 빨리 수렴한다.
from sklearn.linear_model import LogisticRegression

model = LogisticRegression(solver='lbfgs')
model.fit(x_train, y_train)

y_pred = model.predict(x_test)
y_prob = model.predict_proba(x_test)[:, 1]

no_pen = LogisticRegression(penalty=None, solver='lbfgs',
                            max_iter=10000).fit(x_train, y_train)

rows = [
    ("직접구현 10만",  model_theta := np.array([-6.315610, 0.193328])),
    ("직접구현 100만", np.array([-7.578611, 0.237038])),
    ("sklearn C=1",   np.r_[model.intercept_[0], model.coef_[0][0]]),
    ("penalty=None",  np.r_[no_pen.intercept_[0], no_pen.coef_[0][0]]),
]
print(f"{'방법':14s} {'절편':>10s} {'기울기':>9s} {'NLL':>9s} "
      f"{'벌점':>8s} {'목적함수':>9s}")
for name, t in rows:
    pen = 0.5 * t[1] ** 2
    print(f"{name:14s} {t[0]:10.5f} {t[1]:9.6f} {nll(t):9.6f} "
          f"{pen:8.5f} {nll(t) + pen:9.6f}")

print(f"\n결정경계   sklearn {-model.intercept_[0] / model.coef_[0][0]:.4f}   "
      f"무벌점 {-no_pen.intercept_[0] / no_pen.coef_[0][0]:.4f}")
print(f"시험 정확도 sklearn {model.score(x_test, y_test):.4f}   "
      f"무벌점 {no_pen.score(x_test, y_test):.4f}")
print(f"두 모형의 시험 예측이 모두 같은가: "
      f"{bool(np.all(model.predict(x_test) == no_pen.predict(x_test)))}")

방법                     절편       기울기       NLL       벌점      목적함수
직접구현 10만         -6.31561  0.193328  2.871060  0.01869  2.889748
직접구현 100만        -7.57861  0.237038  2.833261  0.02809  2.861355
sklearn C=1      -7.36199  0.228922  2.834167  0.02620  2.860369
penalty=None     -7.58041  0.237104  2.833261  0.02811  2.861370

결정경계   sklearn 32.1594   무벌점 31.9709
시험 정확도 sklearn 0.8571   무벌점 0.8571
두 모형의 시험 예측이 모두 같은가: True

from sklearn.linear_model import LogisticRegression
import numpy as np

X = np.array([[1], [2], [3], [4]])
y = np.array([0, 0, 1, 1])

model = LogisticRegression(penalty=None, solver='lbfgs')
model.fit(X, y)
print(f"beta0 = {model.intercept_[0]:.4f}")
print(f"beta1 = {model.coef_[0][0]:.4f}")

beta0 = -47.5544
beta1 = 18.8115

def loss(self):
    A = self.design_matrix(self.x)
    z = (A @ self.theta).reshape((-1,))
    # -log sigma(z) = softplus(-z),  -log(1 - sigma(z)) = softplus(z)
    return np.mean(np.logaddexp(0, z) - self.y * z)