비율 신뢰구간의 실제 포함률¶

(1) 왈드의 천장. \(\hat p = 0\)이면 추정 표준오차가 \(\sqrt{0 \times 1/n} = 0\)이라 구간이 \([0, 0]\) 한 점으로 붕괴한다. \(p = 0.01 > 0\)이므로 그런 표본에서는 절대 참값을 담지 못한다. 그 일이 일어날 확률이

이므로 포함률은

을 넘을 수 없다. 표본의 \(67\%\)를 버리고 시작하는 구간이다.

윌슨은 붕괴하지 않는다. \(\hat p = 0\)을 윌슨의 부등식에 넣으면

이므로 구간이 \(\left[0,\ \frac{1.96^2}{40+1.96^2}\right] = [0,\ 0.08762]\)다. 폭이 남아 있고 \(p = 0.01\)을 담는다. 분모에 \(\hat p\)이 아니라 \(p\)를 둔 것 하나가 이 차이를 만든다.

(2) 모의실험.

import matplotlib.pyplot as plt
import numpy as np

rng = np.random.default_rng(1)
z, REP = 1.96, 100_000


def coverage(n, p):
    """명목 95% Wald 구간과 Wilson 구간의 실제 포함률을 모의실험으로 구한다."""
    phat = rng.binomial(n, p, size=REP) / n

    # Wald: 표준오차에 p 대신 p^ 를 그대로 꽂는다. p^ = 0 이면 구간이 [0, 0].
    half = z * np.sqrt(phat * (1 - phat) / n)
    wald = np.mean((phat - half <= p) & (p <= phat + half))

    # Wilson: |p - p^| <= z sqrt(p(1-p)/n) 을 p 에 대해 풀어 얻는다.
    c = 1 + z**2 / n
    center = (phat + z**2 / (2 * n)) / c
    half_w = z * np.sqrt(phat * (1 - phat) / n + z**2 / (4 * n**2)) / c
    wilson = np.mean((center - half_w <= p) & (p <= center + half_w))

    return wald, wilson


fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

# (좌) n 을 40 으로 고정하고 참값 p 를 훑는다.
ps = np.arange(0.01, 0.5001, 0.005)
cov = np.array([coverage(40, q) for q in ps])
ax1.plot(ps, cov[:, 0], "-", color="tab:red", lw=1.5, label="Wald")
ax1.plot(ps, cov[:, 1], "-", color="tab:blue", lw=1.5, label="Wilson")
ax1.axhline(0.95, color="black", ls="--", lw=1, label="nominal 0.95")
ax1.set_xlabel("true $p$")
ax1.set_ylabel("actual coverage")
ax1.set_title("n = 40 fixed")
ax1.set_ylim(0.25, 1.0)
ax1.legend(fontsize=9, loc="lower right")

# (우) p 를 0.1 로 고정하고 표본크기를 키운다. 톱니 진동을 보려면 n 을 1 씩 늘려야 한다.
ns = np.arange(10, 201)
cov2 = np.array([coverage(int(m), 0.1) for m in ns])
ax2.plot(ns, cov2[:, 0], "-", color="tab:red", lw=1.2, label="Wald")
ax2.plot(ns, cov2[:, 1], "-", color="tab:blue", lw=1.2, label="Wilson")
ax2.axhline(0.95, color="black", ls="--", lw=1, label="nominal 0.95")
ax2.set_xlabel("sample size $n$")
ax2.set_ylabel("actual coverage")
ax2.set_title("p = 0.1 fixed")
ax2.set_ylim(0.6, 1.0)
ax2.legend(fontsize=9, loc="lower right")

plt.tight_layout()
plt.show()

# 두 패널을 수로 요약한다.
print(f"n = 40, p = 0.01:  P(p-hat = 0) = {0.99 ** 40:.4f}  ->  Wald 포함률의 상한 {1 - 0.99 ** 40:.4f},  모의 {cov[0, 0]:.4f}")
print(f"  p-hat = 0 일 때 Wilson 구간 = [0, {z**2 / (40 + z**2):.4f}]  (p = 0.01 을 담는다)")
print(f"왼쪽 패널 {len(ps)} 개 p 에 대해   Wald  최소 {cov[:, 0].min():.4f}  평균 {cov[:, 0].mean():.4f}  0.95 이상인 비율 {(cov[:, 0] >= 0.95).mean():.3f}")
print(f"                              Wilson 최소 {cov[:, 1].min():.4f}  평균 {cov[:, 1].mean():.4f}  0.95 이상인 비율 {(cov[:, 1] >= 0.95).mean():.3f}")
i29, i30 = list(ns).index(29), list(ns).index(30)
print(f"오른쪽 패널  n = 29 -> 30:  Wald {cov2[i29, 0]:.4f} -> {cov2[i30, 0]:.4f}  (낙차 {cov2[i29, 0] - cov2[i30, 0]:.4f})")
print(f"             n >= 100 에서  Wald  최소 {cov2[90:, 0].min():.4f}  최대 {cov2[90:, 0].max():.4f}")
print(f"                            Wilson 최소 {cov2[90:, 1].min():.4f}  최대 {cov2[90:, 1].max():.4f}")

출력:

n = 40, p = 0.01:  P(p-hat = 0) = 0.6690  ->  Wald 포함률의 상한 0.3310,  모의 0.3313
  p-hat = 0 일 때 Wilson 구간 = [0, 0.0876]  (p = 0.01 을 담는다)
왼쪽 패널 99 개 p 에 대해   Wald  최소 0.3313  평균 0.9017  0.95 이상인 비율 0.030
                              Wilson 최소 0.9216  평균 0.9522  0.95 이상인 비율 0.545
오른쪽 패널  n = 29 -> 30:  Wald 0.9465 -> 0.8093  (낙차 0.1371)
             n >= 100 에서  Wald  최소 0.9107  최대 0.9538
                            Wilson 최소 0.9352  최대 0.9651

천장이 그대로 달성된다. 모의 포함률 \(0.3313\)이 (1)에서 구한 상한 \(0.33103\)과 비율의 몬테카를로 오차 \(\sqrt{0.33 \times 0.67/10^5} = 0.0015\) 안에서 같다. 상한이 "넘을 수 없다"가 아니라 "정확히 그만큼"이라는 뜻이고, \(\hat p \ge 1/40\)인 표본에서는 왈드 구간이 전부 \(p = 0.01\)을 담는다는 것이다. 왈드가 놓치는 표본과 \(\hat p = 0\)인 표본이 같은 집합이다.

왼쪽 패널의 절벽. \(p\)가 작은 쪽에서 빨간 선이 수직으로 떨어지는 것은 \((1-p)^{40}\)이 급격히 커지기 때문이다. 같은 이유로 \(p\)가 조금만 커져도 절벽이 사라진다. \(p = 0.05\)에서 \((1-p)^{40} = 0.1285\), \(p = 0.1\)에서 \(0.0148\)이다.

그런데 절벽을 지나서도 왈드는 명목 수준에 닿지 못한다. 훑어본 \(99\)개 \(p\) 가운데 포함률이 \(0.95\) 이상인 것이 \(3\%\)뿐이고 평균이 \(0.9017\)이다. 윌슨은 최소가 \(0.9216\), 평균이 \(0.9522\), \(0.95\) 이상인 비율이 \(54.5\%\)다. 왈드의 문제는 극단적인 \(p\)에만 있는 것이 아니라 전 구간에 있다.

오른쪽 패널의 톱니. \(n = 29\)에서 \(0.9465\)였던 왈드 포함률이 \(n = 30\)에서 \(0.8093\)으로 \(0.1371\) 떨어진다. 표본을 하나 더 얻고 포함률 \(14\) 퍼센트포인트를 잃었다. 까닭은 이산성이다. \(k = 1\)의 왈드 상한이 \(n = 29\)에서 \(0.10089\)였다가 \(n = 30\)에서 \(0.09757\)로 내려가 \(p = 0.1\)을 넘지 못하게 되고, 그 순간 \(P(X = 1) = 0.1413\)이라는 점질량이 통째로 빠진다. 낙차 \(0.1371\)이 이 점질량과 거의 같다.

톱니는 잡음이 아니다. 되풀이를 10만 번 했으므로 한 점의 몬테카를로 오차가 \(0.0015\) 남짓인데 톱니의 골은 그보다 100배 깊다. 같은 계산을 모의실험 없이 정확히 해도 같은 톱니가 나오며, 그것을 다음 보기에서 확인한다. \(n \ge 100\)에서도 왈드가 \(0.9107\)에서 \(0.9538\) 사이를 오르내리고 윌슨은 \(0.9352\)에서 \(0.9651\) 사이다. 어느 쪽도 한 값으로 수렴하지 않는다. 이항 자료에서 포함률이 \(n\)의 매끄러운 함수가 될 수 없기 때문이다.

(1) 정확히 계산된다. \(\hat p\)이 가질 수 있는 값은 \(k/n\) 꼴의 \(n+1\)개뿐이고 각 값의 확률은 이항 PMF가 준다. 구간이 참값을 덮는지는 \(k\)만 알면 결정되므로, 포함률은 덮는 \(k\)들의 점질량을 더한 유한합이다.

여기서 \(I(\hat p)\)가 왈드든 윌슨이든 무엇이든 상관없다. 모의실험은 확인용이고, 참값은 이 합이다. 앞 보기의 톱니가 잡음이 아니라 구조라는 것도 이 식에서 바로 보인다. \(n\)이 하나 늘면 지시함수가 \(1\)에서 \(0\)으로 뒤집히는 \(k\)가 생기고, 그 \(k\)의 점질량이 합에서 통째로 사라진다.

되풀이 10만 번의 모의 포함률은 이 참값을 \(\sqrt{C(1-C)/10^5} \le 0.0016\)의 오차로 추정한다. 그러므로 소수 셋째 자리까지는 믿을 수 있고 넷째 자리는 믿을 수 없다.

(2) 두 값을 나란히.

import numpy as np
from scipy import stats

z, REP = 1.96, 100_000
rng = np.random.default_rng(1)


def coverage(n, p):
    """보기 1 과 같다. 모의실험으로 두 구간의 포함률을 구한다."""
    phat = rng.binomial(n, p, size=REP) / n
    half = z * np.sqrt(phat * (1 - phat) / n)
    wald = np.mean((phat - half <= p) & (p <= phat + half))
    c = 1 + z**2 / n
    center = (phat + z**2 / (2 * n)) / c
    half_w = z * np.sqrt(phat * (1 - phat) / n + z**2 / (4 * n**2)) / c
    wilson = np.mean((center - half_w <= p) & (p <= center + half_w))
    return wald, wilson


def exact_coverage(n, p):
    """포함률은 유한한 점질량의 합이므로 모의실험 없이 정확히 계산된다."""
    k = np.arange(n + 1)
    phat = k / n
    pmf = stats.binom(n, p).pmf(k)
    half = z * np.sqrt(phat * (1 - phat) / n)
    wald = pmf[(phat - half <= p) & (p <= phat + half)].sum()
    c = 1 + z**2 / n
    center = (phat + z**2 / (2 * n)) / c
    half_w = z * np.sqrt(phat * (1 - phat) / n + z**2 / (4 * n**2)) / c
    wilson = pmf[(center - half_w <= p) & (p <= center + half_w)].sum()
    return wald, wilson


print("n = 40 고정          Wald            Wilson")
print("    p      정확    모의     정확    모의    P(p^=0)")
for p in (0.01, 0.02, 0.05, 0.10, 0.20, 0.30, 0.50):
    w, s = coverage(40, p)
    we, se_ = exact_coverage(40, p)
    print(f"{p:>5.2f}   {we:.4f}  {w:.4f}   {se_:.4f}  {s:.4f}   {(1-p)**40:>7.4f}")

print()
print("p = 0.1 고정         Wald            Wilson")
print("   n      정확    모의     정확    모의")
for n in (10, 20, 29, 30, 40, 50, 100, 200):
    w, s = coverage(n, 0.1)
    we, se_ = exact_coverage(n, 0.1)
    print(f"{n:>4}   {we:.4f}  {w:.4f}   {se_:.4f}  {s:.4f}")

출력:

n = 40 고정          Wald            Wilson
    p      정확    모의     정확    모의    P(p^=0)
 0.01   0.3310  0.3313   0.9393  0.9393    0.6690
 0.02   0.5531  0.5525   0.9543  0.9549    0.4457
 0.05   0.8681  0.8670   0.9520  0.9528    0.1285
 0.10   0.9145  0.9133   0.9433  0.9432    0.0148
 0.20   0.9047  0.9036   0.9283  0.9275    0.0001
 0.30   0.9299  0.9300   0.9443  0.9442    0.0000
 0.50   0.9193  0.9191   0.9615  0.9614    0.0000

p = 0.1 고정         Wald            Wilson
   n      정확    모의     정확    모의
  10   0.6497  0.6496   0.9298  0.9294
  20   0.8760  0.8774   0.9568  0.9572
  29   0.9467  0.9466   0.9784  0.9783
  30   0.8085  0.8093   0.9742  0.9741
  40   0.9145  0.9134   0.9433  0.9431
  50   0.8789  0.8779   0.9703  0.9702
 100   0.9324  0.9325   0.9364  0.9375
 200   0.9271  0.9267   0.9561  0.9551

정확값과 모의값이 모두 맞는다. 서른 쌍 가운데 가장 크게 벌어진 것이 \(n = 20\)의 왈드(\(0.8760\) 대 \(0.8774\))와 \(n = 100\)의 윌슨(\(0.9364\) 대 \(0.9375\))으로 둘 다 \(0.0011\)–\(0.0014\)이고, 그 자리의 몬테카를로 오차 \(0.0010\)–\(0.0008\)의 \(1.3\)–\(1.4\)배다. 나머지는 소수 셋째 자리까지 그대로 겹친다. 앞 보기의 톱니는 모의실험이 만든 것이 아니다.

왈드가 놓치는 표본의 정체. 위쪽 표에서 정확한 포함률과 마지막 열을 더해 보라.

\(p = 0.01\)에서는 소수 넷째 자리까지 정확히 \(1\)이다. 왈드가 놓치는 표본이 곧 \(\hat p = 0\)인 표본이라는 것이며, 그 밖의 표본에서는 하나도 놓치지 않는다. \(p = 0.02\)에서 \(0.9988\)로 조금 모자라는 것은 \(k \ge 5\)인 표본 때문이다. \(k = 5\)면 \(\hat p = 0.125\)이고 구간이 \([0.0225,\ 0.2275]\)가 되어 하한이 \(p = 0.02\)를 넘어선다. \(k \ge 5\)의 점질량을 모두 더하면 \(0.0012\)로, 모자라는 양과 정확히 맞는다.

\(p\)가 극단이 아닐 때도 모자란다. \(p = 0.2\)에서 정확한 왈드 포함률이 \(0.9047\), \(p = 0.5\)에서 \(0.9193\)이다. \(\hat p = 0\)이 사실상 일어나지 않는 자리인데도 \(0.95\)에 못 미친다. 이것은 붕괴와 다른 원인, 곧 \(\hat p(1-\hat p)\)이 \(p(1-p)\)를 평균적으로 작게 추정한다는 사실에서 온다(\(E[\hat p(1-\hat p)] = \frac{n-1}{n}pq\)). 구간이 체계적으로 좁다.

아래쪽 표. \(n = 29 \to 30\)의 정확한 낙차가 \(0.9467 \to 0.8085\)로 \(0.1382\)이고, \(n = 30\)에서 \(P(X=1) = 0.1413\)인 점질량이 빠진 것과 맞는다. 그리고 \(n = 200\)에서도 왈드가 \(0.9271\)이다. \(n\)을 스무 배로 늘려도 \(0.95\)에 닿지 못한다. 윌슨은 여덟 자리 모두 \(0.93\) 위이고 평균적으로 \(0.95\) 근처다.

왈드 구간은

으로 한 점이다. \(p = 0.01 > 0\)이므로 이 구간은 참값을 담지 못한다.

\(\hat p = 0\)이 될 확률은

이므로 포함률은 나머지 표본에서 전부 성공한다 해도

을 넘지 못한다. 모의실험이 0.3313을 주었으므로 이 상한이 실제로 달성된다.

명목 95%라고 적힌 구간이 33%를 담는다. 신뢰구간이 틀렸다기보다, 근사의 전제가 무너진 자리에서 공식을 기계적으로 쓴 결과다.

윌슨 구간의 정의는 \(|\hat p - p| \le z\sqrt{p(1-p)/n}\)을 만족하는 \(p\)의 집합이다. \(\hat p = 0\)을 넣고 양변을 제곱하면

이다. \(p = 0\)은 항상 성립하므로 \(p > 0\)에서 양변을 \(p\)로 나눈다.

따라서 구간은 \(\left[0,\ \dfrac{z^2}{n+z^2}\right]\)이다. \(\square\)

상한이 대략 \(z^2/n \approx 3.84/n\)이다. 앞 페이지에서 본 삼의 법칙 \(3/n\)과 같은 차수이며, 윌슨 구간이 양측이라 조금 더 크다.

\(\hat p\)이 \(n+1\)개의 값만 가지므로 모든 경우를 열거하면 된다.

import numpy as np
from scipy import stats

def coverage_exact(n, p, z=1.96):
    k = np.arange(n + 1)
    phat = k / n
    pmf = stats.binom(n, p).pmf(k)          # 각 경우의 확률
    half = z * np.sqrt(phat * (1 - phat) / n)
    return pmf[(phat - half <= p) & (p <= phat + half)].sum()

\(n = 40\), \(p = 0.01\)에서 이 함수는 \(0.3310\)을 준다. 모의실험이 준 0.3313과 0.0003 차이이며, 10만 회 반복의 표준오차 \(\sqrt{0.33 \times 0.67/100000} = 0.0015\) 안에 있다.

정확 계산이 더 나은 점. 오차가 없고, 계산이 \(O(n)\)으로 빠르며, 그림의 톱니가 진짜 톱니인지 모의실험 잡음인지 구별해 준다.

그럼에도 모의실험이 쓸모 있는 이유. 구간이 \(\hat p\)의 닫힌 함수로 적히지 않으면 열거할 수가 없다. 붓스트랩 구간이나 베이즈 신용구간이 그렇고, 표본추출 구조가 복잡해지면(층화, 군집, 무응답) 더욱 그렇다. 게다가 "많이 반복해서 몇 번 담는지 센다"는 코드가 포함률의 정의 그 자체다.

실무 지침은 이렇다. 정확히 계산할 수 있으면 하고, 없으면 모의실험한다. 이 페이지에서 두 방법이 일치하는 것을 확인했으므로 뒤이어 나올 복잡한 구간에서도 모의실험을 믿을 수 있다.

\(k = 1\)이 탈락한다. 왈드 구간의 상한을 계산하면

이다. \(n = 29\)에서는 상한이 0.1을 아슬아슬하게 넘지만 \(n = 30\)에서는 넘지 못한다.

빠지는 확률의 크기. \(n = 30\), \(p = 0.1\)에서

이 통째로 포함률에서 빠진다. 실제 낙차는 \(0.9466 - 0.8093 = 0.1373\)으로 이 값에 가깝다(\(k = 0\)의 탈락 여부 등 다른 변화가 조금 섞여 있다).

일반적인 그림. \(k = 1\)의 상한은 \(n\)에 대해 \(\approx 2.96/n\)으로 연속적으로 감소하는데, 포함률은 이 값이 \(p\)를 지나는 순간 불연속적으로 떨어진다. \(k = 2, 3, \ldots\)도 각각 자기 지점에서 같은 일을 일으키므로 톱니가 여러 개 생긴다. \(n\)이 커지면 개별 점질량이 작아져 진폭이 줄지만 그 감소가 \(1/\sqrt n\) 속도로 느리다.

양변을 제곱하고 정리한다.

\(p\)에 대한 이차부등식으로 모은다.

\(a = 1 + z^2/n\), \(b = -(2\hat p + z^2/n)\), \(c = \hat p^2\)으로 두면 위로 열린 포물선이므로 두 근 사이가 해집합이다. 근의 공식에서

이고 근호 안을 전개하면

이 되어 근호가 \(z\sqrt{\dfrac{\hat p(1-\hat p)}{n} + \dfrac{z^2}{4n^2}}\)이다. 정리하면 개요의 공식이다. \(\square\)

구조를 읽는다. 중심 \(\dfrac{\hat p + z^2/(2n)}{1 + z^2/n}\)은 \(\hat p\)과 \(1/2\)의 가중평균이고, \(n\)이 크면 \(\hat p\)에 수렴한다. 반폭도 \(z\sqrt{\hat p(1-\hat p)/n}\)에 수렴하므로 \(n\)이 크면 윌슨과 왈드가 같아진다. 두 구간의 차이는 오직 작은 \(n\)과 극단적인 \(\hat p\)에서 나타나는데, 그것이 바로 왈드가 실패하는 자리다.

점근적으로는 참이다. \(p\)를 고정하고 \(n \to \infty\)를 보내면 \(\hat p \overset{p}{\to} p\)이고 \(\sqrt{\hat p(1-\hat p)} \to \sqrt{p(1-p)}\)이므로 슬러츠키 정리에 의해

이고 포함률이 0.95로 간다. 문제는 이 수렴이 \(p\)에 대해 균일하지 않다는 것이다.

반례의 구성. 주어진 \(n\)에 대해 \(p = \dfrac{1}{2n}\)으로 잡는다. 그러면

이므로 \(\hat p = 0\)일 확률이 \(n\)과 무관하게 60% 남는다. 왈드 구간은 그때 붕괴하므로 포함률이

를 넘지 못한다. 정확히 계산하면 다음과 같다.

\(n = 10{,}000\)에서도 0.393이다. \(\square\)

사실 \(p \to 0\)으로 더 보내면 \(\inf_p C(n,p) \to 0\)이다. 점별 수렴은 있으나 균일 수렴이 없다.

왜 이것이 실무적으로 중요한가. 우리는 \(p\)를 모른다. 신뢰구간을 쓰는 목적이 바로 \(p\)를 모르는 상태에서 판단하는 것이다. "\(p\)가 0.2와 0.8 사이면 왈드도 괜찮다"는 조건부 보장은, 그 조건을 확인할 수 없는 상황에서 쓸모가 없다.

대안의 비교. \(n = 40\)에서 몇 가지 방법의 정확한 포함률을 계산하면 이렇다.

\(p \in [0.01, 0.5]\)에서 평균 포함률은 왈드 0.9018, 윌슨 0.9522이고 최솟값은 왈드 0.3310, 윌슨 0.9221이다. 클로퍼-피어슨은 모든 \(p\)에서 0.95 이상을 보장하지만 그 대가로 구간이 넓다(보수적). 윌슨은 0.95를 조금 밑도는 곳이 있으나 평균적으로 가장 0.95에 가깝다. 왈드는 어느 기준에서도 선택할 이유가 없다. 계산이 쉽다는 것이 유일한 장점이고, 그 장점은 컴퓨터가 있는 시대에 장점이 아니다.

import numpy as np
from scipy import stats

z = 1.959964

def ints(k, n):
    ph = k / n
    wald = (ph - z*np.sqrt(ph*(1-ph)/n), ph + z*np.sqrt(ph*(1-ph)/n))
    ctr = (ph + z**2/(2*n)) / (1 + z**2/n)
    half = z*np.sqrt(ph*(1-ph)/n + z**2/(4*n**2)) / (1 + z**2/n)
    wil = (ctr - half, ctr + half)
    nt = n + 4; pt = (k + 2) / nt
    ac = (pt - z*np.sqrt(pt*(1-pt)/nt), pt + z*np.sqrt(pt*(1-pt)/nt))
    lo = 0 if k == 0 else stats.beta.ppf(0.025, k, n-k+1)
    hi = 1 if k == n else stats.beta.ppf(0.975, k+1, n-k)
    return wald, wil, ac, (lo, hi)

def cover_width(n, p, i):
    k = np.arange(n+1); pm = stats.binom.pmf(k, n, p)
    cov = w = 0.0
    for kk, pr in zip(k, pm):
        iv = ints(kk, n)[i]
        if iv[0] <= p <= iv[1]: cov += pr
        w += pr * (iv[1] - iv[0])
    return cov, w

ps = np.linspace(0.01, 0.99, 99)
print(f"{'방법':>9}{'평균 포함률':>13}{'최소 포함률':>13}{'평균 폭':>11}")
for i, name in enumerate(["Wald", "Wilson", "AC", "CP"]):
    res = [cover_width(40, p, i) for p in ps]
    cov = np.mean([r[0] for r in res]); mn = min(r[0] for r in res)
    wid = np.mean([r[1] for r in res])
    print(f"{name:>9}{cov:>13.4f}{mn:>13.4f}{wid:>11.4f}")

출력:

       방법       평균 포함률       최소 포함률       평균 폭
     Wald       0.8998       0.3310     0.2388
   Wilson       0.9508       0.9283     0.2382
       AC       0.9576       0.9382     0.2451
       CP       0.9706       0.9519     0.2601

윌슨이 왈드보다 오히려 좁다(\(0.2382\) 대 \(0.2388\)). 그러면서 평균 포함률은 \(0.90\)에서 \(0.95\)로 오르고 최소 포함률은 \(0.33\)에서 \(0.93\)으로 뛴다. 아무 대가 없이 개선된다.

이것이 왈드를 쓸 이유가 없다는 말의 가장 강한 형태다. 보통 "정확도를 높이려면 구간을 넓혀야 한다"고 생각하지만, 왈드는 넓어서 부정확한 것이 아니라 엉뚱한 자리에 놓여서 부정확하다. \(\hat p\)가 극단이면 \(\sqrt{\hat p(1-\hat p)}\)가 작아져 구간이 좁아지는데, 하필 그때 \(\hat p\)가 \(p\)에서 멀다.

통념이 맞는 것은 클로퍼–피어슨이다. 폭이 \(0.2601\)로 윌슨보다 \(9\%\) 넓고, 그 대가로 최소 포함률 \(0.95\) 이상을 보장한다. 여기서는 "보수적 = 넓다"가 성립한다.

선택 기준. 보장이 필요한 규제 문서에는 CP, 일반적인 보고에는 윌슨, 손으로 계산해야 하면 AC다. 왈드는 어느 칸에도 들어가지 못한다.

중심을 정리한다.

분자에 \(z^2/2\), 분모에 \(z^2\)이 더해졌다. \(z^2 = 1.96^2 = 3.84\)이므로 성공 \(1.92\)개와 실패 \(1.92\)개를 가상으로 추가한 것과 같다.

아그레스티–쿨은 이 \(1.92\)를 \(2\)로 반올림한 것이다. \(\tilde p = (k+2)/(n+4)\)로 두고 왈드 공식을 쓴다. 계산이 훨씬 쉬우면서 윌슨과 비슷한 성능을 낸다(연습문제 7에서 폭 \(+2.9\%\), 포함률 \(0.958\)).

베이즈로 보면 사전분포를 더한 것이다. \(\text{Beta}(a,b)\) 사전분포에 이항 자료를 곱하면 사후분포가 \(\text{Beta}(k+a,\ n-k+b)\)이고 사후평균이

이다. \(a=b=2\)면 아그레스티–쿨, \(a=b=1.92\)면 윌슨 중심이다.

왜 이것이 도움이 되는가. 가상 관측을 더하면 \(\hat p\)가 \(0\)이나 \(1\)이 되는 일이 사라진다. \(k=0\)이어도 \(\tilde p = 1.92/(n+3.84) > 0\)이므로 \(\sqrt{\tilde p(1-\tilde p)}\)가 \(0\)이 되지 않고, 구간이 붕괴하지 않는다. 왈드의 치명적 결함(연습문제 1)이 이 한 수로 제거된다.

제프리스가 특별한 이유. \(\text{Beta}(0.5,0.5)\)는 이항 모형의 피셔 정보량에서 자동으로 유도되는 사전분포다. 모수를 어떻게 다시 매개화해도 같은 답을 주는 불변성이 있어, "정보 없음"을 표현하는 원칙적인 방법으로 여겨진다. 구간 성능도 윌슨에 견줄 만하다.

빈도주의와 베이즈가 같은 자리에서 만난다. 윌슨은 순수하게 빈도주의적 유도(연습문제 5)에서 나왔는데, 결과가 베이즈 사후평균의 모양을 띤다. 좋은 절차는 두 철학에서 비슷한 곳에 도달하는 일이 많으며, 6장에서 이 현상을 다시 본다.

구간과 검정의 대응. "\(p_0\)를 유의수준 \(\alpha\)에서 기각하지 못하는 \(p_0\)의 집합"이 \(1-\alpha\) 신뢰구간이다. 어떤 검정통계량을 쓰느냐가 어떤 구간이 나오는지를 정한다.

핵심은 둘째 줄이다. 점수 검정은 표준오차를 \(\hat p\)가 아니라 \(p_0\)로 계산한다.

귀무가설 아래에서 분모가 상수이므로 \(Z\)의 분포가 정확히 표준정규에 가깝다. 반면 왈드는 분모가 확률변수라 분포가 어긋난다.

이것이 왈드가 극단에서 무너지는 이유의 정확한 설명이다. \(\hat p = 0\)이면 왈드의 분모가 \(0\)이 되어 통계량이 정의되지 않지만, 점수 검정의 분모 \(\sqrt{p_0(1-p_0)/n}\)은 \(p_0 > 0\)인 한 멀쩡하다. 그래서 윌슨 구간은 \(k=0\)에서도 \(\left[0,\ \frac{z^2}{n+z^2}\right]\)이라는 제대로 된 답을 준다(연습문제 2).

연습문제 5의 유도가 바로 이것이었다. \(|\hat p - p| \le z\sqrt{p(1-p)/n}\)에서 출발해 \(p\)에 대해 풀었는데, 이 부등식의 우변에 \(\hat p\)가 아니라 \(p\)가 들어 있다는 점이 점수 검정을 뒤집은 것과 같다. 이차부등식을 풀어야 해서 공식이 복잡해지지만, 그 복잡함이 곧 정확성의 근원이다.

일반 원리. 점수 검정은 귀무가설 아래의 정보만 쓰므로 모형이 잘 정의된 곳에서 통계량을 계산한다. 왈드 검정은 추정값 근처의 정보를 쓰므로 추정이 경계에 가까울 때 무너진다. 이 대비는 비율에 국한되지 않는다.

• 분산 성분이 \(0\) 경계에 있을 때
• 로지스틱 회귀에서 완전분리가 일어날 때(19장)
• 오즈비가 \(0\)이나 \(\infty\)에 가까울 때

모두 왈드형 통계량이 실패하고 점수형이나 우도비형이 살아남는다. "추정값에서 재지 말고 가설값에서 재라"가 기억할 규칙이다.

import numpy as np
# (연습문제 7의 ints, cover_width 를 그대로 쓴다)

ps = np.linspace(0.01, 0.99, 99)
print(f"{'방법':>9}{'평균':>10}{'최소':>10}{'RMSE(0.95 대비)':>18}")
rows = []
for i, name in enumerate(["Wald", "Wilson", "AC", "CP"]):
    cov = np.array([cover_width(40, p, i)[0] for p in ps])
    rows.append((name, cov.mean(), cov.min(),
                 np.sqrt(np.mean((cov - 0.95)**2))))
for name, a, b, c in rows:
    print(f"{name:>9}{a:>10.4f}{b:>10.4f}{c:>18.4f}")

출력:

       방법        평균        최소     RMSE(0.95 대비)
     Wald    0.8998    0.3310            0.1169
   Wilson    0.9508    0.9283            0.0102
       AC    0.9576    0.9382            0.0141
       CP    0.9706    0.9519            0.0231

기준에 따라 1위가 달라진다.

왈드는 세 기준 모두에서 꼴찌다. RMSE가 윌슨의 \(11\)배이고 최소 포함률은 \(0.33\)이다. 어떤 기준을 들이대도 선택할 이유가 없다는 본문의 결론이 정량적으로 확인된다.

CP와 윌슨의 선택은 진짜 맞바꿈이다.

• CP는 "절대 \(0.95\) 아래로 내려가지 않는다"를 보장하지만 평균이 \(0.971\)로 과하게 보수적이고, 구간이 \(9\%\) 넓다(연습문제 7).
• 윌슨은 평균적으로 정확하지만 어떤 \(p\)에서는 \(0.928\)까지 내려간다.

어느 쪽이 옳은지는 손실함수가 정한다. 규제 기관처럼 "명목 신뢰수준을 어긴 사례가 하나라도 있으면 안 된다"면 CP다. 연구 보고처럼 "전반적으로 정직한 불확실성 표현"이 목적이면 윌슨이다.

이산분포에서는 정확히 \(0.95\)가 불가능하다는 점도 기억하라. 포함률은 \(p\)에 대해 톱니 모양으로 오르내리므로(연습문제 4), 모든 \(p\)에서 정확히 \(0.95\)인 구간은 존재하지 않는다. 따라서 "어느 방향으로 틀릴 것인가"를 고르는 문제이며, 그 선택은 통계적이 아니라 실무적이다. \(\square\)