확률변수로서의 통계량¶

(1) 해석적으로. 불편성은 대칭성만으로 나온다. 공에 적힌 수를 뒤집는 사상 \(b \mapsto 32 - b\)는 모집단 \(\{0, 1, \dots, 32\}\)를 자기 자신으로 옮기는 일대일 대응이므로, 어떤 표본 \(S\)를 뽑을 확률과 뒤집은 표본 \(32 - S\)를 뽑을 확률이 같다. 그런데 중앙값은 이 뒤집기를 그대로 따라간다.

따라서 \(\text{Med}\)와 \(32 - \text{Med}\)는 같은 분포를 갖고, 양변의 기댓값을 취하면 \(E[\text{Med}] = 32 - E[\text{Med}]\), 곧

이다. 근사가 아니라 정확한 등식이며, 중앙값이라는 사실조차 쓰지 않았다. 뒤집기와 함께 움직이는 통계량이면 무엇이든 같은 논증이 통한다.

퍼짐은 분포를 직접 적어야 나온다. 중앙값이 \(m\)이려면 뽑힌 다섯 개 가운데 \(m\)보다 작은 것이 둘, 큰 것이 둘이어야 하므로

이다. 이 분포로 계산하면 \(\operatorname{Var}(\text{Med}) = 34\)가 딱 떨어져

이다. 견주어 둘 것이 하나 있다. 같은 표본에서 표본평균을 쓰면 모분산이 \(\sigma^2 = (33^2-1)/12 = 90.67\)이고 비복원이므로 유한모집단 수정이 붙어

이다. 둘 다 불편인데 중앙값 쪽이 \(1.46\)배 더 흔들린다. 불편성은 추정량을 고르는 기준이 되지 못한다는 것을 여기서 이미 볼 수 있다.

(2) 수치적으로. 정확한 분포를 완전열거로 얻고, 그 옆에 50번 모의실험을 나란히 둔다.

import matplotlib.pyplot as plt
import numpy as np
from math import comb

np.random.seed(0)
num_samples = 50

def exact_median_pmf():
    """크기 5짜리 비복원표본에서 표본중앙값의 **정확한** 분포.

    중앙값이 m이려면 m보다 작은 공 2개와 큰 공 2개가 함께 뽑혀야 하므로
    확률이 C(m,2)C(32-m,2)/C(33,5)이다. 모의실험이 아니라 완전열거다.
    """
    total = comb(33, 5)
    return {m: comb(m, 2) * comb(32 - m, 2) / total for m in range(2, 31)}

def main():
    # 모집단: 0부터 32까지 번호가 붙은 공 33개.
    # 개수가 홀수이므로 참 중앙값이 정확히 16으로 딱 떨어진다.
    balls = np.arange(33)
    print(f"모집단 중앙값 = {np.median(balls)}")

    pmf = exact_median_pmf()
    exact_mean = sum(m * q for m, q in pmf.items())
    exact_sd = sum((m - exact_mean) ** 2 * q for m, q in pmf.items()) ** 0.5
    print(f"정확한 E[Med] = {exact_mean:.6f},  sd(Med) = {exact_sd:.6f}")

    # 크기 5짜리 표본을 50번 뽑아 그때마다 표본중앙값을 기록한다.
    # 표본이 달라지면 중앙값도 달라진다는 것,
    # 즉 **통계량이 확률변수라는 것**이 이 보기의 전부다.
    data = []
    for _ in range(num_samples):
        sample = np.random.choice(balls, size=5, replace=False)
        data.append(np.median(sample))

    se = exact_sd / np.sqrt(num_samples)
    print(f"50회 표본중앙값의 평균 = {np.mean(data):.2f}")
    print(f"  그 평균의 표준오차 = {se:.4f},  z = {(np.mean(data) - 16) / se:+.3f}")

    # 되풀이를 20만 번으로 늘려 본다. 같은 표본에서 표본평균도 함께 기록해
    # 중앙값과 평균 가운데 어느 쪽이 덜 흔들리는지 견준다.
    meds, means = [], []
    for _ in range(200_000):
        sample = np.random.choice(balls, size=5, replace=False)
        meds.append(np.median(sample))
        means.append(sample.mean())
    var_pop = (33 ** 2 - 1) / 12                            # 0..32 의 모분산
    sd_mean = (var_pop / 5 * (33 - 5) / (33 - 1)) ** 0.5    # 유한모집단 수정 포함
    print(f"20만회 표본중앙값의 평균   = {np.mean(meds):.6f}")
    print(f"20만회 표본중앙값의 표준편차 = {np.std(meds, ddof=1):.4f}  (정확값 {exact_sd:.4f})")
    print(f"20만회 표본평균의 표준편차   = {np.std(means, ddof=1):.4f}  (정확값 {sd_mean:.4f})")

    # 값마다 몇 번 나왔는지 센다. 표본이 50개뿐이라 히스토그램보다
    # 점그림이 낫다(2.5절에서 본 대로 자료가 적을 때의 선택이다).
    data_dict = {}
    for num in data:
        data_dict[num] = data_dict.get(num, 0) + 1

    fig, ax = plt.subplots(figsize=(12, 3))
    # 같은 값을 세로로 쌓아 점그림을 만든다
    for num, freq in data_dict.items():
        ax.plot([num] * freq, range(1, freq + 1), 'ok')
    # 참 중앙값 16. 점들이 이 선 주위에 흩어지는지 확인한다.
    ax.plot([16, 16], [0, 5], "--r", alpha=0.3, label="True median")
    ax.legend()
    ax.set_title('Simulation-Based Distribution of Sample Median')
    ax.set_xlabel('Sample Median')
    ax.set_ylabel('Number of Samples')
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_position("zero")
    plt.show()

if __name__ == "__main__":
    main()

출력:

모집단 중앙값 = 16.0
정확한 E[Med] = 16.000000,  sd(Med) = 5.830952
50회 표본중앙값의 평균 = 16.44
  그 평균의 표준오차 = 0.8246,  z = +0.534
20만회 표본중앙값의 평균   = 16.004935
20만회 표본중앙값의 표준편차 = 5.8225  (정확값 5.8310)
20만회 표본평균의 표준편차   = 3.9830  (정확값 3.9833)

완전열거가 준 \(E[\text{Med}] = 16.000000\)과 \(\operatorname{sd}(\text{Med}) = 5.830952 = \sqrt{34}\)가 (1)의 유도와 정확히 맞는다. 20만 번 모의실험의 표준편차 \(5.8225\)와 \(3.9830\)도 정확값 \(5.8310\), \(3.9833\)과 맞는다.

16.44는 편향의 증거가 아니다. 되풀이 50번으로 얻은 평균의 표준오차가 \(\sqrt{34}/\sqrt{50} = 0.8246\)이므로 16.44는 16에서 겨우 \(0.53\) 표준오차 떨어져 있다. 50번을 20만 번으로 늘리면 \(16.0049\)로 내려앉는데, 이때의 표준오차는 \(\sqrt{34}/\sqrt{200000} = 0.0130\)이니 여전히 참값과 \(0.4\) 표준오차 안이다. 되풀이 횟수를 늘려 줄어드는 것은 편향이 아니라 몬테카를로 오차다. 편향이었다면 20만 번을 돌려도 그 자리에 남아 있었을 것이다.

그림에서 한 가지 더 눈에 띄는 것이 있다. 점들이 가로축의 아무 데나 찍히지 않고 정수 자리에만 찍힌다. 공에 정수만 적혀 있고 다섯 개 중 가운데 값을 고르므로 표본중앙값도 정수일 수밖에 없다. 통계량의 분포는 이렇게 모집단의 성격을 물려받는다.

(1) 해석적으로. \(0 < p < 1\)에서 미분하면

이고, \(\ell'(p) = 0\)은 \(k(1-p) = (m-k)p\), 곧 \(k = mp\)를 준다. 따라서

이다. 이 정류점이 최대임은 이계도함수가 보여 준다.

\(\ell\)이 \((0,1)\)에서 위로 오목하므로 정류점은 하나뿐이고 그것이 최대다. 끝점도 따로 볼 필요가 없다. \(k = 0\)이면 \(\ell(p) = m\log(1-p)\)가 감소함수라 \(\hat p = 0\), \(k = m\)이면 증가함수라 \(\hat p = 1\)이어서 \(\hat p = k/m\)이 그대로 성립한다.

최대가능도추정량이 표본비율이다. 너무 당연해 보이는 답이지만, 그 당연함이 원리 하나에서 유도되었다는 것이 요점이다.

(2) 수치적으로. \(p\)의 후보를 격자에 늘어놓고 각각의 로그가능도를 재어 가장 큰 곳을 찾는다.

import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

np.random.seed(1)
p_true, m = 0.7, 100

# 참 p = 0.7 인 동전을 100번 던진다. 물론 실제로는 이 값을 모른다.
coins = np.random.binomial(n=1, p=p_true, size=m)
k = coins.sum()

# 로그가능도는 k log p + (m-k) log(1-p) 로 줄어든다. 확률을 100번 곱하면
# 0 으로 언더플로되므로 로그를 취해 합으로 바꾼다.
ps = np.linspace(0.01, 0.99, 100)
loglik = k * np.log(ps) + (m - k) * np.log(1 - ps)

grid_mle = ps[np.argmax(loglik)]   # 격자가 고른 최대점
exact_mle = k / m                  # (1) 에서 해석적으로 구한 값

print(f"앞면 k = {k},  m = {m}")
print(f"해석적  p-hat = k/m = {exact_mle:.4f}")
print(f"격자 최대점         = {grid_mle:.4f}  (격자 간격 {ps[1] - ps[0]:.4f})")

fig, (ax, az) = plt.subplots(1, 2, figsize=(12, 3.2),
                             gridspec_kw={"width_ratios": [1.6, 1]})

for a in (ax, az):
    a.plot(ps, loglik, color="#1565C0", lw=1.8, zorder=1)
    a.axvline(exact_mle, color="#D32F2F", ls="--", lw=1.6, zorder=3)
    a.axvline(grid_mle, color="#E65100", ls="-", lw=1.6, zorder=2)
    a.set_xlabel("$p$")
    a.spines[["top", "right"]].set_visible(False)

ax.set_ylabel("로그가능도 $\\ell(p)$")
ax.set_title("로그가능도 곡선 전체", fontsize=11)

# 오른쪽: 봉우리를 확대해 격자가 띄엄띄엄하다는 것을 보인다.
az.plot(ps, loglik, "o", ms=5, color="#1565C0", zorder=1)
az.set_xlim(exact_mle - 0.035, exact_mle + 0.035)
sel = np.abs(ps - exact_mle) < 0.045
az.set_ylim(loglik[sel].min(), loglik[sel].max() + 0.05)
az.set_title("봉우리 확대 — 격자점은 띄엄띄엄하다", fontsize=11)
az.plot([], [], color="#D32F2F", ls="--", lw=1.6,
        label=f"해석적 $\\hat p = k/m = {exact_mle:.4f}$")
az.plot([], [], color="#E65100", ls="-", lw=1.6,
        label=f"격자 최대점 $= {grid_mle:.4f}$")
az.legend(loc="lower center", fontsize=9)

plt.tight_layout()
plt.show()

출력:

앞면 k = 74,  m = 100
해석적  p-hat = k/m = 0.7400
격자 최대점         = 0.7425  (격자 간격 0.0099)

왼쪽 곡선의 봉우리가 해석적 답과 같은 자리에 있다. 오른쪽은 그 봉우리를 확대한 것인데, 격자가 고른 \(0.7425\)가 해석적 답 \(0.7400\)에서 한 칸 비껴나 있다. 격자 간격이 \(0.0099\)여서 \(0.74\)가 후보에 아예 없기 때문이다. 해석적으로 푼 답은 정확하고, 격자 탐색은 격자만큼만 정확하다.

모수가 하나뿐인 이 문제에서는 격자를 쓸 이유가 없다. 그러나 미분해서 손으로 풀 수 없는 모형에서는 이런 수치 탐색이 유일한 길이 되며, 그때는 격자 대신 기울기를 따라 올라가는 방법을 쓴다. 6장의 최대가능도 최적화가 그 이야기다.

(1) 해석적으로. \(N\)이 정수라 미분할 수 없으므로 이웃한 두 값의 비를 본다. \(L(N) = \binom{M}{m}\binom{N-M}{n-m}\big/\binom{N}{n}\)에서 \(\binom{M}{m}\)은 \(N\)과 무관하므로 약분되고

이다. \(\binom{a}{k}\big/\binom{a-1}{k} = a/(a-k)\)를 두 번 쓰면

을 얻는다. 이 비가 \(1\) 이상인 조건을 정리한다. 분모가 양수인 범위에서

인데, 양변을 펼치면 \(N^2\), \(-MN\), \(-Nn\) 항이 모두 지워지고

만 남는다. 즉 \(L\)은 \(N \le Mn/m\)까지 오르다가 그 뒤로는 내려간다. 따라서

이다. 비례식이 주는 \(Mn/m\)에 바닥함수가 붙는다는 것이 요점이고, 그 까닭은 \(N\)이 정수라는 것 하나다.

\(Mn/m\)이 정수일 때는 특별한 일이 일어난다. \(N = Mn/m\)에서 위 부등식이 등호가 되므로 비가 정확히 \(1\), 곧

이다. 최대가능도추정값이 둘이고 가능도만으로는 어느 쪽도 고를 수 없다.

(2) 수치적으로. \(M = 50\), \(n = 40\), \(m = 10\)은 \(Mn/m = 200\)으로 나누어떨어지는 바로 그 경우다. 그러므로 \(L(199) = L(200)\)이고, 코드가 199를 찍는 것은 격자가 성겨서가 아니라 동점이기 때문이다.

import matplotlib.pyplot as plt
from fractions import Fraction
from math import comb
from scipy import special

def prob(n, c, r, t):
    """포획-재포획의 초기하확률.

    n: 전체 개체수(우리가 추정하려는 미지수)
    c: 1차에서 잡아 표시한 수
    r: 2차에서 잡은 수
    t: 2차에서 잡힌 것 중 표시가 있던 수

    2차 표본 r마리를 고르는 모든 방법 중,
    표시된 것 t마리와 안 된 것 r-t마리를 고르는 방법의 비율이다.
    """
    return special.comb(n - c, r - t) * special.comb(c, t) / special.comb(n, r)

def exact_prob(n, c, r, t):
    """같은 확률을 분수로 계산한다. 반올림이 끼어들지 않는다."""
    return Fraction(comb(n - c, r - t) * comb(c, t), comb(n, r))

def capture_recapture(c=50, r=40, t=10):
    # 가능한 최소 개체수. 표시된 50마리와 2차에서 새로 잡힌 30마리는
    # 서로 다른 개체이므로 최소 50 + 40 - 10 = 80마리는 있어야 한다.
    min_n = c + r - t
    ns = range(min_n, 10 * min_n)

    # n을 바꿔 가며 가능도를 계산한다.
    # **n은 모수이지 확률변수가 아니다.** 자료 (c, r, t)는 고정해 두고
    # "어떤 n이 이 자료를 가장 그럴듯하게 만드는가"를 묻는 것이다.
    probs = [prob(n, c, r, t) for n in ns]

    mle_idx = probs.index(max(probs))
    mle_n = mle_idx + min_n
    # (1) 에서 유도한 floor(c*r/t) 와 비교해 보라.
    print(f"c={c}, r={r}, t={t}:  격자 최대점 = {mle_n},  floor(c*r/t) = {c * r // t}")
    return list(ns), probs, mle_n

ns, probs, mle_n = capture_recapture()

# 199 와 200 의 가능도는 **정확히 같다.** 분수로 재면 드러난다.
print(f"정확한 L(199) == L(200) ? {exact_prob(199, 50, 40, 10) == exact_prob(200, 50, 40, 10)}")
print(f"부동소수점 L(199) - L(200) = {prob(199, 50, 40, 10) - prob(200, 50, 40, 10):+.3e}")

# t를 11로 바꾸면 c*r/t = 181.8... 이라 나누어떨어지지 않고 최대점이 하나뿐이다.
capture_recapture(c=50, r=40, t=11)

# 봉우리가 얼마나 평평한지 재 본다.
for n in (150, 200, 300, 400):
    print(f"L({n}) / L(200) = {prob(n, 50, 40, 10) / prob(200, 50, 40, 10):.3f}")

fig, ax = plt.subplots(figsize=(12, 3))
ax.plot(ns, probs, label='Likelihood')
ax.axvline(mle_n, color='r', linestyle='--', label=f'MLE: N = {mle_n}')
ax.set_xlabel('Population Size (N)')
ax.set_ylabel('Probability')
ax.set_title('Capture–Recapture: Likelihood vs Population Size')
ax.legend()
plt.show()

출력:

c=50, r=40, t=10:  격자 최대점 = 199,  floor(c*r/t) = 200
정확한 L(199) == L(200) ? True
부동소수점 L(199) - L(200) = +2.665e-15
c=50, r=40, t=11:  격자 최대점 = 181,  floor(c*r/t) = 181
L(150) / L(200) = 0.424
L(200) / L(200) = 1.000
L(300) / L(200) = 0.347
L(400) / L(200) = 0.071

199가 나온 까닭이 드러났다. 분수로 재면 \(L(199)\)와 \(L(200)\)이 한 치도 다르지 않다. 그런데 부동소수점으로는 \(199\) 쪽이 \(2.7 \times 10^{-15}\)만큼 크게 계산되었고, 설령 비트까지 같았더라도 max()와 index()는 앞쪽 값을 고르므로 결과는 역시 199였을 것이다. 이 어긋남은 수치의 오차가 아니라 동점이며, 동점이 생긴다는 사실 자체가 (1)의 유도가 예측한 바다.

\(m\)을 11로 바꾸면 \(Mn/m = 181.82\)라 나누어떨어지지 않고, 격자가 찾은 답 181이 \(\lfloor 181.82 \rfloor = 181\)과 정확히 맞는다. 바닥함수가 제 몫을 하는 경우다.

곡선이 봉우리 주위에서 아주 평평하다는 점도 중요하다. \(N\)을 150으로 줄이거나 300으로 늘려도 가능도가 최댓값의 \(0.42\)배, \(0.35\)배로 남는다. 표시된 물고기 10마리라는 적은 정보로는 개체수를 정밀하게 못 맞힌다는 뜻이고, 실제 생태 조사에서 재포획 수를 늘리려 애쓰는 이유다.

(a) \(P(m \mid N) = \binom{M}{m} \binom{N - M}{n - m} / \binom{N}{n}\)이며 초기하분포이다. 주어진 값을 넣으면 \(P(N) = \binom{50}{10}\binom{N-50}{30}/\binom{N}{40}\).

(b) \(N\)은 정수 모수라 미분할 수 없다. 가능도비 \(L(N)/L(N-1)\)을 따지면 \(N \le Mn/m\)까지 오르다가 내려가므로 \(\hat N = \lfloor Mn/m \rfloor\)이다(보기 3). 여기서는 \(Mn/m = 50 \cdot 40 / 10 = 200\)이 정수라 \(L(199) = L(200)\)인 동점이고, 둘 다 최대가능도추정값이다. 이 추정량은 "(표시된 수) × (잡은 총 수) / (잡힌 표시된 수)"라는 비례 추론에 해당한다.

포획–재포획은 야생동물 개체수 추정의 기본 방법이다. 변형(폐쇄/개방 모집단, 여러 번의 재포획, 표지 손실)을 통해 풍부한 추정량 계열이 만들어진다.

(a) \(L(p) = \binom{100}{40} p^{40}(1-p)^{60} \propto p^{40}(1-p)^{60}\) (상수 배수는 최대화에 영향을 주지 않는다).

(b) 로그가능도: \(\ell(p) = 40\ln p + 60\ln(1-p)\). 미분하면 \(40/p - 60/(1-p) = 0 \Rightarrow p = 0.4\).

\(\hat p_{\text{MLE}} = 0.4\)로 표본비율과 같다. 일반적으로 \(X \sim \mathrm{Binomial}(n, p)\)에 대해 MLE는 \(\hat p = X/n\)이다.

로그가능도: \(\ell(\mu, \sigma^2) = -(n/2)\ln(2\pi\sigma^2) - (1/(2\sigma^2))\sum(X_i - \mu)^2\).

\(\partial \ell/\partial \mu = (1/\sigma^2)\sum(X_i - \mu) = 0 \Rightarrow \hat\mu = \bar X\).

\(\partial \ell/\partial \sigma^2 = -n/(2\sigma^2) + (1/(2\sigma^4))\sum(X_i - \hat\mu)^2 = 0 \Rightarrow \hat\sigma^2 = (1/n)\sum(X_i - \bar X)^2\).

참고: MLE는 \(n - 1\)이 아니라 \(n\)으로 나눈다. 따라서 \(\hat\sigma^2_{\text{MLE}}\)는 편향되어 있다: \(\mathbb{E}[\hat\sigma^2] = ((n-1)/n)\sigma^2\). 불편추정을 하려면 \(s^2 = \sum(X_i - \bar X)^2/(n - 1)\)을 사용한다(Bessel 수정).

MLE와 불편추정량의 구별은 반복해서 나타나는 주제이다. MLE는 점근적으로 최적이지만 유한표본에서는 편향될 수 있다.

Fisher-Neyman 인수분해 정리: \(T(X)\)가 \(\theta\)에 대해 충분일 필요충분조건은 결합밀도가 다음과 같이 인수분해되는 것이다.

여기서 \(g\)는 \(T(x)\)를 통해서만 \(\theta\)에 의존하고 \(h\)는 \(\theta\)에 의존하지 않는다.

포아송의 경우: \(f(x_1, \ldots, x_n \mid \lambda) = \prod_i \frac{e^{-\lambda} \lambda^{x_i}}{x_i!} = e^{-n\lambda} \lambda^{\sum x_i} / \prod_i x_i!\).

가능도가 \(g(\sum x_i, \lambda) \cdot h(x) = (e^{-n\lambda} \lambda^{\sum x_i}) \cdot (1/\prod x_i!)\)로 인수분해된다. 따라서 \(T(X) = \sum X_i\)는 충분통계량이다.

의의: \(\lambda\)에 관해 \(X_1, \ldots, X_n\)이 담고 있는 정보가 모두 \(\sum X_i\)에 집약되어 있다. MLE는 자료에 오직 \(T\)를 통해서만 의존하며, 추론에 자료 전체가 필요하지 않다. 이것이 Rao-Blackwell 정리를 통한 효율적 추정의 토대이다.

점수함수: \(\partial \log f/\partial \mu = (x - \mu)/\sigma^2\).

Fisher 정보량: \(I(\mu) = \mathbb{E}\!\left[(\partial \log f / \partial \mu)^2\right] = \mathbb{E}[(X - \mu)^2/\sigma^4] = \sigma^2/\sigma^4 = 1/\sigma^2\).

크기 \(n\)인 i.i.d. 표본에 대해서는 \(I_n(\mu) = n/\sigma^2\)이다.

왜 중요한가: Cramér-Rao 하한. \(\mu\)의 임의의 불편추정량의 분산은 적어도 \(1/I_n(\mu) = \sigma^2/n\)이다. \(\mathrm{Var}(\bar X) = \sigma^2/n\)이 정확히 성립하므로 \(\bar X\)는 Cramér-Rao 하한을 달성하며 효율적이다. 어떤 불편추정량도 이보다 나을 수 없다.

Fisher 정보량은 모수에 관해 표본이 담은 "정보량"을 정량화하고 추정량 분산의 하한을 준다. MLE의 점근이론, 실험설계, 정보기하학에서 쓰인다.

Poisson(\(\lambda\))에서 \(\hat\lambda_{\text{MLE}} = \bar X\)이다. 관측값 하나당 Fisher 정보량은 \(I(\lambda) = 1/\lambda\)이다(\(\partial \log f/\partial \lambda = X/\lambda - 1\)이고 \(\mathbb{E}[(X/\lambda - 1)^2] = \mathrm{Var}(X)/\lambda^2 = 1/\lambda\)이므로).

점근분포: \(\sqrt n(\bar X - \lambda) \xrightarrow{d} N(0, \lambda)\)이며, 이는 \(1/I(\lambda) = \lambda\)와 일치한다.

중심극한정리로 직접 확인: \(\mathrm{Var}(X_i) = \lambda\)인 \(\bar X = (1/n)\sum X_i\)이므로 중심극한정리에 의해 \(\sqrt n(\bar X - \lambda) \to N(0, \lambda)\). ✓

일반적 의의: MLE는 점근적으로 정규분포를 따르며 그 분산은 관측값 하나당 Fisher 정보량의 역수이다. 이로부터 다음을 얻는다:

• 점근 신뢰구간: \(\hat\theta \pm 1.96/\sqrt{n I(\hat\theta)}\) (\(I(\theta)\)의 대입추정값으로 \(I(\hat\theta)\)를 사용).
• 점근 효율성: MLE는 점근적으로 Cramér-Rao 하한을 달성한다.

이 결과들이 현대 통계학에서 가능도 기반 추론이 중심적 위치를 차지하는 이유이다.

감마분포의 적률은 \(E[X] = k\theta\), \(\operatorname{Var}(X) = k\theta^2\)이다. 표본적률과 맞추면

이고, 둘째 식을 첫째 식으로 나누면

이다.

장점. 계산이 한 줄이다. 감마분포의 최대가능도방정식은

꼴로 디감마함수가 들어 있어 수치적으로 풀어야 한다. 적률법 추정값은 그 반복의 좋은 출발점이 된다. 또 분포 전체를 가정하지 않고 적률만 맞추므로 모형이 조금 어긋나도 크게 망가지지 않는다.

단점. 일반적으로 최대가능도추정보다 비효율적이다. 자료의 정보를 처음 몇 개의 적률로만 요약하기 때문이다. 고차 적률을 쓰면 표집변동이 커져 더 불안정해지고, 추정값이 모수공간 밖으로 나가는 일도 있다(예: 분산 추정이 음수). 또 충분통계량을 쓰지 않으므로 정보 손실이 생긴다.

실무에서는 적률법을 초기값이나 빠른 점검으로 쓰고 최종 추정은 최대가능도로 하는 조합이 흔하다.

\(W = \sum_i(X_i-\bar X)^2\)로 두면 \(W/\sigma^2 \sim \chi^2_{n-1}\)이므로

이다. 따라서

이다. \(c\)로 미분해 0으로 두면

을 얻는다.

세 추정량. MSE를 \(\sigma^4\) 단위로 적으면

\(n=10\)이면 각각 \(0.2222\), \(0.1900\), \(0.1818\)로, 불편추정량이 셋 중 가장 나쁘다.

뜻. 불편성은 좋은 성질이지만 최적성의 기준은 아니다. 편향을 조금 받아들이고 분산을 더 줄이면 전체 오차가 작아질 수 있으며, 이것이 편향-분산 맞바꿈이다. 능형회귀, 라소, 축소추정, 정칙화가 모두 같은 거래를 한다.

그런데도 실무에서 \(n-1\)을 쓰는 이유가 있다. 첫째, 불편성이 여러 표본을 결합할 때 좋은 성질을 준다(분산분석에서 제곱평균들을 더하고 나눌 때 편향이 누적되지 않는다). 둘째, 최소 MSE인 \(c=1/(n+1)\)은 모집단이 정규일 때만 최적이라 일반성이 없다. 셋째, \(n\)이 크면 셋의 차이가 \(O(1/n^2)\)로 사라진다.

불변성. \(g\)가 일대일이면 모수를 \(\eta = g(\theta)\)로 바꿔 쓴 가능도가 \(L^*(\eta) = L(g^{-1}(\eta))\)이므로, \(L\)을 최대로 하는 \(\hat\theta\)에서 \(L^*\)가 최대가 되고 그 위치가 \(\hat\eta = g(\hat\theta)\)이다. \(g\)가 일대일이 아니면 유도가능도 \(L^*(\eta) = \sup_{\theta:\,g(\theta)=\eta}L(\theta)\)로 정의하면 같은 결론이 나온다.

오즈의 MLE. \(\hat p = \bar X\)이므로

이다. 100번 중 40번 성공이면 \(0.4/0.6 = 2/3\)이다. 따로 최적화할 필요가 없다.

불편성에는 없다. \(E[\hat\theta] = \theta\)라 해도 일반적으로 \(E[g(\hat\theta)] \ne g(\theta)\)이다. 옌센 부등식에 따라 \(g\)가 볼록이면 \(E[g(\hat\theta)] \ge g(E[\hat\theta]) = g(\theta)\)로 위로 치우친다.

위 예에서도 \(x/(1-x)\)가 \((0,1)\)에서 볼록이므로 오즈의 MLE는 오즈를 과대추정한다. 게다가 \(\bar X = 1\)이면 값이 무한대가 되어 버린다. 그래서 로그오즈를 다룰 때는 \(\bar X\) 대신 \((k+0.5)/(n+1)\)처럼 살짝 보정한 값(하딘-피터스 보정)을 쓰기도 한다.

이 대비가 두 성질의 성격을 잘 보여 준다. 불변성은 "어떤 모수화로 문제를 적든 답이 같다"는 뜻이고, 불편성은 특정 모수화에 묶여 있다. 표준편차의 불편추정량이 분산의 불편추정량의 제곱근이 아니라는 사실도 같은 이야기다.

(가) 불편이지만 일치가 아닌 경우. \(X_1,\dots,X_n \sim N(\mu,\sigma^2)\)에서 \(\hat\mu = X_1\)(첫 관측값만 쓴다)을 생각하자. \(E[X_1] = \mu\)로 불편이지만, \(n\)이 아무리 커져도 분포가 \(N(\mu,\sigma^2)\) 그대로라 \(\mu\)로 수렴하지 않는다. 일치가 아니다.

자료를 버리는 극단적인 예지만 요점은 분명하다. 불편성은 표본크기가 커지는 것과 아무 상관이 없는, 한 표본크기에서의 성질이다.

(나) 편향되었지만 일치인 경우. 같은 설정에서 최대가능도 분산추정량

은 \(E[\hat\sigma^2_{\text{MLE}}] = \frac{n-1}{n}\sigma^2\)로 편향되어 있다. 그러나 편향이 \(-\sigma^2/n \to 0\)이고 분산도 0으로 가므로 \(\hat\sigma^2_{\text{MLE}} \xrightarrow{p} \sigma^2\)이다. 일치추정량이다.

다른 예로 포획-재포획 추정량 \(\hat N = Mn/m\)은 일치이지만 기댓값이 아예 존재하지 않는다. \(P(m = 0) > 0\)이라 \(Mn/m\)이 정의되지 않는 표본이 양의 확률로 나오기 때문이다(\(M = 50\), \(n = 40\), \(N = 200\)에서 \(P(m=0) \approx 2.2 \times 10^{-6}\)). 편향을 따지기 전에 기댓값부터 없는 셈이며, 분모에 \(1\)을 더한 채프먼 추정량 \(\hat N_C = (M+1)(n+1)/(m+1) - 1\)을 쓰는 까닭이 그것이다.

정리. 두 성질은 서로 독립적이다.

• 불편성: \(E[\hat\theta_n] = \theta\). 고정된 \(n\)에서의 성질이며, "평균적으로 맞다"는 뜻이다.
• 일치성: \(\hat\theta_n \xrightarrow{p} \theta\). \(n\to\infty\)에서의 성질이며, "자료를 모으면 결국 맞는다"는 뜻이다.

실무에서 더 중요한 쪽은 일치성이다. 일치가 아닌 추정량은 자료를 아무리 모아도 참값에 다가가지 않으므로 쓸 수 없다. 편향은 크기가 작고 \(n\)과 함께 사라지면 대개 감수할 만하며, 연습문제 8에서 보았듯 일부러 편향을 들여 오차를 줄이기도 한다.

충분조건 하나를 기억해 두면 편하다. 편향과 분산이 모두 0으로 가면 일치이다(MSE 수렴이 확률수렴을 함의하므로).