1 이 포스트의 위치
로지스틱 회귀: 모형에서 GLM 3요소, 로지스틱 함수 성질, MLE 스코어 방정식, 피셔 정보행렬을 다루었다. 이 포스트는 “스코어 방정식을 어떻게 수치적으로 푸는가”에 집중한다.
핵심 질문: 비선형 스코어 방정식을 왜 반복 가중 최소제곱(IRLS)으로 풀 수 있는가?
순서: 로그 우도의 순오목성 증명 → Newton-Raphson 유도 → IRLS 재표현 → 다중 변수 확장 → 이탈도 → 수렴 진단
2 로그 우도의 순오목성
2.1 헤시안 행렬 유도
이전 포스트에서 스코어 벡터를 유도하였다. MLE의 유일성을 보장하려면 로그 우도가 순오목(strictly concave)임을 확인해야 한다. 이를 위해 헤시안 행렬(Hessian matrix)을 계산한다.
로지스틱 링크 \(F(w) = e^w/(1+e^w)\) 에서 \(f_i/[F_i(1-F_i)] = 1\) 이므로 스코어는:
\[\frac{\partial \log L}{\partial \alpha} = \sum_{i=1}^n (y_i - \pi_i), \qquad \frac{\partial \log L}{\partial \beta} = \sum_{i=1}^n (y_i - \pi_i) x_i\]
여기서 \(\pi_i = F(\alpha + \beta x_i)\).
2차 도함수를 계산한다. \(\partial \pi_i / \partial \alpha = \pi_i(1-\pi_i)\), \(\partial \pi_i / \partial \beta = \pi_i(1-\pi_i)x_i\) 이므로:
\[\frac{\partial^2 \log L}{\partial \alpha^2} = -\sum_{i=1}^n \pi_i(1-\pi_i)\]
\[\frac{\partial^2 \log L}{\partial \alpha \partial \beta} = -\sum_{i=1}^n \pi_i(1-\pi_i) x_i\]
\[\frac{\partial^2 \log L}{\partial \beta^2} = -\sum_{i=1}^n \pi_i(1-\pi_i) x_i^2\]
헤시안 행렬:
\[H(\alpha, \beta) = \begin{pmatrix} \dfrac{\partial^2 \log L}{\partial \alpha^2} & \dfrac{\partial^2 \log L}{\partial \alpha \partial \beta} \\[6pt] \dfrac{\partial^2 \log L}{\partial \beta \partial \alpha} & \dfrac{\partial^2 \log L}{\partial \beta^2} \end{pmatrix} = -\begin{pmatrix} \sum w_i & \sum w_i x_i \\ \sum w_i x_i & \sum w_i x_i^2 \end{pmatrix}\]
여기서 로지스틱 가중치 \(w_i = \pi_i(1-\pi_i) > 0\) 이다.
\(-H = \mathbf{X}^T \mathbf{W} \mathbf{X}\) 로 쓸 수 있는데, 이것이 바로 피셔 정보행렬 \(I(\alpha, \beta)\) 와 같다. 즉 관측 정보행렬 = 기대 정보행렬 (로지스틱 링크의 특수 성질).
2.2 순오목성 증명
로지스틱 회귀의 로그 우도 \(\log L(\alpha, \beta | \mathbf{y})\) 는 \((\alpha, \beta) \in \mathbb{R}^2\) 에서 순오목이다.
즉, 모든 영 벡터가 아닌 \(\mathbf{v} = (v_1, v_2)^T\) 에 대해 \(\mathbf{v}^T H \mathbf{v} < 0\) 이다.
증명:
\[\mathbf{v}^T H \mathbf{v} = -\mathbf{v}^T (\mathbf{X}^T \mathbf{W} \mathbf{X}) \mathbf{v} = -\sum_{i=1}^n w_i (v_1 + v_2 x_i)^2\]
\(w_i = \pi_i(1-\pi_i) > 0\) 이고, \(v_1 + v_2 x_i = 0\) 이 모든 \(i\) 에 대해 동시에 성립하려면 \(x_i\) 들이 모두 같아야 한다 (완전 공선). 예측변수가 상수가 아니면 이 합은 음수이다.
따라서 \(x_i\) 들이 모두 같지 않으면 헤시안은 음정치(negative definite)이고 로그 우도는 순오목이다.
귀결: 스코어 방정식의 해가 존재하면 그 해는 유일하며, 이것이 전역 MLE이다 (Casella & Berger, 2002, §12.3.2). 완전 분리(complete separation) 상황에서는 해가 존재하지 않지만, 이 경우 로지스틱 모형이 데이터에 적합하지 않다는 신호이다.
3 Newton-Raphson 알고리즘
3.1 기본 원리
비선형 방정식 \(\nabla \log L(\boldsymbol{\theta}) = \mathbf{0}\) 을 풀기 위해 현재 추정값 \(\boldsymbol{\theta}^{(t)} = (\alpha^{(t)}, \beta^{(t)})^T\) 에서 테일러 전개를 사용한다:
\[\nabla \log L(\boldsymbol{\theta}) \approx \nabla \log L(\boldsymbol{\theta}^{(t)}) + H(\boldsymbol{\theta}^{(t)})(\boldsymbol{\theta} - \boldsymbol{\theta}^{(t)})\]
이를 \(\mathbf{0}\) 으로 설정하고 \(\boldsymbol{\theta}\) 에 대해 풀면 Newton-Raphson 업데이트:
\[\boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} - [H(\boldsymbol{\theta}^{(t)})]^{-1} \nabla \log L(\boldsymbol{\theta}^{(t)})\]
직관: \(-[H]^{-1}\) 은 “얼마나 어느 방향으로 움직일지”를 결정한다. \(\nabla \log L\) 은 우도가 증가하는 방향(기울기)을 알려주고, \(-H^{-1}\) 은 그 방향의 곡률(curvature)로 스텝 크기를 조절한다. 곡률이 크면(우도 곡면이 가파르면) 스텝이 작아지고, 곡률이 작으면(완만하면) 스텝이 크다. 1차원 최대화 문제 \(\ell''(\theta^{(t)}) < 0\) 일 때 \(\theta^{(t+1)} = \theta^{(t)} - \ell'(\theta^{(t)})/\ell''(\theta^{(t)})\) 와 완전히 동일한 구조이다.
로지스틱 회귀에서 \(H = -\mathbf{X}^T\mathbf{W}\mathbf{X}\) 이므로:
\[\boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} + [\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X}]^{-1} \mathbf{X}^T (\mathbf{y} - \boldsymbol{\pi}^{(t)})\]
여기서 \(\mathbf{W}^{(t)} = \mathrm{diag}(\pi_i^{(t)}(1-\pi_i^{(t)}))\) 은 현재 추정값에서의 가중행렬이다.
3.2 Fisher Scoring과의 동치
Fisher scoring은 실제 헤시안 대신 기대 헤시안(피셔 정보행렬의 음수)을 사용하는 변형이다:
\[\boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} + [I(\boldsymbol{\theta}^{(t)})]^{-1} \nabla \log L(\boldsymbol{\theta}^{(t)})\]
로지스틱 링크에서는 \(-H = I(\boldsymbol{\theta})\) 이므로 (관측 정보 = 기대 정보) Newton-Raphson과 Fisher scoring이 동치이다 (McCullagh & Nelder, 1989, §2.5). 이것은 정준 링크(canonical link)의 특수 성질이다.
4 IRLS: 반복 재가중 최소제곱
4.1 조정 반응변수 (Working Response)
Newton-Raphson 업데이트를 WLS 문제로 변환하면 계산 구조가 명확해진다.
\(t\) 번째 반복에서 조정 반응변수 \(z_i^{(t)}\) 를 다음과 같이 정의한다 (McCullagh & Nelder, 1989, §2.5):
\[z_i^{(t)} = \hat{\eta}_i^{(t)} + (y_i - \hat{\pi}_i^{(t)}) \cdot \frac{d\eta}{d\mu}\bigg|_{\hat{\pi}_i^{(t)}}\]
로지스틱 링크에서 \(\eta = \log[\pi/(1-\pi)]\) 이므로:
\[\frac{d\eta}{d\mu} = \frac{d}{d\pi}\log\frac{\pi}{1-\pi} = \frac{1}{\pi} + \frac{1}{1-\pi} = \frac{1}{\pi(1-\pi)}\]
따라서:
\[z_i^{(t)} = \hat{\eta}_i^{(t)} + \frac{y_i - \hat{\pi}_i^{(t)}}{\hat{\pi}_i^{(t)}(1-\hat{\pi}_i^{(t)})}\]
\(z_i^{(t)}\) 는 “현재 선형 예측자에서, 잔차 \((y_i - \hat{\pi}_i)\) 를 링크 도함수로 스케일한 값을 더한 것”이다. 이것은 링크 함수를 데이터 \(y_i\) 에 1차 테일러 전개한 것과 같다 (McCullagh & Nelder, 1989, §2.5):
\[g(y_i) \approx g(\hat{\mu}_i) + (y_i - \hat{\mu}_i)g'(\hat{\mu}_i) = \hat{\eta}_i + (y_i - \hat{\pi}_i)\frac{d\eta}{d\mu}\]
4.2 IRLS 업데이트 공식
가중행렬 \(\mathbf{W}^{(t)} = \mathrm{diag}(w_i^{(t)})\) 에서 \(w_i^{(t)} = \hat{\pi}_i^{(t)}(1-\hat{\pi}_i^{(t)})\) 이고, 조정 반응변수 벡터 \(\mathbf{z}^{(t)} = (z_1^{(t)}, \ldots, z_n^{(t)})^T\) 를 정의하면 Newton-Raphson 업데이트는 아래와 같이 WLS 방정식과 동치이다 (McCullagh & Nelder, 1989, §4.4):
\[\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X}\, \hat{\boldsymbol{\beta}}^{(t+1)} = \mathbf{X}^T \mathbf{W}^{(t)} \mathbf{z}^{(t)}\]
명시적 해:
\[\hat{\boldsymbol{\beta}}^{(t+1)} = (\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W}^{(t)} \mathbf{z}^{(t)}\]
이것이 바로 IRLS이다.
동치 증명: Newton-Raphson 업데이트를 WLS 형태로 변환한다.
\[\hat{\boldsymbol{\beta}}^{(t+1)} = \hat{\boldsymbol{\beta}}^{(t)} + (\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})^{-1}\mathbf{X}^T(\mathbf{y} - \hat{\boldsymbol{\pi}}^{(t)})\]
양변에 \((\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})\) 를 곱하면:
\[(\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})\hat{\boldsymbol{\beta}}^{(t+1)} = (\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})\hat{\boldsymbol{\beta}}^{(t)} + \mathbf{X}^T(\mathbf{y} - \hat{\boldsymbol{\pi}}^{(t)})\]
우변:
\[= \mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X}\hat{\boldsymbol{\beta}}^{(t)} + \mathbf{X}^T\mathbf{W}^{(t)}\mathbf{W}^{(t)^{-1}}(\mathbf{y} - \hat{\boldsymbol{\pi}}^{(t)})\]
\[= \mathbf{X}^T\mathbf{W}^{(t)}\left[\mathbf{X}\hat{\boldsymbol{\beta}}^{(t)} + \mathbf{W}^{(t)^{-1}}(\mathbf{y} - \hat{\boldsymbol{\pi}}^{(t)})\right]\]
\[= \mathbf{X}^T\mathbf{W}^{(t)}\left[\hat{\boldsymbol{\eta}}^{(t)} + \frac{\mathbf{y} - \hat{\boldsymbol{\pi}}^{(t)}}{w^{(t)}}\right]\]
\[= \mathbf{X}^T\mathbf{W}^{(t)}\mathbf{z}^{(t)} \quad \blacksquare\]
4.3 IRLS 알고리즘 절차
초기화: \(\hat{\boldsymbol{\pi}}^{(0)} = \mathbf{y} + \epsilon\) (또는 전체 성공률 \(\bar{y}\)로 상수 초기화)
각 반복 \(t = 0, 1, 2, \ldots\) 에서:
- 선형 예측자: \(\hat{\eta}_i^{(t)} = \log[\hat{\pi}_i^{(t)}/(1-\hat{\pi}_i^{(t)})]\)
- 가중치: \(w_i^{(t)} = \hat{\pi}_i^{(t)}(1-\hat{\pi}_i^{(t)})\)
- 조정 반응: \(z_i^{(t)} = \hat{\eta}_i^{(t)} + (y_i - \hat{\pi}_i^{(t)})/w_i^{(t)}\)
- WLS 풀기: \(\hat{\boldsymbol{\beta}}^{(t+1)} = (\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{z}^{(t)}\)
- 새 확률: \(\hat{\pi}_i^{(t+1)} = \text{logistic}(\mathbf{x}_i^T\hat{\boldsymbol{\beta}}^{(t+1)})\)
- 수렴 확인: \(\|\hat{\boldsymbol{\beta}}^{(t+1)} - \hat{\boldsymbol{\beta}}^{(t)}\| < \epsilon_{\text{tol}}\)
수렴 후: \(\hat{\boldsymbol{\beta}} = \hat{\boldsymbol{\beta}}^{(t+1)}\)
IRLS의 장점: 각 반복이 표준 WLS 문제이므로 QR 분해나 Cholesky 분해 등 안정적인 선형대수 루틴을 재사용할 수 있다.
5 다중 로지스틱 회귀로의 확장
단일 예측변수 \(x_i\) 를 \(p\) 개 예측변수 벡터 \(\mathbf{x}_i = (1, x_{i1}, \ldots, x_{i,p-1})^T\) 로 확장한다.
모형:
\[\log\!\left(\frac{\pi_i}{1-\pi_i}\right) = \mathbf{x}_i^T \boldsymbol{\beta}, \quad \boldsymbol{\beta} = (\beta_0, \beta_1, \ldots, \beta_{p-1})^T\]
로그 우도, 스코어, 헤시안 모두 벡터 형태로 동일하게 성립한다:
| 항목 | 단변수 표현 | 행렬 표현 |
|---|---|---|
| 선형 예측자 | \(\hat{\eta}_i = \alpha + \beta x_i\) | \(\hat{\boldsymbol{\eta}} = \mathbf{X}\boldsymbol{\beta}\) |
| 스코어 | \(\sum(y_i - \pi_i)\), \(\sum(y_i-\pi_i)x_i\) | \(\mathbf{X}^T(\mathbf{y} - \boldsymbol{\pi})\) |
| 헤시안 (음수) | \(\begin{pmatrix}\sum w_i & \sum w_i x_i \\ \sum w_i x_i & \sum w_i x_i^2\end{pmatrix}\) | \(\mathbf{X}^T\mathbf{W}\mathbf{X}\) |
| 정보행렬 | \(I(\alpha,\beta)\) | \(I(\boldsymbol{\beta}) = \mathbf{X}^T\mathbf{W}\mathbf{X}\) |
| IRLS 업데이트 | \(\hat{\boldsymbol{\beta}}^{(t+1)} = (\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{z}^{(t)}\) | 동일 |
MLE의 점근 분포:
\[\sqrt{n}(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}) \xrightarrow{d} N(\mathbf{0},\; I(\boldsymbol{\beta})^{-1})\]
각 계수의 표준오차는 \(I(\hat{\boldsymbol{\beta}})^{-1}\) 의 대각 원소의 제곱근이다.
모수 해석 (다중 변수):
- \(e^{\beta_k}\): 다른 모든 예측변수를 고정했을 때, \(x_k\) 가 1 단위 증가할 때의 조건부 오즈비(conditional odds ratio)
- 이것이 로지스틱 회귀의 핵심 해석이다 — 다른 변수 통제 후의 효과
6 이탈도 (Deviance)
IRLS로 MLE \(\hat{\boldsymbol{\beta}}\) 를 구한 후, “이 모형이 데이터를 얼마나 잘 설명하는가?”를 측정해야 한다. 잔차 제곱합이 선형 회귀의 적합도 척도라면, 로지스틱 회귀에서는 이탈도(deviance)가 그 역할을 한다. 이탈도는 모형 비교(LRT), 적합도 검정, 예측변수 선택에 모두 사용된다.
6.1 정의
이탈도(deviance)는 포화 모형(saturated model, 각 \(\pi_i = y_i\))과 적합 모형의 로그 우도 차이의 2배로 정의된다 (McCullagh & Nelder, 1989, §4.4):
\[D(\mathbf{y}; \hat{\boldsymbol{\pi}}) = 2\left[\ell(\tilde{\boldsymbol{\pi}}; \mathbf{y}) - \ell(\hat{\boldsymbol{\pi}}; \mathbf{y})\right]\]
여기서 \(\tilde{\pi}_i = y_i/m_i\) (포화 모형의 MLE), \(m_i\) 는 \(x_i\) 에서의 관측 수이다.
포화 모형을 기준으로 쓰는 이유: 포화 모형은 “각 관측점에 하나의 모수”를 할당하여 데이터를 완벽하게 적합한다 — 즉 달성 가능한 최대 로그 우도이다. 이탈도는 “우리 모형이 이 완벽한 적합에서 얼마나 떨어져 있는가”를 측정한다. 절대적인 로그 우도 값은 해석하기 어렵지만, 포화 모형과의 차이(이탈도)는 “설명하지 못한 정보”의 양을 나타내어 모형 간 비교가 가능하다.
이항(binomial) 케이스에서 명시적으로:
\[D(\mathbf{y}; \hat{\boldsymbol{\pi}}) = 2\sum_{j=1}^J \left\{y_j^* \log\!\left(\frac{y_j^*}{\hat{\mu}_j}\right) + (m_j - y_j^*)\log\!\left(\frac{m_j - y_j^*}{m_j - \hat{\mu}_j}\right)\right\}\]
여기서 \(\hat{\mu}_j = m_j \hat{\pi}(x_j)\) 이다.
이 공식이 말하는 것: 각 항 \(y_j^*\log(y_j^*/\hat{\mu}_j)\) 는 관측된 성공 횟수 \(y_j^*\) 와 모형 예측 성공 기대치 \(\hat{\mu}_j\) 의 로그비 가중합이다. 이것은 쿨백-라이블러(KL) 발산의 이항 형태이다 — 관측 분포와 적합 분포가 얼마나 다른지를 측정한다. 관측과 예측이 완벽히 일치하면(\(y_j^* = \hat{\mu}_j\)) 각 항이 0이 되어 \(D = 0\), 예측이 나쁠수록 \(D\) 가 커진다. 인수 2를 곱하는 것은 LRT의 \(\chi^2\) 극한 분포를 맞추기 위해서이다.
6.2 이탈도와 LRT의 관계
두 모형의 이탈도 차이는 우도비 검정 통계량이다:
\[-2\log\lambda = D(\mathbf{y}; \hat{\boldsymbol{\pi}}_0) - D(\mathbf{y}; \hat{\boldsymbol{\pi}})\]
귀무 모형이 \(p_0\) 개 모수, 대안 모형이 \(p_1\) 개 모수이면 (\(p_1 > p_0\)):
\[D(\mathbf{y}; \hat{\boldsymbol{\pi}}_0) - D(\mathbf{y}; \hat{\boldsymbol{\pi}}) \xrightarrow{d} \chi^2_{p_1 - p_0}\]
이것이 LRT의 점근 분포이다.
특수 케이스: \(H_0: \beta = 0\) 검정에서 \(p_0 = 1\) (절편만), \(p_1 = 2\) (절편+기울기):
\[-2\log\lambda = 2[\ell(\hat{\alpha},\hat{\beta}|\mathbf{y}) - \ell(\hat{\alpha}_0,0|\mathbf{y})] \xrightarrow{d} \chi^2_1\]
이것이 Casella & Berger(2002) §12.3.2의 LRT 통계량과 동일하다.
6.3 이탈도 기반 적합도 검정
\(J\) 개 고유 \(x\) 값에서 \(m_j\) 개씩 관측한 경우, 적합 모형 자체의 이탈도가 적합도 통계량으로 사용된다:
\[D(\mathbf{y}; \hat{\boldsymbol{\pi}}) \xrightarrow{d} \chi^2_{J - p}\]
이는 \(m_j \to \infty\) 인 극한에서 성립한다. \(m_j\) 가 작은 경우(개별 이진 관측) 이 근사는 좋지 않다. 이유: \(\chi^2_{J-p}\) 근사는 각 셀의 기대 빈도 \(m_j\hat{\pi}_j\) 가 충분히 커야(통상 \(\geq 5\)) 성립한다. 개별 이진 관측(\(m_j = 1\))에서는 \(\hat{\mu}_j = \hat{\pi}(x_j) \in (0,1)\) 이므로 기대 빈도가 항상 1 미만이고, 이탈도의 표본 분포가 \(\chi^2\) 와 크게 다르다. 이 경우에는 Hosmer-Lemeshow 검정이나 잔차 분석 등 대안적 적합도 진단을 사용한다.
7 수렴 진단
7.1 초기값 선택
McCullagh & Nelder(1989, §4.4)가 권장하는 초기값:
\[\tilde{\mu}_i = \frac{y_i + 0.5}{m_i + 1}\]
이것은 순수한 \(y_i/m_i\) 보다 안전한데, \(y_i = 0\) 이나 \(y_i = m_i\) 일 때 \(\log(0)\) 또는 \(\log(\infty)\) 계산을 피하기 때문이다.
7.2 수렴 기준
두 가지 기준을 병행한다:
- 계수 변화량: \(\|\hat{\boldsymbol{\beta}}^{(t+1)} - \hat{\boldsymbol{\beta}}^{(t)}\|_2 < 10^{-6}\)
- 이탈도 변화량: \(|D^{(t+1)} - D^{(t)}| < 10^{-8}\)
완전 분리 상황에서는 계수가 발산해도 이탈도와 적합 확률 \(\hat{\boldsymbol{\pi}}\) 는 수렴하는 경우가 있다. 따라서 계수 기준만으로는 부족하고 두 기준 모두 확인해야 한다 (McCullagh & Nelder, 1989, §4.4).
이유: 완전 분리 시 \(Y=0\) 군은 \(\hat{\pi}_i \to 0\), \(Y=1\) 군은 \(\hat{\pi}_i \to 1\) 로 수렴한다. 이 극한에서 로그 우도 \(\sum y_i\log\hat{\pi}_i + (1-y_i)\log(1-\hat{\pi}_i) \to 0\) 으로 수렴하고 이탈도도 일정 값에 수렴한다. 그러나 \(\hat{\pi} = \text{logistic}(\mathbf{x}^T\boldsymbol{\beta})\) 가 0 또는 1이 되려면 \(\|\boldsymbol{\beta}\| \to \infty\) 이어야 하므로 계수 자체는 발산한다. 즉 “극한의 적합확률”은 잘 정의되어도 그것을 생성하는 “계수”는 유한하지 않다.
7.3 수렴 실패의 징후
| 증상 | 원인 | 대처 |
|---|---|---|
| \(\|\hat{\boldsymbol{\beta}}\|\) 발산 | 완전 분리 | Firth 교정, 릿지 정규화 |
| 이탈도 수렴, 계수 발산 | 부분 분리 | 계수 해석 불가, 예측값만 신뢰 |
| 수렴 느림 | 예측변수 스케일 불균형 | 표준화 후 재실행 |
| 이탈도 진동 | 초기값 문제 | 다른 초기값 시도 |
대표 시나리오: 이진 분류에서 예측변수 하나가 \(Y=0\)과 \(Y=1\)을 완벽히 분리(예: 검사 점수 50점 미만은 전부 실패, 이상은 전부 합격)하면, 로지스틱 모형은 임계 점수 부근에서 기울기를 무한히 크게 만들려 하므로 \(\|\hat{\boldsymbol{\beta}}\|\) 가 반복마다 증가한다. 이때 IRLS를 50회 이상 돌려도 수렴 판정이 나지 않고, summary() 에서 계수의 표준오차가 매우 크거나 nan으로 나온다. Firth 교정은 로그 우도에 페널티 \(\frac{1}{2}\log\det I(\boldsymbol{\beta})\) 를 더해 완전 분리에서도 유한한 MLE를 보장한다.
8 코드 예시
8.1 Step 1: IRLS 직접 구현
import numpy as np
def logistic(eta):
"""수치적으로 안정한 시그모이드"""
return np.where(eta >= 0,
1 / (1 + np.exp(-eta)),
np.exp(eta) / (1 + np.exp(eta)))
def irls_logistic(X, y, tol=1e-8, max_iter=100, verbose=False):
"""
IRLS로 로지스틱 회귀 MLE를 계산한다.
X: (n, p) 설계 행렬 (절편 열 포함)
y: (n,) 이진 반응 벡터
반환: (beta, n_iter, converged)
"""
n, p = X.shape
# 초기화: μ_0 = (y + 0.5) / (m + 1) 권장 (McCullagh & Nelder, §4.4)
pi = np.clip(y + 0.5, 0.01, 0.99) / (1 + 1)
pi = np.full(n, y.mean().clip(0.01, 0.99))
beta = np.zeros(p)
for t in range(max_iter):
# 1. 선형 예측자
eta = X @ beta
# 2. 적합 확률
pi = logistic(eta)
# 3. 가중치
w = pi * (1 - pi)
W = np.diag(w)
# 4. 조정 반응변수: z_i = η̂_i + (y_i - π̂_i) / w_i
z = eta + (y - pi) / w
# 5. WLS: β^(t+1) = (X^T W X)^{-1} X^T W z
XtWX = X.T @ W @ X
XtWz = X.T @ (w * z)
beta_new = np.linalg.solve(XtWX, XtWz)
# 6. 수렴 확인
delta = np.linalg.norm(beta_new - beta)
if verbose:
ll = np.sum(y * np.log(pi + 1e-15) + (1-y) * np.log(1 - pi + 1e-15))
print(f" iter {t+1}: ‖Δβ‖ = {delta:.2e}, log L = {ll:.4f}")
beta = beta_new
if delta < tol:
return beta, t + 1, True
return beta, max_iter, False
# ──────────────────────────────────────
# 챌린저 데이터 적용
# ──────────────────────────────────────
temp = np.array([53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70,
70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81], dtype=float)
failure = np.array([1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0], dtype=float)
n = len(temp)
X = np.column_stack([np.ones(n), temp]) # 설계 행렬 (절편 + 기온)
beta_hat, n_iter, converged = irls_logistic(X, failure, verbose=True)
print(f"\n=== IRLS 결과 ===")
print(f"α̂ = {beta_hat[0]:.4f} (Casella: 15.043)")
print(f"β̂ = {beta_hat[1]:.4f} (Casella: -0.232)")
print(f"수렴: {converged}, 반복 횟수: {n_iter}")
# 31°F 예측
pi_31 = logistic(beta_hat[0] + beta_hat[1] * 31)
print(f"\n31°F에서 실패 확률: {pi_31:.4f}")
# 피셔 정보행렬과 표준오차
pi_hat = logistic(X @ beta_hat)
w_hat = pi_hat * (1 - pi_hat)
I_hat = X.T @ np.diag(w_hat) @ X
I_inv = np.linalg.inv(I_hat)
se_beta = np.sqrt(np.diag(I_inv))
print(f"\nse(α̂) = {se_beta[0]:.4f}")
print(f"se(β̂) = {se_beta[1]:.4f}")
print(f"\n95% CI for β: [{beta_hat[1]-1.96*se_beta[1]:.4f}, {beta_hat[1]+1.96*se_beta[1]:.4f}]")
# 이탈도 계산
ll_fit = np.sum(failure * np.log(pi_hat + 1e-15) + (1-failure) * np.log(1 - pi_hat + 1e-15))
ll_sat = np.sum(failure * np.log(failure + 1e-15) + (1-failure) * np.log(1 - failure + 1e-15))
deviance = 2 * (ll_sat - ll_fit)
print(f"\n이탈도 D = {deviance:.4f}")결과 해석: IRLS는 통상 5–10회 반복 내에 수렴한다 (converged=True, n_iter≈6). \(\hat{\alpha} \approx 15.04\), \(\hat{\beta} \approx -0.232\) 가 Casella 예시와 일치함을 확인한다. 이탈도 \(D \approx 20.3\) 은 절편만인 귀무 모형의 이탈도와 비교했을 때 기온 변수가 설명하는 “정보 감소량”이다. 계수가 수렴한 직후 스코어 벡터가 \(\approx (0, 0)\) 임을 함께 확인하면 MLE 조건이 충족되었음을 이중 검증할 수 있다.
8.2 Step 2: IRLS 반복 과정 시각화 + statsmodels 검증
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy.stats import chi2
# ──────────────────────────────────────
# IRLS 반복 추적
# ──────────────────────────────────────
def logistic(eta):
return np.where(eta >= 0,
1 / (1 + np.exp(-eta)),
np.exp(eta) / (1 + np.exp(eta)))
temp = np.array([53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70,
70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81], dtype=float)
failure = np.array([1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0], dtype=float)
n = len(temp)
X = np.column_stack([np.ones(n), temp])
# IRLS 반복 과정 추적
beta = np.zeros(2)
beta_trace = [beta.copy()]
ll_trace = []
for t in range(30):
eta = X @ beta
pi = logistic(eta)
w = pi * (1 - pi)
z = eta + (y - pi) / np.maximum(w, 1e-10) if False else eta + (failure - pi) / np.maximum(w, 1e-10)
XtWX = X.T @ np.diag(w) @ X
XtWz = X.T @ (w * z)
beta_new = np.linalg.solve(XtWX, XtWz)
ll = np.sum(failure * np.log(logistic(X @ beta_new) + 1e-15)
+ (1-failure) * np.log(1 - logistic(X @ beta_new) + 1e-15))
ll_trace.append(ll)
beta = beta_new
beta_trace.append(beta.copy())
if np.linalg.norm(beta_trace[-1] - beta_trace[-2]) < 1e-8:
break
beta_trace = np.array(beta_trace)
# 반복 경로 시각화
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
axes[0].plot(range(len(ll_trace)), ll_trace, 'bo-', ms=4)
axes[0].set_xlabel('반복 횟수 (t)')
axes[0].set_ylabel('로그 우도 ℓ(β)')
axes[0].set_title('IRLS: 로그 우도 수렴')
axes[0].grid(True, alpha=0.3)
axes[1].plot(beta_trace[:, 0], beta_trace[:, 1], 'bo-', ms=4, lw=1.5)
axes[1].plot(beta_trace[0, 0], beta_trace[0, 1], 'gs', ms=8, label='초기값')
axes[1].plot(beta_trace[-1, 0], beta_trace[-1, 1], 'r*', ms=12, label='MLE')
axes[1].set_xlabel('α 추정값')
axes[1].set_ylabel('β 추정값')
axes[1].set_title('IRLS: 모수 공간에서의 수렴 경로')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# ──────────────────────────────────────
# statsmodels 검증
# ──────────────────────────────────────
model = sm.Logit(failure, X)
result = model.fit(disp=False)
print(result.summary2())
# 이탈도 기반 LRT (H0: β = 0)
X_null = X[:, [0]] # 절편만
model_null = sm.Logit(failure, X_null)
result_null = model_null.fit(disp=False)
lrt_stat = result_null.deviance - result.deviance
df_diff = result_null.df_resid - result.df_resid
p_lrt = chi2.sf(lrt_stat, df=df_diff)
print(f"\nLRT 통계량 = {lrt_stat:.4f}, df = {df_diff}, p = {p_lrt:.4f}")
# 헤시안 (음의 정보행렬) 직접 확인
pi_hat = result.predict()
w_hat = pi_hat * (1 - pi_hat)
H_neg = X.T @ np.diag(w_hat) @ X
print(f"\n-Hessian = Fisher Info:\n{H_neg.round(4)}")
print(f"(-Hessian)^(-1) (분산 추정):\n{np.linalg.inv(H_neg).round(4)}")
print(f"(표준오차)^2: se(α̂)²={np.linalg.inv(H_neg)[0,0]:.4f}, se(β̂)²={np.linalg.inv(H_neg)[1,1]:.4f}")결과 해석: 좌측 패널(로그 우도 수렴)에서 초기 몇 번의 반복 만에 로그 우도가 가파르게 증가한 뒤 평탄해지는 것을 확인할 수 있다 — 순오목성이 보장하는 빠른 수렴의 시각적 증거이다. 우측 패널(모수 공간 수렴 경로)에서 \((\alpha, \beta)\) 가 \((0, 0)\) 초기값에서 MLE \((15.04, -0.232)\) 를 향해 부드럽게 이동한다. statsmodels의 LRT 통계량 \(\approx 7.95\), \(p \approx 0.005\) 는 기온 변수가 유의미하게 O-링 실패율을 설명함을 확인한다. $-H = $ Fisher Info 결과가 Step 1에서 직접 계산한 \(\mathbf{X}^T\mathbf{W}\mathbf{X}\) 와 일치하면 구현이 정확함을 이중 검증한 것이다.
9 요약
| 항목 | 핵심 내용 |
|---|---|
| 헤시안 | \(H = -\mathbf{X}^T\mathbf{W}\mathbf{X}\), \(w_i = \pi_i(1-\pi_i)\) |
| 순오목성 증명 | \(\mathbf{v}^T H\mathbf{v} = -\sum w_i(v_1+v_2 x_i)^2 < 0\) (예측변수 비상수 시) |
| N-R 업데이트 | \(\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} + (\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})^{-1}\mathbf{X}^T(\mathbf{y}-\hat{\boldsymbol{\pi}}^{(t)})\) |
| 정준 링크 특수성 | 관측 정보 = 기대 정보 → N-R = Fisher scoring |
| 조정 반응변수 | \(z_i = \hat{\eta}_i + (y_i-\hat{\pi}_i)/w_i\) — 링크의 1차 테일러 전개 |
| IRLS | \(\hat{\boldsymbol{\beta}}^{(t+1)} = (\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{W}^{(t)}\mathbf{z}^{(t)}\) |
| N-R = IRLS | 동치 — WLS 풀기의 반복이 N-R 업데이트와 같다 |
| 이탈도 | \(D = 2[\ell(\text{포화}) - \ell(\text{적합})]\), LRT와 직접 연결 |
| 초기값 | \(\tilde{\mu}_i = (y_i+0.5)/(m_i+1)\) 권장 (분리 문제 회피) |
| 수렴 진단 | 계수 변화량 + 이탈도 변화량 병행 — 분리 시 계수 발산해도 이탈도 수렴 |
10 관련 주제
선행 지식
- Logistic Regression: The Model — GLM 3요소, 스코어 방정식, 정보행렬
- 최대우도추정량 — MLE 일반론
- M-추정량 — 손실함수 최적화 관점
후속 주제
- Robust Regression — Casella §12.4, M-추정, IRLS의 로버스트 일반화
- Multiple Linear Regression — 행렬 OLS, Gauss-Markov
- GLM 확장 — 포아송 회귀, 감마 회귀, 과산포(quasi-likelihood)
- 다항 로지스틱 회귀 — McCullagh & Nelder Ch.5
관련 개념
- 지수족 분포 — IRLS 일반 유도의 기반
- 점근 가설 검정 — LRT의 χ² 극한 분포
- EIV: MLE 유도 — 다른 비선형 MLE 사례