최대가능도 최적화 보기¶

(1) 해석적으로. \(\sigma = \hat\sigma\)를 고정하면

인데, 제곱합을 \(\bar x\) 둘레로 펴면

이다(교차항 \(2(\bar x - \mu)\sum_i(x_i - \bar x)\)가 \(0\)이다). 따라서

로 꼭짓점이 \(\bar x\)인 정확한 포물선이다. 근사가 아니라 항등식이고, 그러므로

이다.

격자 오차의 상한. 격자 간격이 \(h = (7-3)/499 = 0.0080160\)이므로 \(\bar x\)에서 가장 가까운 격자점은 많아야 반 칸, 곧 \(h/2 = 0.0040080\) 떨어져 있다. 그때 잃는 로그가능도는

이하다. 격자 탐색은 반 칸보다 정밀할 수 없다는 것이 요점이고, 이는 자료를 아무리 늘려도 나아지지 않는 종류의 오차다.

(2) 수치적으로.

import numpy as np
from scipy import stats

def grid_search_normal_mean(data, mu_grid):
    """mu 후보를 격자로 늘어놓고 로그가능도가 가장 큰 것을 고른다.

    가장 단순하고 가장 확실한 방법이다. 미분도 초기값도 필요 없고
    국소 최적에 갇히지도 않는다. 대신 격자 간격보다 정밀할 수 없고,
    모수가 d개면 격자점이 (격자 수)^d 로 폭발해 d가 3~4만 넘어도 못 쓴다.
    """
    sigma_hat = data.std(ddof=0)     # ddof=0 이 MLE 판본이다(n으로 나눔)

    # logpdf 를 더한다. pdf를 곱한 뒤 로그를 취하면 언더플로가 나므로
    # 처음부터 로그로 계산해 더하는 것이 정석이다.
    log_liks = np.array([
        np.sum(stats.norm.logpdf(data, loc=mu, scale=sigma_hat))
        for mu in mu_grid
    ])
    best_idx = np.argmax(log_liks)
    return mu_grid[best_idx], log_liks

# 보기.
rng = np.random.default_rng(42)
data = rng.normal(loc=5.0, scale=2.0, size=100)
mu_grid = np.linspace(3.0, 7.0, 500)
mu_hat, log_liks = grid_search_normal_mean(data, mu_grid)
print(f"Grid search MLE: mu_hat = {mu_hat:.4f}")
print(f"Closed-form MLE: mu_hat = {data.mean():.4f}")

# (1) 에서 따진 상한과 견준다.
step = mu_grid[1] - mu_grid[0]
mu_exact = data.mean()
sigma_hat = data.std(ddof=0)
ll = lambda m: np.sum(stats.norm.logpdf(data, loc=m, scale=sigma_hat))

print(f"\n격자 간격 = {step:.7f},  반 칸 = {step / 2:.7f}")
print(f"격자 - 닫힌해 = {mu_hat - mu_exact:+.7f}  "
      f"(반 칸의 {abs(mu_hat - mu_exact) / (step / 2):.1%})")
print(f"l 의 손실   = {ll(mu_exact) - ll(mu_hat):.6e}")
print(f"이차식 예측 = {0.5 * len(data) / sigma_hat**2 * (mu_hat - mu_exact)**2:.6e}")
print(f"손실의 상한 = {0.5 * len(data) / sigma_hat**2 * (step / 2)**2:.6e}")
print(f"소수 넷째 자리까지 맞히려면 격자점이 {int((7 - 3) / 2e-4) + 1}개 필요하다")

출력:

Grid search MLE: mu_hat = 4.8998
Closed-form MLE: mu_hat = 4.8995

격자 간격 = 0.0080160,  반 칸 = 0.0040080
격자 - 닫힌해 = +0.0003388  (반 칸의 8.5%)
l 의 손실   = 2.402909e-06
이차식 예측 = 2.402909e-06
손실의 상한 = 3.362425e-04
소수 넷째 자리까지 맞히려면 격자점이 20001개 필요하다

(1)의 주장이 모두 확인된다. 격자와 닫힌 해의 차이가 \(0.00034\)로 반 칸 \(0.0040\)의 \(8.5\%\) 안에 들어 상한을 지켰다. 이번에는 운이 좋아 반 칸보다 훨씬 가까웠지만, \(\bar x\)가 어디 떨어지느냐에 따라 최악에는 반 칸까지 벌어질 수 있다.

더 눈여겨볼 것은 로그가능도의 손실이다. 실제 \(2.402909\times10^{-6}\)과 이차식이 예측한 \(2.402909\times10^{-6}\)이 일곱 자리까지 같다. (1)에서 포물선이 근사가 아니라 항등식이라고 한 그대로이며, 어림이 아니라 정확한 예측이다.

격자의 대가는 여기서 분명해진다. \(\hat\mu\)를 소수 넷째 자리까지 맞히려면 간격을 \(2\times10^{-4}\) 아래로 줄여야 하므로 격자점이 \(20{,}001\)개 필요하다. 지금의 \(500\)점보다 \(40\)배다. 닫힌 해 \(\bar x\)는 덧셈 \(100\)번이면 끝나고 자리 손실도 없다. 닫힌 해가 있으면 격자를 쓸 이유가 없고, 격자는 닫힌 해가 없을 때 또는 곡면의 모양을 눈으로 보고 싶을 때 쓰는 도구다.

모수가 \(d\)개면 격자점이 \((\text{격자 수})^d\)로 폭발한다는 점도 함께 기억해 두자. 지금 쓴 \(500\)점을 모수 넷에 똑같이 적용하면 \(500^4 = 6.25\times10^{10}\)점이 되어 계산이 불가능하다.

(1) 해석적으로. \(g(\sigma^2) = \log\sigma^2\)이 \((0,\infty)\)에서 \(\mathbb{R}\)로 가는 전단사이므로

이고 두 최대점이 \(\hat\phi = \log\hat\sigma^2\)로 맞물린다. 최댓값을 옮겨 놓았을 뿐 바꾸지 않았다. 이것이 최대가능도추정량의 불변성이며, 덤으로 \(e^\phi > 0\)이 저절로 보장되어 제약이 사라진다. 최적화기는 \(\phi\)를 실수 전체에서 자유롭게 움직일 수 있다.

닫힌 해는 잘 알려진 대로

이다. 분모가 \(n-1\)이 아니라 \(n\)이라는 점이 중요하다.

(2) 증거가 아니다. 코드의 출발점 목록을 다시 보라.

for mu0 in [0, data.mean(), np.median(data)]:
    for ls0 in [0, np.log(data.var())]:

np.var의 기본값은 ddof=0이므로 data.var()는 정확히 \(\hat\sigma^2\)이다. 곧 출발점 \((\bar x, \log\hat\sigma^2)\)은 이미 정답 그 자리다. 넬더-미드는 그 자리에서 단체를 만들고 어느 꼭짓점도 더 낫지 않으면 \(x_0\)를 그대로 돌려주므로, 일치는 "소수 넷째 자리"가 아니라 비트까지다. 알고리즘의 정확도와는 아무 상관이 없다.

넬더-미드의 실제 정확도를 보려면 나머지 다섯 출발점을 따로 보아야 한다.

import numpy as np
from scipy import optimize

def mle_normal_numerical(data):
    """재매개변수화를 써서 정규분포의 MLE를 수치적으로 찾는다."""
    def neg_log_lik(params):
        # sigma^2 을 직접 다루지 않고 log(sigma^2) 을 최적화한다.
        # 이유: sigma^2 > 0 이라는 제약을 최적화기에 알려 주기 어려운데,
        #       log를 쓰면 log_sigma2 가 어떤 실수든 exp를 거치며 자동으로 양수가 된다.
        #       제약 없는 최적화 문제로 바뀌므로 Nelder-Mead 같은 단순한 방법도 쓸 수 있다.
        mu, log_sigma2 = params
        sigma2 = np.exp(log_sigma2)
        n = len(data)
        # 음의 로그가능도. 최소화하는 것이 가능도를 최대화하는 것과 같다.
        return 0.5 * n * np.log(2 * np.pi * sigma2) + np.sum((data - mu) ** 2) / (2 * sigma2)

    # 출발점을 여러 개 시도한다.
    # 정규분포의 로그가능도는 볼록해서 사실 한 번이면 충분하지만,
    # 봉우리가 여럿인 문제에서는 이렇게 여러 곳에서 출발해
    # 가장 좋은 것을 골라야 국소 최적에 갇히지 않는다.
    # (numpy 배열에는 .median() 이 없으므로 np.median 을 쓴다.)
    best_result = None
    for mu0 in [0, data.mean(), np.median(data)]:
        for ls0 in [0, np.log(data.var())]:
            result = optimize.minimize(neg_log_lik, x0=[mu0, ls0], method="Nelder-Mead")
            # result.fun 이 그 출발점에서 도달한 최솟값이다
            if best_result is None or result.fun < best_result.fun:
                best_result = result

    mu_hat = best_result.x[0]
    sigma2_hat = np.exp(best_result.x[1])    # log에서 되돌린다
    return mu_hat, sigma2_hat

rng = np.random.default_rng(42)
data = rng.normal(5.0, 2.0, 100)
mu_hat, sigma2_hat = mle_normal_numerical(data)
print(f"Numerical MLE: mu = {mu_hat:.4f}, sigma^2 = {sigma2_hat:.4f}")
print(f"Closed-form:   mu = {data.mean():.4f}, sigma^2 = {np.mean((data - data.mean())**2):.4f}")

# 소수 넷째 자리가 아니라 비트까지 같은지 본다.
mu_closed = data.mean()
sigma2_closed = np.mean((data - mu_closed) ** 2)
print(f"\nmu  비트까지 같은가? {mu_hat == mu_closed}")
print(f"s^2 비트까지 같은가? {sigma2_hat == sigma2_closed}")
print(f"np.var(data) 가 MLE 와 같은가? {data.var() == sigma2_closed}")

# 출발점마다 따로 보면 넬더-미드의 실제 정확도가 드러난다.
def neg_ll(params):
    mu, ls = params
    s2 = np.exp(ls)
    return 0.5 * len(data) * np.log(2 * np.pi * s2) + np.sum((data - mu) ** 2) / (2 * s2)

print("\n출발점별 오차")
for mu0, nm in [(0.0, "0"), (mu_closed, "mean"), (float(np.median(data)), "median")]:
    for ls0, ln in [(0.0, "0"), (float(np.log(data.var())), "log var")]:
        r = optimize.minimize(neg_ll, x0=[mu0, ls0], method="Nelder-Mead")
        print(f"  ({nm:>6}, {ln:>7}):  mu 오차 {r.x[0] - mu_closed:>9.2e},"
              f"  s^2 오차 {np.exp(r.x[1]) - sigma2_closed:>9.2e},  호출 {r.nfev:>3}")

출력:

Numerical MLE: mu = 4.8995, sigma^2 = 2.3888
Closed-form:   mu = 4.8995, sigma^2 = 2.3888

mu  비트까지 같은가? True
s^2 비트까지 같은가? True
np.var(data) 가 MLE 와 같은가? True

출발점별 오차
  (     0,       0):  mu 오차  2.17e-06,  s^2 오차  7.58e-05,  호출 152
  (     0, log var):  mu 오차 -2.75e-05,  s^2 오차  8.64e-05,  호출 153
  (  mean,       0):  mu 오차  2.15e-05,  s^2 오차  5.11e-05,  호출 105
  (  mean, log var):  mu 오차  0.00e+00,  s^2 오차  0.00e+00,  호출  45
  (median,       0):  mu 오차  1.16e-05,  s^2 오차  5.96e-05,  호출 102
  (median, log var):  mu 오차  1.62e-05,  s^2 오차 -7.70e-05,  호출  49

짐작이 맞았다. \(\hat\mu\)와 \(\hat\sigma^2\)이 닫힌 해와 비트까지 같고, data.var()가 MLE와 같다는 것도 True다. 표를 보면 사정이 분명해진다. 출발점이 (mean, log var)인 줄만 오차가 정확히 \(0\)이고 함수 호출이 \(45\)번으로 가장 적다. 그 자리에서 출발해 아무 데도 가지 않은 것이다. best_result가 최솟값으로 그 줄을 고르니 최종 답이 닫힌 해와 같아진다.

나머지 다섯 줄이 넬더-미드의 참모습이다. \(\hat\mu\)의 오차가 \(10^{-6}\sim10^{-5}\), \(\hat\sigma^2\)의 오차가 \(5\times10^{-5}\sim9\times10^{-5}\)다. 기본 허용오차 xatol = fatol = 1e-4에 걸맞은 정확도이며, 소수 넷째 자리까지만 믿을 수 있다는 뜻이다. 처음 출력이 넷째 자리까지 맞아 보인 것은 운이 아니라 반칙이었다.

여기서 배울 것은 두 가지다. 첫째, 수치 해를 닫힌 해와 견줄 때는 출발점이 답을 알고 있지 않은지 먼저 확인해야 한다. 둘째, 알고리즘의 정확도를 알고 싶으면 일치하는 자릿수를 세지 말고 허용오차 설정을 보라. \(10^{-4}\)를 넣고 \(10^{-10}\)의 정확도를 기대할 수는 없다.

(1) 해석적으로. 국소 최적이 아니다 — 라벨 바꿈이다. 혼합밀도

에서 치환

을 하면 두 항이 자리만 바꿀 뿐 합이 한 치도 달라지지 않는다. 모든 \(x\)에서 \(f\)가 같으므로 가능도도 정확히 같다. 성분이 \(K\)개면 이런 치환이 \(K!\)가지 있고, 여기서는 \(2! = 2\)다. 곧 전역 최대가 둘이고 둘은 같은 해를 다르게 이름 붙인 것에 지나지 않는다.

출발 2는 \(\mu_1 = \mu_2 = 2\)라는 대칭 자리에서 시작했다. 두 봉우리에서 꼭 같은 거리에 있으므로 어느 쪽으로 가느냐는 단체가 처음 어떻게 펴지느냐가 정할 뿐이고, 어느 쪽으로 가든 같은 높이에 닿는다. 국소 최적에 갇혔다면 음의 로그가능도가 더 컸어야 하는데 셋이 모두 \(402.58\)이다.

\(\hat\pi\) 사이에는 치환이 그대로

을 요구한다. 출력의 \(0.374\)와 \(0.626\)이 정확히 그 관계다.

(2) 수치적으로.

import numpy as np
from scipy import optimize, stats

def mixture_log_likelihood(params, data):
    """성분 둘짜리 정규 혼합모형의 음의 로그가능도."""
    pi, mu1, mu2, sigma = params[0], params[1], params[2], np.exp(params[3])
    pi = 1 / (1 + np.exp(-pi))  # sigmoid transform for mixing weight
    ll = np.sum(np.log(
        pi * stats.norm.pdf(data, mu1, sigma) +
        (1 - pi) * stats.norm.pdf(data, mu2, sigma)
    ))
    return -ll

# 봉우리 둘짜리 혼합분포에서 자료를 만든다.
rng = np.random.default_rng(42)
n = 200
z = rng.binomial(1, 0.4, n)
data = np.where(z, rng.normal(0, 1, n), rng.normal(4, 1, n))

# 출발값을 바꿔 가며 어디로 수렴하는지 본다.
starts = [[0, -1, 5, 0], [0, 2, 2, 0], [0, 0, 3, 0.5]]
for i, x0 in enumerate(starts):
    result = optimize.minimize(mixture_log_likelihood, x0, args=(data,), method="Nelder-Mead")
    pi_hat = 1 / (1 + np.exp(-result.x[0]))
    print(f"Start {i+1}: pi={pi_hat:.3f}, mu1={result.x[1]:.3f}, "
          f"mu2={result.x[2]:.3f}, nll={result.fun:.2f}")

# (1) 이 말한 불변성을 직접 확인한다.
best = optimize.minimize(mixture_log_likelihood, starts[0], args=(data,),
                         method="Nelder-Mead").x
swapped = np.array([-best[0], best[2], best[1], best[3]])   # pi -> 1-pi, mu1 <-> mu2
a = mixture_log_likelihood(best, data)
b = mixture_log_likelihood(swapped, data)

print(f"\nnll(theta)   = {a:.12f}")
print(f"nll(swapped) = {b:.12f}")
print(f"비트까지 같은가? {a == b}")

sig = lambda t: 1 / (1 + np.exp(-t))
print(f"pi + (1-pi) 로 맞물리는가: {sig(best[0]):.6f} + {sig(-best[0]):.6f} "
      f"= {sig(best[0]) + sig(-best[0]):.6f}")
print(f"sigma-hat = {np.exp(best[3]):.4f}  (참값 1)")
print(f"pi-hat = {sig(best[0]):.4f},  표본의 실제 섞임 비율 = {z.mean():.4f},  참값 0.4")

출력:

Start 1: pi=0.374, mu1=0.154, mu2=3.900, nll=402.58
Start 2: pi=0.626, mu1=3.900, mu2=0.154, nll=402.58
Start 3: pi=0.374, mu1=0.154, mu2=3.899, nll=402.58

nll(theta)   = 402.578187662901
nll(swapped) = 402.578187662901
비트까지 같은가? True
pi + (1-pi) 로 맞물리는가: 0.374353 + 0.625647 = 1.000000
sigma-hat = 1.0101  (참값 1)
pi-hat = 0.3744,  표본의 실제 섞임 비율 = 0.3850,  참값 0.4

불변성이 비트까지 확인되었다. 이름을 맞바꾼 모수에서 음의 로그가능도를 다시 재면 \(402.578187662901\)로 자릿수 하나 다르지 않다. 반올림 오차 수준으로 같은 것이 아니라 정확히 같은 수인데, 두 식이 같은 항들을 다른 순서로 더한 것이 아니라 아예 같은 곱의 합이기 때문이다. \(\hat\pi\)도 \(0.374353 + 0.625647 = 1.000000\)으로 맞물린다.

따라서 세 결과는 모두 전역 최대다. 출발 2가 다른 답을 준 것이 아니라 같은 답에 다른 이름을 붙였을 뿐이고, 국소 최적에 갇힌 것과는 성질이 전혀 다르다. 세 nll을 소수 아래까지 펼쳐 보면 \(5\times10^{-8}\) 안에서 같은데, 이는 넬더-미드가 멈춘 자리의 차이일 뿐 봉우리가 달라서 생긴 차이가 아니다.

그렇다고 라벨 바꿈이 무해하지는 않다. 같은 자료에 여러 번 적합하거나 MCMC로 사후표본을 모으면 \(\mu_1\)이 어떤 때는 \(0.15\), 어떤 때는 \(3.90\)으로 튀어 평균을 내면 \(2\) 언저리라는 엉뚱한 값이 나온다. 그래서 실무에서는 \(\mu_1 < \mu_2\) 같은 순서 제약을 걸거나, 적합이 끝난 뒤 성분을 평균순으로 정렬해 이름을 통일한다.

나머지 추정값도 제자리다. \(\hat\sigma = 1.0101\)이 참값 \(1\)에 가깝고, \(\hat\pi = 0.3744\)는 참값 \(0.4\)보다 이 표본에서 실제로 뽑힌 비율 \(0.3850\)에 더 가깝다. 추정량이 겨냥하는 것은 모집단의 참값이지만 손에 쥔 자료는 표본이므로 이것이 당연한 모습이다.

로그가능도는:

\(\ell'(\lambda) = n/\lambda - \sum x_i = 0\)으로 두면 \(\hat{\lambda}_{\text{MLE}} = n/\sum x_i = 1/\bar{X}\)이다.

import numpy as np

rng = np.random.default_rng(42)
data = rng.exponential(scale=2.0, size=50)  # true lambda = 0.5
lam_grid = np.linspace(0.1, 5, 1000)
ll = len(data) * np.log(lam_grid) - lam_grid * data.sum()
lam_hat_grid = lam_grid[np.argmax(ll)]
lam_hat_exact = 1 / data.mean()
print(f"Grid MLE:    {lam_hat_grid:.4f}")
print(f"Analytic MLE: {lam_hat_exact:.4f}")

출력:

Grid MLE:    0.6003
Analytic MLE: 0.5982

두 값이 거의 일치해야 한다. \(\square\)

지수함수 \(\sigma^2 = e^\phi\)는 \(\phi \in \mathbb{R}\)를 \(\sigma^2 \in (0, \infty)\)로 보낸다. 모든 실수 \(\phi\)에 대해 \(e^\phi > 0\)이므로 최적화기는 유효하지 않은(양수가 아닌) 분산을 만들 걱정 없이 \(\mathbb{R}\) 전체를 탐색할 수 있다. 이 재모수화가 없으면 기울기 단계가 \(\sigma^2\)을 0 아래로 밀어낼 수 있고, (정규 로그가능도에 \(\log(\sigma^2)\)이 나오므로) 로그가능도가 정의되지 않게 된다. 이 변환은 곡률을 더 고르게 만들어 최적화 지형도 개선한다. \(\square\)

(상수를 무시한) 로그가능도는:

점수: \(\ell'(p) = x/p - (n-x)/(1-p)\).

Hessian: \(\ell''(p) = -x/p^2 - (n-x)/(1-p)^2\).

Newton-Raphson 갱신: \(p^{(t+1)} = p^{(t)} - \ell'(p^{(t)})/\ell''(p^{(t)})\).

\(n = 20, x = 14, p^{(0)} = 0.5\)일 때:

• \(\ell'(0.5) = 14/0.5 - 6/0.5 = 28 - 12 = 16\)
• \(\ell''(0.5) = -14/0.25 - 6/0.25 = -56 - 24 = -80\)
• \(p^{(1)} = 0.5 - 16/(-80) = 0.5 + 0.2 = 0.7\)

\(p^{(1)} = 0.7\)에서:

• \(\ell'(0.7) = 14/0.7 - 6/0.3 = 20 - 20 = 0\)

따라서 \(p^{(2)} = 0.7\)이며, 이는 이미 MLE \(\hat{p} = x/n = 14/20 = 0.7\)이다. Newton-Raphson이 한 단계 만에 수렴했다. \(\square\)

혼합밀도 \(\pi \cdot N(x_1, \sigma_1^2) + (1-\pi) \cdot N(\mu_2, \sigma_2^2)\)을 생각하자. \(\mu_1 = x_1\)(어떤 자료점)로 두고 \(\sigma_1 \to 0\)으로 보내면 \(x_1\)에서 첫 성분의 밀도가 \(1/\sigma_1 \to \infty\)로 발산하여 로그가능도가 유계가 아니게 된다.

실무에서 이것이 문제가 되지 않는 이유는:

1. 이런 퇴화된 해는 자료점 하나에 과적합한 것이어서 통계적으로 무의미하다.
2. 혼합모형의 표준 방법인 EM 알고리즘은 합리적인 출발값에서 그런 퇴화된 해에 도달하지 못한다.
3. 실무자는 최소 분산 제약을 두거나 벌점 가능도를 쓴다.
4. 유용한 MLE는 전역 상한이 아니라 가능도의 국소 최댓값이다. \(\square\)

import numpy as np

rng = np.random.default_rng(42)
n = 50
data = rng.binomial(1, 0.3, n)
x_sum = data.sum()

# 뉴턴-랩슨: 관측정보량을 쓴다.
p_nr = 0.5
for i in range(10):
    score = x_sum / p_nr - (n - x_sum) / (1 - p_nr)
    hessian = -x_sum / p_nr**2 - (n - x_sum) / (1 - p_nr)**2
    p_nr = p_nr - score / hessian
    print(f"NR  iter {i+1}: p = {p_nr:.8f}")

# 피셔 스코어링: 기대정보량을 쓴다.
p_fs = 0.5
for i in range(10):
    score = x_sum / p_fs - (n - x_sum) / (1 - p_fs)
    fisher_info = n / (p_fs * (1 - p_fs))
    p_fs = p_fs + score / fisher_info
    print(f"FS  iter {i+1}: p = {p_fs:.8f}")

출력:

NR  iter 1: p = 0.36000000
NR  iter 2: p = 0.36000000
NR  iter 3: p = 0.36000000
NR  iter 4: p = 0.36000000
NR  iter 5: p = 0.36000000
NR  iter 6: p = 0.36000000
NR  iter 7: p = 0.36000000
NR  iter 8: p = 0.36000000
NR  iter 9: p = 0.36000000
NR  iter 10: p = 0.36000000
FS  iter 1: p = 0.36000000
FS  iter 2: p = 0.36000000
FS  iter 3: p = 0.36000000
FS  iter 4: p = 0.36000000
FS  iter 5: p = 0.36000000
FS  iter 6: p = 0.36000000
FS  iter 7: p = 0.36000000
FS  iter 8: p = 0.36000000
FS  iter 9: p = 0.36000000
FS  iter 10: p = 0.36000000

둘 다 \(\hat{p} = x_{\text{sum}}/n\)으로 수렴한다. 베르누이에서는 관측 정보량과 기대 정보량이 밀접하게 연결되어 있어 수렴 속도가 거의 같다. 일반적으로는 관측 Hessian의 조건수가 나쁠 때 Fisher 점수법이 더 안정적일 수 있다. \(\square\)

(1) 목적함수의 변화. \(|\ell^{(k+1)}-\ell^{(k)}| < \varepsilon\).

• 실패: 목적함수가 평평한 고원을 지날 때 변화가 작아 조기에 멈춘다. 실제로는 아직 봉우리에서 멀 수 있다.
• 상대 기준 \(|\Delta\ell|/(|\ell|+1) < \varepsilon\)을 쓰는 편이 척도에 덜 민감하다.

(2) 모수의 변화. \(\|\boldsymbol\theta^{(k+1)}-\boldsymbol\theta^{(k)}\| < \varepsilon\).

• 실패: 보폭이 작게 잡혀 천천히 기어갈 때도 만족된다. 감쇠가 과한 알고리즘에서 흔하다.
• 모수마다 척도가 다르면 상대 변화로 재야 한다.

(3) 기울기의 크기. \(\|\nabla\ell(\boldsymbol\theta^{(k)})\| < \varepsilon\).

• 가장 원리적이다. 내부 최댓값의 필요조건이 기울기 0이기 때문이다.
• 실패: 최댓값이 경계에 있으면 기울기가 0이 아니어서 영원히 만족되지 않는다. 또 기울기의 크기는 목적함수의 척도에 의존하므로 절대 기준을 정하기 어렵다.

권고. 세 기준을 모두 요구하고, 여기에 다음을 더한다.

• 헤시안 확인. 수렴한 점에서 헤시안이 음정부호인지 본다. 그렇지 않으면 안장점이거나 수렴하지 않은 것이다. 표준오차를 구하려면 어차피 계산해야 한다.
• 여러 초기값. 서로 다른 출발점에서 같은 곳으로 가는지 본다. 다봉성을 잡아내는 가장 확실한 방법이다.
• 최대 반복 횟수 도달 여부. 수렴 플래그를 반드시 확인한다. 많은 라이브러리가 조용히 최대 횟수에서 멈추고 결과를 돌려준다.

실무의 사고. scipy.optimize.minimize의 결과에서 res.success를 확인하지 않고 res.x만 쓰는 것이 흔한 실수다. 수렴하지 않은 값을 추정값으로 보고하게 된다.

(1) 재모수화. 제약을 자동으로 만족하는 변환을 쓴다.

• 장점: 제약 없는 최적화기를 그대로 쓸 수 있다. 로그가능도가 새 척도에서 오히려 더 이차식에 가까워지는 경우도 많다.
• 단점: 최적해가 경계에 있으면 \(\phi \to -\infty\)로 발산한다. 표준오차를 원래 척도로 되돌리려면 델타 방법이 필요하다.

(2) 경계 제약 최적화. L-BFGS-B처럼 상자 제약을 지원하는 알고리즘을 쓴다.

• 장점: 직관적이고, 경계해가 나와도 그 값을 그대로 준다.
• 단점: 선형·비선형 제약(예: \(\sum\pi_k=1\))은 상자로 표현할 수 없다. 경계에서 멈추면 표준오차의 정규근사가 성립하지 않는다.

(3) 벌점이나 사전분포. 목적함수에 \(-\gamma/\sigma^2\) 같은 항을 더하거나 베이즈 사전분포를 둔다.

• 장점: 경계로 가는 퇴화를 근본적으로 막는다. 혼합모형의 분산 퇴화를 다루는 표준적인 방법이다.
• 단점: 목적함수가 달라지므로 더 이상 순수한 MLE가 아니다. 벌점의 세기를 정해야 한다.

권고. 대부분의 경우 재모수화가 가장 깔끔하다. 최적해가 경계에 있을 가능성이 실질적이면(분산성분이 0일 수 있는 혼합효과 모형 등) 경계 제약 최적화를 쓰고, 경계해가 나왔을 때 표준적인 추론을 쓰면 안 된다는 점을 기억한다. 그 경우 우도비 통계량의 극한분포가 \(\chi^2\)가 아니라 \(\chi^2\)들의 혼합이 된다.

확인.

• 여러 초기값에서 돌린다. 수렴한 \(\ell\) 값들을 히스토그램으로 보면 봉우리가 몇 개인지 드러난다. 값이 여러 무리로 갈리면 다봉성이다.
• 모수가 1~2개면 격자로 그린다. 가장 확실하다. 3개 이상이면 두 개씩 골라 프로파일 곡면을 그린다.
• 수렴점 사이를 잇는 경로를 따라 \(\ell\)을 계산한다. 두 해 사이에 골짜기가 있으면 서로 다른 봉우리다.

전역 최댓값을 찾는 전략.

• 다중 출발(multistart). 모수공간에서 무작위로(또는 라틴 초입방 설계로) 초기값을 여러 개 뽑아 각각 국소 최적화하고 가장 좋은 것을 고른다. 가장 널리 쓰이고 병렬화가 쉽다.
• 좋은 초기값을 설계한다. 적률법 추정값, 더 단순한 모형의 해, 자료를 나눠 각각 적합한 값 등. 혼합모형이라면 \(k\)-평균 군집 결과로 초기화하는 것이 표준이다.
• 전역 최적화 알고리즘. basinhopping, differential_evolution, 시뮬레이티드 어닐링. 느리지만 넓게 훑는다.
• 점진적 접근. 단순한 모형에서 시작해 복잡도를 조금씩 올리며 직전 해를 초기값으로 쓴다. 성분 수를 1부터 늘려 가는 혼합모형 적합이 그 예다.
• 결정론적 어닐링. 목적함수를 처음에는 매끄럽게 만들었다가 점차 원래대로 되돌린다.

보고할 때. 전역 최댓값을 찾았다고 증명할 수는 없다. 정직하게 하려면 몇 개의 초기값에서 출발했고 몇 개의 서로 다른 국소해를 찾았는지 기록한다. 국소해들의 \(\ell\) 값이 비슷하면서 모수가 크게 다르다면, 그것은 모형이 식별되지 않는다는 신호일 수 있으므로 최적화 문제로만 볼 것이 아니다.

대체로 그대로 쓰면 안 된다.

이유 1 — 준뉴턴 근사는 정확하지 않다. BFGS가 유지하는 hess_inv는 최적화를 안내하기 위한 근사이며, 수렴 경로에서 본 기울기 변화로 쌓아 올린 것이다. 봉우리에서의 참 곡률과는 다를 수 있고, 특히 반복 횟수가 적으면 초기 근사(보통 단위행렬)의 흔적이 남는다.

이유 2 — 척도 문제. 재모수화해서 최적화했다면 헤시안도 그 척도의 것이다. 원래 모수의 표준오차를 얻으려면 야코비안으로 변환해야 한다.

이것을 잊는 것이 가장 흔한 실수다.

이유 3 — 목적함수의 부호와 상수. 최소화를 위해 음의 로그가능도를 넘겼다면 그 헤시안이 곧 관측정보량이다. 부호를 한 번 더 뒤집으면 음수 분산이 나온다. 또 \(-2\ell\)을 최소화했다면 헤시안이 2배이므로 표준오차가 \(\sqrt2\)배 틀린다.

권장 절차.

1. 최적화는 아무 방법으로나 한다.
2. 수렴한 \(\hat{\boldsymbol\theta}\)에서 헤시안을 다시 계산한다. 자동미분이 가장 좋고, 없으면 numdifftools 같은 고정밀 수치 미분을 쓴다.
3. 그 헤시안이 양정부호인지 확인한다(음의 로그가능도 기준). 아니면 수렴하지 않았거나 안장점이다.
4. 역행렬의 대각 제곱근이 표준오차다. 필요하면 델타 방법으로 원 척도로 옮긴다.
5. 가능하면 프로파일 가능도나 부트스트랩으로 교차 확인한다. 왈드 표준오차와 크게 다르면 이차 근사가 나쁘다는 뜻이다.

(1) 확률적 경사법. 매 단계에서 자료 전체가 아니라 무작위로 고른 작은 묶음(minibatch)으로 기울기를 근사한다.

• 한 단계의 비용이 \(n\)과 무관하다.
• 학습률 \(\eta_k\)가 \(\sum\eta_k=\infty\), \(\sum\eta_k^2<\infty\)를 만족하면 수렴이 보장된다(로빈스-먼로).
• 기계학습의 표준이며, Adam 같은 적응적 변형이 널리 쓰인다.
• 단점: 표준오차를 곧바로 주지 않고, 수렴 판정이 까다롭다.

(2) 충분통계량 활용. 지수족이면 자료 전체가 아니라 저차원 요약만 있으면 된다. 정규분포는 \((\sum x_i, \sum x_i^2)\), 포아송은 \(\sum x_i\), 로지스틱 회귀는 \(X^\top\mathbf{y}\)와 \(X^\top X\)면 충분하다.

• 자료를 한 번만 훑으며 통계량을 누적하면 그 뒤로는 \(n\)과 무관한 비용으로 최적화한다.
• 분산 처리도 쉽다. 각 노드에서 통계량을 계산해 더하기만 하면 된다.
• 단점: 지수족이 아니면 쓸 수 없다.

(3) 표본추출·분할정복.

• 부분표본. \(n=10^7\)이면 \(10^5\)만 써도 표준오차가 \(\sqrt{100}=10\)배 커질 뿐이다. 이미 충분히 정밀하다면 나머지를 쓸 이유가 없다. 정밀도가 이미 실무 요구를 넘었는지 먼저 따져 볼 일이다.
• 분할정복. 자료를 \(K\)조각으로 나눠 각각 \(\hat\theta_j\)를 구하고 평균한다. 앞서 본 정보량 가중을 쓰면 된다. 일차적으로는 전체 MLE와 점근적으로 동등하다.
• 중요도 가중 부분표본. 영향력이 큰 관측을 더 자주 뽑고 가중치로 보정한다. 희귀사건 자료에서 특히 효과적이다.

고르는 기준. 모형이 지수족이면 (2)가 압도적으로 좋다. 그렇지 않고 정밀도가 이미 충분하면 (3)의 부분표본이 가장 간단하다. 복잡한 모형을 전체 자료로 적합해야 한다면 (1)을 쓰되, 표준오차는 부트스트랩이나 샌드위치로 따로 구한다.