베르누이분포의 MLE¶

(1) 해석적으로. \(k = \sum_i x^{(i)} = 74\)라 두면

이고 미분해 \(0\)으로 두면

이다. 한 번 더 미분하면

이 \((0,1)\) 전체에서 성립하므로 \(J\)가 엄밀히 볼록이고 정류점이 유일한 전역 최소다. 가능도 쪽에서 보면 \(\ell = -J\)가 엄밀히 오목이라는 같은 말이다.

최솟값이 무엇인가. \(\hat p = k/m\)을 넣으면 \(k = m\hat p\)이므로

곧 관측 하나당 비용의 바닥이 \(\text{Bernoulli}(\hat p)\)의 엔트로피다. \(\hat p = 0.74\)에서

이므로 \(J(\hat p) = 57.3057\), \(\ell(\hat p) = -57.3057\)이다.

이것이 교차엔트로피 손실을 쓸 때 꼭 알아 두어야 할 사실이다. 손실은 아무리 잘 맞춰도 \(0\)으로 내려가지 않는다. 바닥은 자료 자체의 엔트로피이며, 그보다 낮은 값이 나왔다면 자료를 외웠다는 뜻이다. 결정론적인 자료(\(\hat p = 0\) 또는 \(1\))일 때만 \(H = 0\)이 되어 손실이 \(0\)에 닿는다.

어긋남의 예고. 격자 간격이 \((0.99 - 0.01)/99 = 0.0098990\)이고 \((0.74 - 0.01)/0.0098990 = 73.745\)가 정수가 아니므로 \(0.74\)는 격자에 없다. 가장 가까운 격자점은 \(j = 74\)인 \(0.7425253\)이고, \(\lvert\ell''(\hat p)\rvert = m/(\hat p(1-\hat p)) = 519.75\)이므로 로그가능도가

만큼 낮게 나와야 한다. 곧 \(-57.3057 - 0.0017 = -57.3074\)다. 출력에 찍힌 값과 같다.

(2) 수치적으로.

import numpy as np
import matplotlib.pyplot as plt

# 같은 결과가 다시 나오도록 난수 씨앗을 고정한다.
seed = 1
np.random.seed(seed)

# 앞면 확률과 던지는 횟수를 정한다.
p = 0.7
n_samples = 100

def load_data():
    """
    Simulate coin flips based on a binomial distribution.

    Returns:
    - numpy array: Array of coin flips (1 for heads, 0 for tails).
    """
    return np.random.binomial(n=1, p=p, size=(n_samples,))  # Shape (100,)

def compute_prob(coin, p):
    """
    Compute the probability of a single coin flip outcome.

    Parameters:
    - coin: Outcome of the coin flip (1 for heads, 0 for tails).
    - p: Probability of heads.

    Returns:
    - float: Probability of observing the outcome.
    """
    return p**coin * (1 - p)**(1 - coin)

def compute_log_prob(coin, p):
    """
    Compute the log-probability of a single coin flip outcome.

    Parameters:
    - coin: Outcome of the coin flip (1 for heads, 0 for tails).
    - p: Probability of heads.

    Returns:
    - float: Log-probability of observing the outcome.
    """
    return coin * np.log(p) + (1 - coin) * np.log(1 - p)

def compute_likelihood(coins, p):
    """
    Compute the joint probability of all coin flips for a given probability.

    Parameters:
    - coins: Array of coin flip outcomes.
    - p: Probability of heads.

    Returns:
    - float: Joint probability of observing all outcomes.
    """
    joint_prob = 1.0
    for coin in coins:
        joint_prob *= compute_prob(coin, p)
    return joint_prob

def compute_log_likelihood(coins, p):
    """
    Compute the log-likelihood of all coin flips for a given probability.

    Parameters:
    - coins: Array of coin flip outcomes.
    - p: Probability of heads.

    Returns:
    - float: Log-likelihood of observing all outcomes.
    """
    log_joint_prob = 0.0
    for coin in coins:
        log_joint_prob += compute_log_prob(coin, p)
    return log_joint_prob

# 동전을 n번 던진다. 참 p 는 우리가 모르는 값이라고 둔다.
coins = load_data()

# p 후보를 격자로 늘어놓는다. 이 중 로그가능도가 가장 큰 것을 고른다.
ps = np.linspace(0.01, 0.99, 100)

# 후보마다 로그가능도를 계산한다.
log_likelihood_list = [compute_log_likelihood(coins, p) for p in ps]
log_likelihood = np.array(log_likelihood_list)

# 로그가능도가 가장 큰 후보가 최대가능도추정값이다.
idx = np.argmax(log_likelihood)
mle_p = ps[idx]
log_likelihood_max = log_likelihood[idx]
print(f"MLE index: {idx}")
print(f"MLE probability (p): {mle_p:.4f}")
print(f"Max log-likelihood: {log_likelihood_max:.4f}\n")

# 로그가능도 곡선을 그리고 최댓값 자리를 표시한다.
fig, ax = plt.subplots(figsize=(12, 3))
ax.plot(ps, log_likelihood, label="Log-likelihood")
ax.plot([mle_p, mle_p], [0, log_likelihood_max], '--or', label="MLE")
ax.legend(loc="lower right")

# 축 이름과 범례를 다듬는다.
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_position("zero")
ax.spines['left'].set_position("zero")
ax.set_xlabel("Probability (p)")
ax.set_ylabel("Log-likelihood")
plt.show()

# (1) 에서 따진 것을 확인한다.
k = coins.sum()
m = len(coins)
p_exact = k / m
exact_ll = k * np.log(p_exact) + (m - k) * np.log(1 - p_exact)
entropy = -(p_exact * np.log(p_exact) + (1 - p_exact) * np.log(1 - p_exact))
step = ps[1] - ps[0]

print(f"k = {k},  m = {m},  해석적 p-hat = k/m = {p_exact}")
print(f"격자 간격 = {step:.7f},  (0.74 - 0.01)/간격 = {(0.74 - 0.01) / step:.3f}")
print(f"해석적 l(p-hat) = {exact_ll:.6f}")
print(f"격자   l({mle_p:.7f}) = {log_likelihood_max:.6f}")
print(f"어긋남 = {exact_ll - log_likelihood_max:.6f}   "
      f"(예측 {0.5 * m / (p_exact * (1 - p_exact)) * (mle_p - p_exact)**2:.6f})")
print(f"J(p-hat)/m = {-exact_ll / m:.7f},  H(p-hat) = {entropy:.7f}")

# 가능도를 곱으로 그대로 계산해도 로그를 취하면 같은 값이 나온다.
raw = compute_likelihood(coins, p_exact)
print(f"곱으로 계산한 L(p-hat) = {raw:.6e},  log = {np.log(raw):.6f}")

출력:

MLE index: 74
MLE probability (p): 0.7425
Max log-likelihood: -57.3074

k = 74,  m = 100,  해석적 p-hat = k/m = 0.74
격자 간격 = 0.0098990,  (0.74 - 0.01)/간격 = 73.745
해석적 l(p-hat) = -57.305692
격자   l(0.7425253) = -57.307356
어긋남 = 0.001664   (예측 0.001657)
J(p-hat)/m = 0.5730569,  H(p-hat) = 0.5730569
곱으로 계산한 L(p-hat) = 1.295550e-25,  log = -57.305692

세 가지가 모두 맞는다.

첫째, \(J(\hat p)/m = 0.5730569\)와 \(H(\hat p) = 0.5730569\)가 소수 일곱째 자리까지 같다. (1)에서 유도한 "교차엔트로피 손실의 바닥은 경험분포의 엔트로피"가 수로 확인되었다.

둘째, 격자와 해석적 답의 어긋남이 \(0.001664\)인데 포물선 근사로 예고한 \(0.001657\)과 셋째 자리까지 맞는다. 남는 \(7\times10^{-6}\)은 포물선이 근사일 뿐이라서 생기는 삼차 이상의 항이다. 격자 탐색이 \(0.7425\)를 준 것은 \(0.74\)가 격자에 없기 때문이며, \((0.74-0.01)/0.0098990 = 73.745\)가 정수가 아니라는 것이 그 전부다.

셋째, 가능도를 \(100\)개의 곱으로 그대로 계산한 \(1.2956\times10^{-25}\)에 로그를 취하면 \(-57.305692\)로 로그가능도와 정확히 같다. \(m = 100\)에서는 아직 언더플로가 없어 두 길이 같은 답을 준다. 그러나 이 값이 이미 \(10^{-25}\)인 데서 보듯, \(m\)을 열 배로 늘리면 \(10^{-250}\) 자리가 되고 그 뒤로는 곱 쪽이 무너진다. compute_likelihood와 compute_log_likelihood를 나란히 정의해 둔 까닭이 이것이다.

\(n = 20\)번 중 \(k = 13\)번이 앞면이라 하자. 로그가능도는:

도함수를 0으로 두면:

2계도함수가 \(-13/p^2 - 7/(1-p)^2 < 0\)이므로 최댓값임이 확인된다.

\(k = \sum x_i\)번 성공한 \(n\)번의 독립 베르누이 시행에 대해 로그가능도는:

미분하여 0으로 두면:

이는 \(0 \leq k \leq n\)인 임의의 \(k\)와 \(n\)에 대해 성립한다. \(\square\)

Bernoulli\((p)\) 관측값 하나에 대해 로그가능도는 \(\ell(p) = x\log p + (1-x)\log(1-p)\)이다. 2계도함수는:

(\(E[X] = p\)를 사용하여) 기댓값에 음수를 취하면:

\(n\)개의 관측값에 기반한 \(\hat{p}\)의 점근분산은:

표본비율의 분산에 대한 익숙한 공식이다.

MLE \(\hat{p} = 0\)은 이 동전에서 앞면이 결코 나올 수 없다는 뜻인데, 관측값 10개만으로 내리기에는 극단적인 결론이다. 문제는 특히 소표본에서 MLE가 비합리적인 경계값 추정치를 낼 수 있다는 것이다.

한 가지 대안은 라플라스 평활(또는 균등 사전분포를 쓴 베이즈 접근)이다. 자료에 "가상의 성공" 하나와 "가상의 실패" 하나를 더하여 \(\hat{p}_{\text{Laplace}} = (0+1)/(10+2) = 1/12 \approx 0.083\)을 얻는다. 자료에 근거하면서도 0이라는 추정값을 피한다. 더 형식적으로 이는 Beta\((1,1)\)(균등) 사전분포 아래의 사후평균에 해당하며 \(\hat{p}_{\text{Bayes}} = (k+1)/(n+2)\)를 준다.

\(\hat p = 0.65\)이다.

왈드 검정. MLE에서 평가한 정보량을 쓴다.

점수 검정. 귀무가설 아래의 정보량을 쓴다.

우도비 검정.

\(\Lambda \sim \chi^2_1\)이므로 \(p\text{-값} = 0.176\)이다.

비교. 셋 다 5%에서 기각하지 못하고 \(p\)-값이 0.16~0.18로 비슷하다. 참고로 정확 이항검정은 0.263으로 꽤 다르다.

세 검정의 관계. 모두 \(n\to\infty\)에서 점근적으로 동등하며 \(H_0\) 아래에서 \(\chi^2_1\)을 따른다. 차이는 로그가능도 곡선을 어디서 재느냐에 있다.

실무적 권고는 우도비 검정이다. 이유는 두 가지다. 첫째, 모수화에 불변이다. \(p\)로 검정하든 로그오즈로 검정하든 같은 값이 나온다. 왈드는 그렇지 않아서 \(p\)로 하느냐 \(\ln\{p/(1-p)\}\)로 하느냐에 따라 결과가 달라진다. 둘째, 로그가능도가 이차식에서 멀어질수록 왈드가 부정확해지는데(경계 근처에서 특히 심하다) 우도비는 훨씬 잘 버틴다.

점수 검정은 \(\hat\theta\)를 구하지 않아도 계산할 수 있다는 장점이 있어, 모형을 적합하기 전에 변수 추가의 효과를 미리 보는 데 쓰인다.

프로파일 가능도 구간. 우도비 검정을 뒤집어

을 구간으로 삼는다. \(k=0\)이면 \(\ell(p) = 10\ln(1-p)\)이고 \(\hat p = 0\)에서 \(\ell(\hat p) = 0\)이므로

이다. 구간은 \((0,\ 0.175)\)다.

비교.

읽는 법. 프로파일 가능도 구간이 가장 좁고 클로퍼-피어슨이 가장 넓다. 클로퍼-피어슨은 포함확률이 명목값 이상임을 보장하느라 보수적이고, 프로파일 구간은 카이제곱 근사에 기대므로 이 극단적인 상황에서는 다소 낙관적이다.

핵심은 왈드만 빼고 셋 다 쓸 만하다는 점이다. 왈드가 무너지는 이유는 로그가능도를 \(\hat p\) 근처에서 이차함수로 근사하는데, \(\hat p = 0\)이 모수공간의 경계라 그 근사가 성립하지 않기 때문이다. 실제로 \(k=0\)일 때 로그가능도 \(10\ln(1-p)\)는 \(p=0\)에서 기울기가 \(-10\)으로 0이 아니다. 봉우리가 아니라 절벽 끝인 셈이다.

일반 원리. 추정값이 모수공간의 경계에 있거나 로그가능도가 심하게 비대칭이면 왈드 방법을 쓰면 안 된다. 프로파일 가능도나 정확 방법으로 옮겨야 한다.

MLE. 불변성에 따라

표준오차. \(g(p) = \ln\{p/(1-p)\}\)의 도함수가 \(g'(p) = 1/\{p(1-p)\}\)이므로

두 구간.

(\(p\) 척도에서 바로)

(\(\theta\) 척도에서 만들어 되돌림)

이고 로지스틱 변환 \(p = e^\theta/(1+e^\theta)\)로 되돌리면

다르다. 두 번째 구간은 비대칭이며(\(\hat p\)에서 아래로 0.224, 위로 0.173) 언제나 \((0,1)\) 안에 머문다. 첫 번째는 대칭이고 \(\hat p\)가 극단적이면 \([0,1]\)을 벗어난다.

어느 쪽이 나은가. 일반적으로 로그오즈 척도에서 만드는 쪽이 낫다. 로그가능도가 \(\theta\) 척도에서 이차함수에 훨씬 가깝기 때문이고, 이것이 왈드 근사의 전제다. 모수공간이 \(\mathbb{R}\) 전체로 펴져 경계 문제도 사라진다.

일반 교훈. 왈드 구간은 모수화에 의존한다. 분산, 비율, 위험비처럼 경계가 있거나 치우친 모수는 로그나 로짓 척도에서 구간을 만들고 되돌리는 것이 표준적인 관행이다. 반면 우도비 구간은 단조변환에 불변이라 이런 고민이 없다.

(가) 포화모형. 각 병원이 독립이므로 로그가능도가 분리되고

이다. 모수가 3개다.

(나) 축소모형. 공통 \(p\)에 대해

이므로

로 전체를 합친 비율이다. 모수가 1개다.

우도비 통계량.

자유도는 모수 개수의 차이인 \(3-1 = 2\)이므로 \(\Lambda \sim \chi^2_2\)이다.

알아 둘 점.

• 이 \(\Lambda\)를 이탈도(deviance) 라 부르며, 일반화선형모형에서 내포된 모형을 비교하는 표준 도구다. 로지스틱 회귀에서 "병원"을 요인으로 넣은 모형과 절편만 있는 모형을 비교하는 것과 같다.
• 같은 자료에 카이제곱 동질성 검정을 적용해도 되고, 두 통계량은 점근적으로 동등하다. 우도비 쪽을 \(G^2\), 피어슨 쪽을 \(X^2\)이라 부른다.
• 기대도수가 작으면 두 근사 모두 나빠진다. 그때는 정확검정이나 순열검정을 쓴다.

사후분포는 \(\text{Beta}(a+k,\ b+n-k)\)이다.

사후평균.

MAP. \(\text{Beta}(\alpha,\beta)\)의 최빈값이 \(\alpha,\beta>1\)일 때 \((\alpha-1)/(\alpha+\beta-2)\)이므로

MLE. \(\hat p = k/n\).

관계.

• \(a=b=1\)(균등 사전분포)이면 \(\hat p_{\text{MAP}} = k/n = \hat p_{\text{MLE}}\)다. 균등 사전분포의 MAP는 MLE와 같다. 사후분포가 가능도에 비례하므로 당연하다. 반면 사후평균은 \((k+1)/(n+2)\)로 다르다(라플라스 평활).
• \(a=b=1/2\)(제프리스 사전분포)이면 \(\hat p_{\text{MAP}} = (k-0.5)/(n-1)\), 사후평균은 \((k+0.5)/(n+1)\)이다.
• \(n\to\infty\)이면 셋 모두 \(k/n\)으로 수렴한다. 사전분포의 영향이 \(O(1/n)\)로 사라진다.

세 값이 모두 일치하는 조건. MAP와 사후평균이 같으려면 사후분포가 대칭이어야 하고, 베타분포는 두 모수가 같을 때만 대칭이므로 \(a+k = b+n-k\)가 필요하다. 여기에 MLE까지 맞추려면 \(a=b=1\)이어야 하고, 그러면 \(k = n-k\), 즉 \(k = n/2\)다. 즉 균등 사전분포에서 정확히 절반이 성공했을 때만 셋이 모두 \(1/2\)로 일치한다.

해석. MAP는 "사후분포의 봉우리", 사후평균은 "사후분포의 무게중심"이다. 사후분포가 치우쳐 있으면 둘이 갈라지며, 치우침이 심한 경계 근처에서 차이가 크다. \(k=0\), \(n=10\), 균등 사전분포라면 MAP는 0이지만 사후평균은 \(1/12\)다. MAP는 MLE의 문제(경계값)를 물려받지만 사후평균은 그렇지 않다.

어느 쪽을 쓸지는 손실함수가 정한다. 제곱오차 손실이면 사후평균, 0-1 손실이면 MAP, 절대오차 손실이면 사후중앙값이 최적이다.

왜 존재하지 않는가. 로지스틱 회귀에서 \(P(Y=1\mid x) = \sigma(\beta_0+\beta_1x)\)이다. 완전분리가 있으면 기울기 \(\beta_1\)을 키울수록 모든 관측의 예측확률이 관측된 값(0 또는 1)에 가까워지므로 가능도가 계속 커진다.

그러나 이 상한은 달성되지 않는다. 로그가능도가 위로 유계이지만 최댓값이 무한대에서만 접근되므로 MLE가 존재하지 않는다. 수치 최적화는 수렴하지 않고 계수와 표준오차가 발산한다.

증상. 소프트웨어가 "완전분리" 경고를 내거나, 계수 추정값이 \(\pm20\) 이상으로 터무니없이 크고 표준오차가 더 크게 나온다(\(\hat\beta = 15\), \(\operatorname{SE} = 3000\) 식). 왈드 검정의 \(z = \hat\beta/\operatorname{SE}\)가 0으로 가서 무한히 강한 효과인데 유의하지 않다고 나오는 역설적 결과가 생긴다.

대처.

• 펄스 보정(Firth 방법). 가능도에 제프리스 사전분포에 해당하는 벌점 \(|I(\boldsymbol\beta)|^{1/2}\)을 곱한다. 유한한 추정값이 언제나 존재하고, 덤으로 소표본 편향도 줄여 준다. 현재 표준 권고이며 R의 logistf, Python의 여러 구현이 있다.
• 약한 정보 사전분포. 계수에 \(t_1(0, 2.5)\) 같은 사전분포를 두는 베이즈 접근이다. 실질적으로 펄스 방법과 비슷한 역할을 한다.
• 벌점회귀. 능형(L2)이나 라소(L1)를 쓰면 계수가 유한하게 눌린다.
• 왈드 대신 프로파일 가능도 구간. 계수가 발산해도 프로파일 구간은 한쪽이 무한인 형태로 의미 있는 답을 준다.

하지 말아야 할 것. 문제를 일으키는 변수를 그냥 빼는 것은 최악이다. 완전분리는 그 변수가 결과를 완벽히 예측한다는 뜻이고, 이는 통계적 문제이기 이전에 정보다. 그 변수가 결과의 일부이거나 결과 이후에 측정된 것은 아닌지(정보 누출) 먼저 확인해야 한다.