1 도입 — 왜 수치 최적화인가
Klein 책의 본문 (Ch.1-13) 에서 나온 추정 절차의 거의 전부가 다음 한 가지 문제로 환원된다:
\[ \widehat{\theta} = \arg\max_{\theta} L(\theta) \quad \text{또는} \quad \arg\max_{\theta} \ell(\theta) = \log L(\theta) . \]
| 본문 위치 | 추정 대상 | 우도/score 형태 |
|---|---|---|
| Ch.8 Cox PH | \(\hat{\beta}\) | 부분우도 \(L_p(\beta)\) |
| Ch.10 Aalen 가산 | \(\widehat{B}_k(t)\) | 시점별 OLS (closed-form 있음) |
| Ch.10 Lin-Ying | \(\hat{\alpha}\) | \(\hat{\alpha} = A^{-1} B\) (closed-form) |
| Ch.12 모수 회귀 (Weibull, LL, LN) | \((\hat{\mu}, \hat{\sigma}, \hat{\gamma})\) | full 우도 \(L(\mu, \sigma, \gamma)\) |
| Ch.13 Gamma frailty EM | \((\hat{\theta}, \hat{\beta}, \widehat{H}_0)\) | E + M step |
| Ch.13 Marginal model | \(\hat{\beta}\) | 일반 Cox + sandwich |
이 중 closed-form 이 있는 것은 Aalen, Lin-Ying 두 가지뿐. 나머지 — Cox 부분우도, 모수 회귀 MLE, frailty EM 의 M-step, 일반화 감마, 시간의존 공변량 — 는 모두 반복적 수치 최적화 에 의존한다.
생존분석의 우도는 거의 항상 다음 두 인자를 포함한다:
- 검열 (censoring): 사건 관찰자는 밀도 \(f\), 검열자는 생존함수 \(S\) 로 우도에 들어가 \(L = \prod f^\delta S^{1-\delta}\) 형태. 이 곱 우도는 \(\theta\) 에 대한 비선형 함수.
- 위험집합 합 (risk set sum): Cox 부분우도, frailty marginal 우도 모두 \(\sum_{j \in R(t)} \exp(\beta^t Z_j)\) 형태가 들어가 \(\beta\) 의 비선형 함수.
이 두 비선형성이 합쳐져 score 방정식 \(\partial \ell / \partial \theta = 0\) 의 닫힌 해가 거의 존재하지 않는다. 따라서 score 의 영점을 수치적으로 찾는 것 이 표준 절차이다.
부록 A 는 이 “score 영점 찾기” 의 표준 도구 6 개 (단변수 3 + 다변수 3) 를 제시한다.
본 포스트는 § A.1 단변수 방법과 § A.2 다변수 방법의 알고리즘과 직관, Klein 의 Weibull MLE 두 예제로 본 수렴 비교, 그리고 생존분석 본문에서의 실제 응용을 정리한다.
2 § A.1 Univariate Methods — 단변수 score 영점 찾기
2.1 문제 정의
\(f(x)\) 의 1차원 최대값을 찾는다. 정칙 조건 (\(f\) 미분 가능, 단봉) 하에서:
\[ \hat{x} = \arg\max_x f(x) \;\Leftrightarrow\; f'(\hat{x}) = 0 \;\text{and}\; f''(\hat{x}) < 0 . \]
따라서 score 함수 \(f'(x)\) 의 영점 (root) 을 찾는 문제 로 환원된다. Klein 은 세 가지 방법을 제시한다.
세 방법 모두 \(f'(x) = 0\) 의 해를 찾을 뿐, 그 해가 최대인지 최소인지 변곡점인지 알지 못한다.
수렴 후 반드시 \(f''(\hat{x}) < 0\) (단변수의 음의 정칙) 을 확인하거나, 또는 인접한 두 점에서 \(f\) 값을 비교해 진짜 최대값임을 검증해야 한다. 우도 함수가 다봉형이면 여러 다른 시작값으로 시도해 전역 최대를 보장한다.
2.2 Bisection 이분법 — 가장 견고한 시작점
2.2.1 알고리즘
\(f'(x)\) 의 영점이 구간 \([x_L, x_U]\) 에 있다는 정보가 있다고 가정 (\(f'(x_L) \cdot f'(x_U) < 0\), 즉 양 끝의 score 가 부호가 다름).
While (x_U - x_L) > tol:
x_N = (x_L + x_U) / 2 # 중점
if f'(x_L) * f'(x_N) > 0:
x_L ← x_N # 해가 [x_N, x_U]
else:
x_U ← x_N # 해가 [x_L, x_N]각 step 에서 구간 길이가 정확히 절반으로 줄어든다 → \(k\) step 후 길이는 초기의 \(2^{-k}\). 따라서 \(\log_2(\text{초기 폭}/\text{허용오차})\) step 으로 무조건 수렴.
장점: 초기 구간만 잡으면 반드시 수렴 한다. 함수 모양 무관.
단점: 수렴 속도가 선형 (\(k\) step 에 정확도 \(2^{-k}\)). 다른 방법이 가능하면 더 빠른 것을 선택.
사용 시점: 초기값 정보가 부실하거나 함수가 거칠 때의 안전 시작.
2.3 Secant Method 할선법 — 식 (A.1)
2.3.1 알고리즘
두 초기값 \(x_0, x_1\) 에서 시작 (영점을 사이에 두지 않아도 OK). 각 step:
\[ x_{i+1} = x_i - f'(x_i) \cdot \frac{x_i - x_{i-1}}{f'(x_i) - f'(x_{i-1})} . \tag{식 A.1} \]
\(f'\) 의 영점을 찾는 데 \(f'\) 자체의 두 점 \((x_{i-1}, f'(x_{i-1}))\), \((x_i, f'(x_i))\) 을 잇는 직선 (할선, secant) 의 \(x\) 절편을 다음 후보로 삼는다.
이는 \(f'\) 가 거의 직선 이라는 가정 하에서 합리적이다. 실제 \(f'\) 가 곡선이어도 두 점이 충분히 가까우면 직선 근사가 좋다.
수렴 차수: 황금 비율 \(\phi = (1+\sqrt{5})/2 \approx 1.618\) — bisection (선형, 1) 보다 빠르고 Newton (2) 보다 느림. 그러나 2차 도함수 계산이 불필요 하다는 결정적 장점.
사용 시점: \(f''\) 가 복잡하거나 계산 비용이 큰 경우. Cox 부분우도의 Hessian 처럼 비용이 큰 경우 secant 가 유리.
2.4 Newton-Raphson 방법 — 식 (A.2)
2.4.1 알고리즘
단일 초기값 \(x_0\) 에서 시작:
\[ x_{i+1} = x_i - \frac{f'(x_i)}{f''(x_i)} . \tag{식 A.2} \]
\(x_i\) 근처에서 \(f'\) 의 1차 Taylor 근사: \(f'(x) \approx f'(x_i) + f''(x_i)(x - x_i)\). 이 근사를 0 으로 놓고 \(x\) 에 대해 풀면:
\[ x = x_i - \frac{f'(x_i)}{f''(x_i)} . \]
이게 식 (A.2). 즉 \(f'\) 의 접선 의 \(x\) 절편을 다음 추정값으로.
수렴 차수: 2 차 (quadratic) — 매 step 에서 정확도가 거의 제곱된다. 정확도 \(10^{-2} \to 10^{-4} \to 10^{-8}\) 의 폭발적 가속.
장점: 초기값이 좋으면 매우 빠름. Cox/MLE 의 표준 알고리즘.
단점: - 초기값이 나쁘면 발산 또는 잘못된 영점으로 수렴 가능. - \(f''\) 계산 비용이 큰 경우 (예: 다변수 Hessian) 부담. - \(f''(x_i)\) 가 0 에 가까우면 step 이 폭발 (수치 불안정).
사용 시점: 좋은 초기값 + \(f''\) 계산 가능 + 빠른 수렴 원할 때.
2.5 세 방법 비교 — 식 (A.1), (A.2)
| 방법 | 정보 요구 | 초기값 | 수렴 차수 | 견고성 |
|---|---|---|---|---|
| Bisection | \(f'\) 만 | 부호 다른 두 점 | 1 (선형) | 매우 강 |
| Secant | \(f'\) 만 | 두 점 (자유) | \(\phi \approx 1.618\) | 중간 |
| Newton-Raphson | \(f', f''\) | 한 점 | 2 (이차) | 약 (초기값 의존) |
Klein 도 명시하지 않은 결합 패턴이 표준 실무이다:
1 단계 — Bisection 으로 초기 구간 좁히기: 거친 단계에서 견고한 bisection 으로 영점 부근 찾기 (예: 2-3 step 만).
2 단계 — Newton-Raphson 으로 마무리: 영점 근처에서 빠른 2차 수렴.
이 결합이 R 의 uniroot(), Python scipy.optimize.brentq 등의 내부 알고리즘이다 (Brent’s method = bisection + secant + inverse quadratic 의 결합).
따라서 직접 구현 시는 보통 Newton-Raphson 만 써도 충분 (좋은 초기값 + safeguard) 하지만, 라이브러리는 결합 알고리즘을 사용해 양쪽 장점을 취한다.
2.6 Example A.1 — Weibull 단모수 MLE
2.6.1 데이터
10 개의 비검열 관측 (단순화된 생존시간):
2.57, 0.58, 0.82, 1.02, 0.78, 0.46, 1.04, 0.43, 0.69, 1.37진짜 모형: \(\lambda = 1\) 고정 의 Weibull, \(h(t) = \alpha t^{\alpha-1} e^{-t^\alpha}\). 형상 모수 \(\alpha\) 만 추정.
2.6.2 우도와 score
\[ \ell(\alpha) = n \ln \alpha + (\alpha - 1) \sum \ln t_j - \sum t_j^\alpha . \]
\[ f'(\alpha) = \frac{n}{\alpha} + \sum \ln t_j - \sum t_j^\alpha \ln t_j . \]
\[ f''(\alpha) = -\frac{n}{\alpha^2} - \sum t_j^\alpha [\ln t_j]^2 . \]
2.6.3 세 방법 적용 — Klein 본문 결과
| 방법 | 시작값 | 수렴 step | \(\hat{\alpha}\) |
|---|---|---|---|
| Bisection | \([\alpha_L, \alpha_U] = [1.5, 2]\) | 8 | 1.705 |
| Secant | \(\alpha_0 = 1, \alpha_1 = 1.5\) | 2 | 1.705 |
| Newton-Raphson | \(\alpha_0 = 1.5\) | 2 | 1.705 |
(정지 기준: \(|f'(\alpha)| < 0.01\).)
세 방법 모두 같은 답 (\(\hat{\alpha} = 1.705\)) 에 수렴 — 단일 단봉 우도라 어떤 방법이든 OK.
수렴 속도 비교:
- Bisection: 8 step (선형 수렴, 매 step 절반)
- Secant: 2 step (sub-quadratic, 두 step 만에 거의 정확)
- Newton-Raphson: 2 step (quadratic, 같은 2 step 이지만 더 정확한 끝값)
Newton-Raphson 의 첫 step: \(\alpha_0 = 1.5 \to \alpha_1 = 1.701\). Secant 의 첫 step (\(\alpha_0 = 1, \alpha_1 = 1.5 \to \alpha_2 = 1.671\)) 보다 정답에 더 가까움. 2차 수렴의 위력.
이 예제는 모수 1 개라 Newton-Raphson 의 약점 (Hessian 비용) 이 거의 무. 실제 생존분석 모형은 모수 \(p \geq 5\) 가 흔하므로 Hessian 의 \(p \times p\) 행렬 역수가 부담이 된다 — 그래서 § A.2 다변수 방법이 등장.
3 § A.2 Multivariate Methods — 다변수 score 영점 찾기
3.1 표기
\(f(\mathbf{x})\) 의 \(p\)-차원 최대값을 찾는다. 다음 두 양이 핵심.
Gradient (1차 도함수 벡터):
\[ \mathbf{u}(\mathbf{x}) = [u_1(\mathbf{x}), \ldots, u_p(\mathbf{x})]^t , \quad u_j(\mathbf{x}) = \frac{\partial f(\mathbf{x})}{\partial x_j} . \tag{식 A.3} \]
Hessian (2차 도함수 행렬):
\[ \mathbf{H}(\mathbf{x}) = (H_{ij}(\mathbf{x})) , \quad H_{ij}(\mathbf{x}) = \frac{\partial^2 f(\mathbf{x})}{\partial x_i \partial x_j} . \tag{식 A.4} \]
\(\mathbf{u}(\mathbf{x})\) 는 점 \(\mathbf{x}\) 에서 \(f\) 가 가장 빠르게 증가하는 방향을 가리키는 단위 벡터에 그 비율을 곱한 것이다 (단, normalize 안 한 형태). 영점 \(\mathbf{u}(\mathbf{x}^*) = \mathbf{0}\) 인 점 \(\mathbf{x}^*\) 가 임계점.
Hessian \(\mathbf{H}\) 는 그 점의 곡률 (curvature). \(\mathbf{H}\) 가 음의 정칙 (negative definite) 이면 그 임계점이 최대.
다변수 score 영점 찾기는 단변수와 비슷하지만 방향 + 보폭 을 둘 다 결정해야 한다. 세 방법은 이 결정 방식이 다르다.
3.2 Steepest Ascent — 식 (A.5)
3.2.1 알고리즘
각 step 에서 그래디언트 방향으로 직선 탐색:
\[ \mathbf{x}_{k+1} = \mathbf{x}_k + d_k \, \mathbf{u}(\mathbf{x}_k) , \tag{식 A.5} \]
여기서 \(d_k\) 는 그 방향에서 \(f\) 를 최대화하는 step size (1D 최적화 — bisection/Newton 등으로 1 차원 풀이).
Hessian 정보를 안 쓴다 — gradient 만 보고 그 방향으로 1D 최적화. 따라서:
장점: Hessian 불필요. 초기값 견고 (수렴 보장).
단점: 수렴 매우 느림. 길쭉한 우도 함수 (한 축의 곡률이 다른 축보다 크게 다름) 에서 지그재그로 진동하며 천천히 수렴.
사용 시점: 초기값이 매우 나쁠 때의 첫 몇 step 만. 또는 Hessian 계산이 불가능할 때 (예: 비미분 함수의 근사).
3.3 Newton-Raphson 다변수 — 식 (A.6)
3.3.1 알고리즘
\[ \mathbf{x}_{k+1} = \mathbf{x}_k - \mathbf{H}(\mathbf{x}_k)^{-1} \mathbf{u}(\mathbf{x}_k) . \tag{식 A.6} \]
단변수 식 (A.2) \(x_{i+1} = x_i - f'(x_i)/f''(x_i)\) 에서 \(f'/f'' \to \mathbf{H}^{-1} \mathbf{u}\) 로 대체. Hessian 의 역수가 곡률 보정 역할.
또한 식 (A.6) 의 \(-\mathbf{H}^{-1} \mathbf{u}\) 가 Hessian-orthogonal 좌표계에서의 steepest ascent 와 같다. 즉 Hessian 으로 좌표를 정규화한 후 steepest ascent 를 적용한 것과 동치. 이는 길쭉한 우도 함수의 지그재그 문제를 자동 해결.
수렴 차수: 2 차 (quadratic). 영점 부근에서 매우 빠름.
단점 (단변수와 같음): - 초기값이 나쁘면 발산 또는 잘못된 방향. - \(\mathbf{H}\) 가 singular 또는 ill-conditioned 면 수치 불안정. - \(p \times p\) Hessian 행렬 역수 비용 (\(O(p^3)\)) — \(p\) 클 때 부담.
대처: \(f(\mathbf{x}_k) > f(\mathbf{x}_{k+1})\) (잘못된 방향) 이 검출되면 step 크기를 절반으로 줄임 (\(\mathbf{x}_{k+1} = \mathbf{x}_k - \mathbf{H}^{-1}\mathbf{u}/2\)) — SAS 와 BMDP 의 Cox 회귀가 이 방식 사용.
3.4 Marquardt’s Method — Steepest Ascent 와 Newton-Raphson 의 절충
3.4.1 알고리즘
\(\mathbf{S}_k\) 를 대각 스케일 행렬 (\(S_{ii} = |H_{ii}(\mathbf{x}_k)|^{-1/2}\)). \(\gamma \geq 0\) 을 blending 상수.
\[ \mathbf{x}_{k+1} = \mathbf{x}_k - \mathbf{S}_k \big(\mathbf{S}_k \mathbf{H}(\mathbf{x}_k) \mathbf{S}_k + \gamma \mathbf{I}\big)^{-1} \mathbf{S}_k \mathbf{u}(\mathbf{x}_k) . \]
- \(\gamma = 0\): \((\mathbf{S}_k \mathbf{H} \mathbf{S}_k)^{-1}\) 만 남아 정규화된 Newton-Raphson.
- \(\gamma \to \infty\): 분모의 \(\gamma \mathbf{I}\) 가 지배 → \(\mathbf{S}_k \cdot \frac{\mathbf{S}_k \mathbf{u}}{\gamma}\) → steepest ascent (Hessian 영향 사라짐).
적응 절차:
- 작은 \(\gamma\) 로 시작 (예: \(\gamma = 0.01\)).
- 한 step 후 \(f(\mathbf{x}_{k+1}) > f(\mathbf{x}_k)\) 면 OK, 다음 step.
- \(f(\mathbf{x}_{k+1}) \leq f(\mathbf{x}_k)\) 면 (방향 잘못됨) \(\gamma\) 를 10 배 늘려 재시도.
- 수렴 직전 마지막 step 은 \(\gamma = 0\) (순수 Newton-Raphson) 으로 정확 마무리.
Marquardt (1963) 의 핵심 통찰: 초기에는 Hessian 정보가 부정확할 수 있으므로 (\(\mathbf{H}\) 가 양의 정칙일 수도) steepest ascent 쪽에 가깝게, 영점 근처에서는 Newton-Raphson 의 빠른 수렴 활용.
비유: 처음에는 안전 운전 (steepest ascent), 도착지 가까워지면 액셀 (Newton-Raphson). 비선형 회귀 (nls() in R) 의 default 알고리즘.
3.5 수렴 판정 — 4 가지 기준
다변수 최적화의 수렴 여부 결정은 다음 중 하나 (또는 조합) 로:
- 함수값 변화: \(f(\mathbf{x}_{k+1}) - f(\mathbf{x}_k) < \epsilon\) 또는 상대 변화 \(|[\,]/f(\mathbf{x}_k)| < \epsilon\).
- 그래디언트 노름: \(\sum_j u_j(\mathbf{x}_{k+1})^2 < \epsilon\) 또는 \(\max_j |u_j(\mathbf{x}_{k+1})| < \epsilon\).
- 모수 변화: \(\sum_j (x_{k+1, j} - x_{k, j})^2 < \epsilon\) 또는 \(\max_j |x_{k+1, j} - x_{k, j}| < \epsilon\).
- 반복 한계: \(k \geq k_{\max}\) — 발산 방지 안전장치.
그래디언트 기준: 이론적으로 가장 정직 (\(\mathbf{u} = 0\) 이 최적성 조건). 그러나 \(\mathbf{u}\) 의 척도가 모수 차원마다 달라 단일 \(\epsilon\) 으로 설정하기 까다로움.
모수 변화 기준: 직관적이지만 평탄한 우도에서 가짜 수렴 가능 (step 이 작아져도 진짜 최적이 아닐 수 있음).
함수값 변화 기준: 가장 흔히 사용. 우도가 거의 안 변한다 = 사실상 수렴.
실무: 함수값 + 그래디언트 둘 다 점검. 상대 변화 \(\epsilon = 10^{-6}\) 또는 \(10^{-8}\) 권장.
3.6 Example A.2 — Weibull 두 모수 MLE
같은 10 개 데이터로 \(S(t) = \exp(-\lambda t^\alpha)\) 의 두 모수 추정. 시작값 \((\lambda_0, \alpha_0) = (1.024, 1.000)\) (\(\lambda_0 = 10/\sum t_i\) 의 method-of-moments 추정).
3.6.1 Score 와 Hessian
\[ \begin{aligned} u_\lambda &= \frac{n}{\lambda} - \sum t_i^\alpha , \\ u_\alpha &= \frac{n}{\alpha} + \sum \ln t_i - \lambda \sum t_i^\alpha \ln t_i , \end{aligned} \]
\[ \mathbf{H}(\lambda, \alpha) = \begin{pmatrix} -\frac{n}{\lambda^2} & -\sum t_i^\alpha \ln t_i \\ -\sum t_i^\alpha \ln t_i & -\frac{n}{\alpha^2} - \lambda \sum t_i^\alpha (\ln t_i)^2 \end{pmatrix} . \]
3.6.2 세 방법 비교 — Klein 본문 결과
| 방법 | 수렴 step | \((\hat{\lambda}, \hat{\alpha})\) |
|---|---|---|
| Steepest ascent | 5 | \((0.839, 1.792)\) |
| Newton-Raphson | 3 | \((0.832, 1.796)\) |
| Marquardt (\(\gamma = 0.5\)) | (~3-4) | 수렴 |
Newton-Raphson 의 우월성: 단 3 step 으로 수렴 — Hessian 정보를 활용해 직접 영점으로 점프.
Steepest ascent 의 지그재그: 5 step 필요. 자세히 보면:
Step 0: λ=1.024, α=1.000 → α 방향으로 큰 step (u_α=7.035)
Step 1: λ=1.025, α=1.693 → λ 방향으로 step (u_λ=0.661)
Step 2: λ=0.865, α=1.694 → α 방향 보정
Step 3: λ=0.865, α=1.777
Step 4: λ=0.839, α=1.777
Step 5: λ=0.839, α=1.792 (수렴)각 step 이 한 모수씩 번갈아 큰 변화 — 우도 표면이 길쭉할 때 steepest ascent 의 전형적 비효율 패턴.
Marquardt 의 절충: 초기 step 은 steepest 쪽 (안전), 후반은 Newton 쪽 (빠름). 실무 default.
Newton-Raphson 의 미세 차이: 같은 데이터에 steepest 가 \((0.839, 1.792)\) Newton 이 \((0.832, 1.796)\) 으로 약간 다름 — 수렴 정밀도 차이 (모두 \(|u| < 0.1\) 까지만 풀어 정확한 최적값에 충분히 가까이는 안 옴). \(\epsilon\) 을 더 작게 하면 일치한다.
4 생존분석에서의 응용 — 수치 최적화가 어디에 등장하는가
문제: \(\partial \log L_p(\beta) / \partial \beta = 0\). \(\beta \in \mathbb{R}^p\) 이라 다변수.
Score: \(u_k(\beta) = \sum_j \delta_j [Z_{jk} - \bar{Z}_k(t_j; \beta)]\). (\(\bar{Z}\) 는 위험집합 가중 평균.)
Hessian: \(H_{kl}(\beta) = -\sum_j \delta_j \text{Var}_{R(t_j)}(Z_k, Z_l; \beta)\). (위험집합의 가중 공분산.)
알고리즘: Newton-Raphson (식 A.6). R coxph() 의 표준. 수렴 보통 4-7 iteration.
문제: \(\arg\max_{(\mu, \sigma, \gamma)} L(\mu, \sigma, \gamma; \text{data})\) — Weibull, log-logistic, log-normal 등.
모수 수: \(1 + 1 + p\) (intercept + scale + 회귀 계수).
알고리즘: Newton-Raphson 또는 quasi-Newton. R survreg() 의 표준. 시작값을 OLS 등으로 보강하면 수렴 안정.
문제: M-step 에서 \(\hat{u}_i\) 를 known offset 으로 둔 Cox partial likelihood 최대화 (식 13.3.7) + \(\theta\) 의 1D 최대화.
알고리즘: - \(\beta\) 업데이트: 다변수 Newton-Raphson (Cox 와 동일). - \(\theta\) 업데이트: 단변수 Newton-Raphson 또는 grid search (Nielsen profile EM).
EM 전체는 외부 반복 + 내부 (M-step) 반복의 이중 구조. 외부 수렴 느림 → Nielsen profile EM 으로 가속.
일반화 감마 (Ch.12.4): 3 모수 (\(\mu, \sigma, \theta\)), 어떤 분포든 포괄. Newton-Raphson 으로 풀이 — Marquardt 권장 (초기 hessian 약함).
시간의존 공변량 (Ch.9): \(\beta\) 가 정상이지만 \(Z(t)\) 가 시간 변동. Cox partial likelihood 의 표현 변경, 알고리즘 (Newton-Raphson) 동일.
5 코드 예시
5.1 Step 1 — Python: scipy.optimize 의 통합 API
import numpy as np
from scipy.optimize import minimize_scalar, minimize, newton
# 1. Univariate (Example A.1) — Weibull α MLE
data = np.array([2.57, 0.58, 0.82, 1.02, 0.78, 0.46, 1.04, 0.43, 0.69, 1.37])
n = len(data)
def neg_log_lik(alpha):
return -(n * np.log(alpha) + (alpha - 1) * np.sum(np.log(data))
- np.sum(data ** alpha))
# Brent's method (bisection + secant + inverse quadratic)
res = minimize_scalar(neg_log_lik, bracket=(0.5, 3.0), method='brent')
print(f"Brent: α̂ = {res.x:.4f}, iters = {res.nit}")
# Newton-Raphson (직접 구현 — 식 A.2)
def score(alpha):
return n / alpha + np.sum(np.log(data)) - np.sum(data ** alpha * np.log(data))
def hessian(alpha):
return -n / alpha**2 - np.sum(data ** alpha * np.log(data) ** 2)
alpha_hat = newton(score, x0=1.5, fprime=hessian, tol=1e-6)
print(f"Newton: α̂ = {alpha_hat:.4f}")
# 2. Multivariate (Example A.2) — Weibull (λ, α) MLE
def neg_log_lik_2(params):
lam, alpha = params
if lam <= 0 or alpha <= 0:
return np.inf
return -(n * np.log(lam) + n * np.log(alpha) + (alpha - 1) * np.sum(np.log(data))
- lam * np.sum(data ** alpha))
# BFGS (quasi-Newton, default for unconstrained smooth)
res = minimize(neg_log_lik_2, x0=[1.0, 1.0], method='BFGS')
print(f"BFGS: λ̂ = {res.x[0]:.4f}, α̂ = {res.x[1]:.4f}, iters = {res.nit}")
# Levenberg-Marquardt 변형 (least_squares 통해)
# scipy.optimize.minimize 는 직접 LM 안 지원, least_squares 의 method='lm' 으로5.2 Step 2 — R: optim() 통합 인터페이스
# 데이터
data <- c(2.57, 0.58, 0.82, 1.02, 0.78, 0.46, 1.04, 0.43, 0.69, 1.37)
n <- length(data)
# 음의 log-likelihood (R optim 은 minimize)
neg_log_lik <- function(params) {
lam <- params[1]
alpha <- params[2]
if (lam <= 0 || alpha <= 0) return(Inf)
-(n * log(lam) + n * log(alpha) + (alpha - 1) * sum(log(data))
- lam * sum(data ^ alpha))
}
# BFGS (default for smooth)
fit_bfgs <- optim(par = c(1.0, 1.0), fn = neg_log_lik, method = "BFGS")
cat("BFGS: λ̂ =", fit_bfgs$par[1], "α̂ =", fit_bfgs$par[2], "\n")
# Nelder-Mead (gradient-free, robust for noisy/non-smooth)
fit_nm <- optim(par = c(1.0, 1.0), fn = neg_log_lik, method = "Nelder-Mead")
cat("NM: λ̂ =", fit_nm$par[1], "α̂ =", fit_nm$par[2], "\n")
# L-BFGS-B (with box constraints, lambda > 0, alpha > 0)
fit_lbfgs <- optim(par = c(1.0, 1.0), fn = neg_log_lik, method = "L-BFGS-B",
lower = c(0.01, 0.01))
cat("LBFGS:λ̂ =", fit_lbfgs$par[1], "α̂ =", fit_lbfgs$par[2], "\n")
# Levenberg-Marquardt — minpack.lm 패키지
library(minpack.lm)
# nlsLM() 은 nls() 의 LM 변형 — 회귀 형태에 사용5.3 Step 3 — survreg / coxph 내부 알고리즘 확인
library(survival)
# Cox PH — coxph 는 Newton-Raphson 사용 (with step-halving)
fit_cox <- coxph(Surv(time, status) ~ x, data = ...)
fit_cox$iter # 수렴 step 수
fit_cox$loglik # 초기 + 수렴 시 log-likelihood
# Weibull AFT — survreg 는 Newton-Raphson
fit_wei <- survreg(Surv(time, status) ~ x, dist = "weibull", data = ...)
fit_wei$iter
fit_wei$loglik6 핵심 요약
생존분석의 거의 모든 추정 (Cox · 모수 회귀 · frailty EM) 은 score 방정식 \(\partial \ell / \partial \theta = 0\) 의 수치 영점 찾기로 환원된다. 단변수에서는 bisection (견고·느림) → secant (식 A.1, 빠름·\(f''\) 불필요) → Newton-Raphson (식 A.2, 가장 빠름) 의 trade-off. 다변수에서는 steepest ascent (식 A.5, 견고·느림) → Newton-Raphson (식 A.6, 빠름·Hessian 비용) → Marquardt 절충 (\(\gamma\) 다이얼로 두 방법 결합) 의 trade-off. 실무에서는 R
optim의 BFGS, Pythonscipy.optimize의 Brent / BFGS / LM 등이 이 알고리즘들의 결합 형태로 자동 사용된다.
| 절 | 방법 | 핵심 식 | 수렴 차수 | 견고성 |
|---|---|---|---|---|
| A.1 | Bisection | 절반씩 좁힘 | 1 (선형) | 매우 강 |
| A.1 | Secant | 식 (A.1) | \(\phi \approx 1.618\) | 중간 |
| A.1 | Newton-Raphson | 식 (A.2) | 2 (이차) | 약 |
| A.2 | Steepest Ascent | 식 (A.5) | 선형 | 매우 강 |
| A.2 | Newton-Raphson 다변수 | 식 (A.6) | 2 (이차) | 약 |
| A.2 | Marquardt | \(\gamma\) blending | 적응 | 중간-강 |
1. 좋은 시작값이 가장 중요: Newton-Raphson 의 빠른 수렴은 시작값이 영점 근처일 때만. Method-of-moments, OLS 등으로 사전 추정 후 시작.
2. Hessian 비용이 크면 Quasi-Newton (BFGS): \(\mathbf{H}\) 를 score 변화로 근사. R optim, Python scipy.optimize.minimize 의 default. Cox 보다 더 큰 모형에서 표준.
3. 비선형 회귀에는 Marquardt (LM): 초기에 안전, 수렴 부근에 빠름. R nls(), minpack.lm::nlsLM(), Python scipy.optimize.least_squares(method='lm').
4. Step-halving 으로 발산 방지: \(f(\mathbf{x}_{k+1}) < f(\mathbf{x}_k)\) 면 step 절반으로. SAS Cox 의 표준.
5. 수렴 후 검증: \(\mathbf{H}(\hat{\mathbf{x}})\) 가 음의 정칙인지 (eigenvalue < 0) 확인. 안 그러면 saddle point 또는 local min — 다른 시작값으로 재시도.
7 관련 주제
Appendix A 시리즈
- § A.1 — Univariate Methods — Bisection·Secant·Newton-Raphson 의 알고리즘·수렴 차수 유도·발산 사례·Example A.1 Weibull α MLE step trace
- § A.2 — Multivariate Methods — Steepest Ascent·Newton-Raphson 다변수·Marquardt 절충, 길쭉한 우도의 지그재그·Hessian 정규화 직관, Example A.2 Weibull (λ, α) MLE step trace, BFGS quasi-Newton
Klein 본문 시리즈 (수치 최적화 응용 위치)
- Ch.8 — Cox 비례위험 모형 — 부분우도 Newton-Raphson
- Ch.12 — 모수적 회귀 모형 — Weibull/LL/LN MLE
- § 13.3-13.4 — Gamma Frailty EM & Marginal — EM + Newton-Raphson 이중 구조
관련 부록
- Appendix B — Large-Sample Tests Based on Likelihood Theory — LR · Wald · Score 세 가능도 검정의 식·직관·점근 동등성 (수치 최적화로 구한 MLE 의 가설 검정 도구)
관련 개념
8 참고 문헌
- Klein, J. P., & Moeschberger, M. L. (2003). Survival Analysis: Techniques for Censored and Truncated Data (2nd ed.). Springer. Appendix A.
- Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters. J. SIAM, 11(2), 431-441. (Marquardt 의 원논문)
- Thisted, R. A. (1988). Elements of Statistical Computing: Numerical Computation. Chapman and Hall. (Klein 이 인용한 통계 수치 계산 표준 교재)
- Nocedal, J., & Wright, S. J. (2006). Numerical Optimization (2nd ed.). Springer. (현대 표준 교재 — BFGS, L-BFGS, trust region 등 포괄)
- Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (2007). Numerical Recipes (3rd ed.). Cambridge. (Brent’s method, Marquardt 의 실용 구현)