p̂의 표본분포 (표본크기의 효과)¶

(1) 해석적으로. \(n\hat p \sim \text{Binomial}(5, 0.3)\)이므로 점질량은 이항 확률 그대로다.

여섯 개를 손으로 적으면

가 된다. \(0.7^5 = 0.16807\)과 \(5(0.3)(0.7)^4 = 0.36015\)가 첫 두 개다. 합이 \(1\)인 여섯 개의 수가 분포 전부이며, 여기에는 "근사"라 할 것이 전혀 없다.

최대 점질량은 국소극한정리에서 어림한다. 최빈값이 \(k \approx np\)에 있고 그 자리에서 정규밀도가 최대 \(\phi(0)/\operatorname{SE}\)이므로

이다. \(1/\sqrt n\)으로 줄어든다. 격자 간격과 표준오차의 비도 같은 차수다.

그래서 "격자가 얼마나 거친가"를 재는 양이 \(\sqrt{npq}\) 하나로 요약된다. \(n = 5\)에서 \(\sqrt{npq} = 1.025\)이고 \(n = 1000\)에서 \(14.49\)다.

(2) 모의실험.

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

p = 0.3

print("n = 5 의 점질량 전부:")
for k in range(6):
    print(f"  P(p-hat = {k / 5:.1f}) = {stats.binom(5, p).pmf(k):.6f}")

fig, axes = plt.subplots(2, 2, figsize=(12, 6))
print("\n  n   격자 1/n      SE   격자/SE   최대 점질량   0.39894/sqrt(npq)")
for ax, n, ms in zip(axes.ravel(), (5, 20, 100, 1000), (7, 5, 3, 1.2)):
    se = np.sqrt(p * (1 - p) / n)

    # p-hat 의 정확한 분포. 이항 PMF 를 격자 k/n 위에 올린 점질량이다.
    k = np.arange(n + 1)
    pmf = stats.binom(n, p).pmf(k)
    x = k / n

    print(f"{n:5d} {1 / n:8.3f} {se:8.5f} {1 / (n * se):8.3f} {pmf.max():13.6f} "
          f"{0.398942 / np.sqrt(n * p * (1 - p)):17.6f}")

    lo, hi = max(0.0, p - 4.5 * se), min(1.0, p + 4.5 * se)
    m = (x >= lo) & (x <= hi)
    ax.vlines(x[m], 0, pmf[m], color="tab:blue", lw=1.0 if n < 500 else 0.4,
              alpha=0.7, label=r"exact PMF of $\hat p$")
    ax.plot(x[m], pmf[m], "o", color="tab:blue", ms=ms)

    # 정규 밀도를 n 으로 나누면 확률의 눈금이 된다(국소극한정리).
    g = np.linspace(lo, hi, 400)
    ax.plot(g, stats.norm(p, se).pdf(g) / n, "--r", lw=1.8,
            label=r"$N(p,\,pq/n)$ / $n$")

    ax.set_title(f"n = {n},  spacing 1/n = {1/n:g},  SE = {se:.4f}", fontsize=10)
    ax.set_xlim(lo, hi)
    ax.set_ylim(bottom=0)
    ax.set_xlabel(r"$\hat p$")
    ax.set_ylabel("Probability")

axes[0, 0].legend(fontsize=8)
plt.tight_layout()
plt.show()

# 국소극한정리를 한 점에서 확인한다.
n, k = 100, 30
exact = stats.binom(n, p).pmf(k)
approx = stats.norm(p, np.sqrt(p * (1 - p) / n)).pdf(k / n) / n
print(f"\n국소극한정리  n = {n}, k = {k}:  정확 {exact:.6f},  (1/n) f(k/n) = {approx:.6f},  "
      f"상대오차 {approx / exact - 1:+.4%}")

출력:

n = 5 의 점질량 전부:
  P(p-hat = 0.0) = 0.168070
  P(p-hat = 0.2) = 0.360150
  P(p-hat = 0.4) = 0.308700
  P(p-hat = 0.6) = 0.132300
  P(p-hat = 0.8) = 0.028350
  P(p-hat = 1.0) = 0.002430

  n   격자 1/n      SE   격자/SE   최대 점질량   0.39894/sqrt(npq)
    5    0.200  0.20494    0.976      0.360150          0.389328
   20    0.050  0.10247    0.488      0.191639          0.194664
  100    0.010  0.04583    0.218      0.086784          0.087056
 1000    0.001  0.01449    0.069      0.027521          0.027530

국소극한정리  n = 100, k = 30:  정확 0.086784,  (1/n) f(k/n) = 0.087056,  상대오차 +0.3140%

점질량은 손계산과 똑같다. 여섯 개가 소수 여섯째 자리까지 일치한다.

최대 점질량의 어림이 빠르게 좋아진다. \(0.39894/\sqrt{npq}\)가 실제 최대 점질량과 견주어 \(n = 5\)에서 \(8\%\) 크고(\(0.389\) 대 \(0.360\)), \(n = 20\)에서 \(1.6\%\), \(n = 100\)에서 \(0.31\%\), \(n = 1000\)에서 \(0.03\%\) 크다. 표본크기에 거의 반비례해 줄어든다(\(1/n\) 차수다). \(n = 5\)에서 어긋나는 것은 그 자리에서 최빈값 \(k = 1\)이 \(np = 1.5\)에서 반 칸 비껴나 있기도 해서다.

국소극한정리. \(n = 100\), \(k = 30\)에서 정확값 \(0.086784\)와 근사값 \(0.087056\)이 \(+0.31\%\) 차이다. 격자 한 칸의 폭 \(1/n\)에 그 자리의 정규밀도를 곱한 것이 그 칸의 확률이라는 해석이 세 자리까지 맞는다는 뜻이다. 이 관계가 있어야 점질량과 밀도곡선을 같은 축에 올릴 수 있다.

네 패널. 넷째 열의 "격자/SE"가 그림을 그대로 설명한다. \(n = 5\)에서 그 값이 \(0.976\)이다. 격자 한 칸이 표준오차만큼 크다는 뜻이므로 분포 전체가 막대 두세 개 폭에 들어가고, 붉은 곡선이 그 위를 지나가기는 하지만 두 대상을 같다고 할 수 없다. \(\hat p = 0\)에 아직 \(0.168\)이 남아 있는 것도 그 때문이다.

\(n = 20\)에서 비가 \(0.488\)로 절반이 되고 곡선이 꼭대기를 잘 따라간다. \(n = 100\)에서 \(0.218\)이면 \(p \pm \operatorname{SE}\) 안에 격자점이 \(2\sqrt{npq} \approx 9\)개 들어가므로 점들이 곡선 위에 얹혀 보인다. \(n = 1000\)에서는 \(0.069\), 격자점 \(29\)개여서 눈으로는 면과 구별되지 않는다.

눈이 속는 지점이 바로 여기다. 그림만 보면 \(n = 100\)에서 근사가 끝난 것 같다. 그러나 격자/SE가 \(0.218\)이라는 것은 CDF가 여전히 한 칸마다 최대 점질량 \(0.087\)만큼 뛰어오른다는 뜻이기도 하다. 다음 보기에서 그 값이 확률의 오차로 얼마가 되는지 잰다.

(1) 두 항을 따로 어림한다.

격자 항. 정확한 CDF는 격자점마다 그 점의 확률만큼 뛰는 계단이고 정규 CDF는 매끄럽다. 계단의 왼쪽 끝에서 재면 정규 CDF가 그 칸의 확률을 거의 절반 세므로, 어긋남은 가장 큰 점질량의 절반 정도다. 앞 보기에서 최대 점질량이 \(0.39894/\sqrt{npq}\)였으니

이다.

왜도 항. 에지워스 전개의 1차 보정은

이다. 보정항의 크기는 \(h(z) = (z^2-1)\phi(z)\)가 정하는데, \(h'(z) = z\phi(z)(3 - z^2)\)이므로 임계점이 \(z = 0, \pm\sqrt3\)이고

에서 \(z = 0\)이 최대다. 따라서

두 항의 비. 둘 다 최대가 \(z = 0\) 근처에서 나므로 수정 없는 오차는 둘의 합, 수정한 오차는 왜도 항만 남는다. 그러므로

이다. \(\sqrt{npq}\)가 통째로 약분되어 \(n\)이 사라진다. \(p = 0.3\)에 넣으면 \(1 + 3/0.4 = 8.5\)다. 연속성 수정이 \(n\)과 무관하게 오차를 \(8.5\)배 줄인다는 예측이고, 이것이 본문에서 "8배 이상"이라 한 수의 출처다.

(2) 수치로.

import numpy as np
from scipy import stats

p, q = 0.3, 0.7

print(" n     1/n    #values      SE     skew   sup|F-Phi|   with c.c.   격자항 예측   왜도항 예측")
for n in (5, 20, 100, 1000):
    se = np.sqrt(p * q / n)
    skew = (1 - 2 * p) / np.sqrt(n * p * q)

    # 격자점 k/n 에서 정확한 CDF 와 두 정규근사를 견준다.
    k = np.arange(n + 1)
    F_exact = stats.binom(n, p).cdf(k)
    F_plain = stats.norm.cdf((k / n - p) / se)                       # 수정 없음
    F_cc = stats.norm.cdf((k + 0.5 - n * p) / np.sqrt(n * p * q))    # 연속성 수정

    d_plain = np.max(np.abs(F_exact - F_plain))
    d_cc = np.max(np.abs(F_exact - F_cc))

    # 에지워스가 예측하는 두 항. phi(0)/2 = 0.19947, phi(0)/6 = 0.06649 이다.
    grid_term = stats.norm.pdf(0) / 2 / np.sqrt(n * p * q)
    skew_term = stats.norm.pdf(0) / 6 * skew
    print(f"{n:>4} {1/n:>7.3f} {n+1:>8} {se:>9.5f} {skew:>8.4f} "
          f"{d_plain:>10.4f} {d_cc:>11.4f} {grid_term:>13.4f} {skew_term:>13.4f}")
print(f"\n두 오차의 비 예측 = 1 + 3/(1-2p) = {1 + 3 / (1 - 2 * p):.2f}")

출력:

 n     1/n    #values      SE     skew   sup|F-Phi|   with c.c.   격자항 예측   왜도항 예측
   5   0.200        6   0.20494   0.3904     0.2154      0.0282        0.1947        0.0260
  20   0.050       21   0.10247   0.1952     0.1080      0.0127        0.0973        0.0130
 100   0.010      101   0.04583   0.0873     0.0491      0.0058        0.0435        0.0058
1000   0.001     1001   0.01449   0.0276     0.0156      0.0018        0.0138        0.0018

두 오차의 비 예측 = 1 + 3/(1-2p) = 8.50

왜도 항이 거의 완벽하게 맞는다. 마지막 열과 "with c.c." 열을 견주면 \(n = 20\)에서 \(0.0130\) 대 \(0.0127\), \(n = 100\)에서 \(0.0058\) 대 \(0.0058\), \(n = 1000\)에서 \(0.0018\) 대 \(0.0018\)이다. \(n = 5\)만 \(0.0260\) 대 \(0.0282\)로 \(8\%\) 작은데, 에지워스의 다음 항(첨도, \(n^{-1}\) 차수)이 그 자리에서는 아직 무시할 수 없기 때문이다.

격자 항도 맞는다. "sup\(\lvert F-\Phi \rvert\)"에서 "with c.c."를 뺀 값이 \(n = 5, 20, 100, 1000\)에서 \(0.1872\), \(0.0953\), \(0.0433\), \(0.0138\)이고, 예측한 격자 항이 \(0.1947\), \(0.0973\), \(0.0435\), \(0.0138\)이다. \(n = 20\)부터는 소수 넷째 자리까지 맞는다.

비도 맞는다. 실제 비가 \(7.63\), \(8.47\), \(8.50\), \(8.50\)이다. 예측 \(8.50\)에 \(n = 20\)부터 붙어 있고, \(n\)이 커져도 더 좋아지지 않는다. 연속성 수정의 효용은 표본크기에 달린 것이 아니라 \(p\)에 달려 있다. \(p\)가 \(1/2\)에 가까우면 \(1-2p\)가 작아져 비가 커지고(\(p = 0.45\)에서 \(31\)배), 극단적인 \(p\)에서는 작아진다.

이 수들을 확률로 읽으면 사정이 분명해진다. 수정하지 않은 정규근사의 최대 오차가 \(n = 100\)에서 \(0.0491\)이다. 그림에서는 완벽해 보였는데 유의수준 \(0.05\)를 다루는 자리에서 오차가 \(0.05\)라는 뜻이다. 같은 \(n\)에서 수정 한 줄을 넣으면 \(0.0058\)이 된다. 수정 없이 \(0.0058\)에 닿으려면 \(n\)을 \(100\)배 가까이 키워야 하므로, 표본을 더 모으는 것보다 공식을 고치는 쪽이 압도적으로 싸다.

가능한 값은 \(0, 0.05, 0.10, \ldots, 0.95, 1\)의 21개다. 격자 간격은 \(1/20 = 0.05\)이고

이므로 격자 한 칸이 표준오차의 \(0.05/0.10247 = 0.488\), 즉 절반쯤 된다. 분포가 퍼진 폭 안에 격자점이 네댓 개밖에 들어가지 않는다는 뜻이다.

\(\hat p = 0.31\)은 불가능하다. \(20 \times 0.31 = 6.2\)는 정수가 아니므로 그런 표본은 존재하지 않는다. \(\hat p\)의 표본분포를 연속으로 취급하면 이런 값에 0이 아닌 확률을 주게 된다. \(n = 20\)인 조사에서 "지지율 31%"라는 보고는 있을 수 없다.

\(\hat p \le 0.4\)는 \(X \le 8\)과 같다(\(X \sim \text{Binomial}(20, 0.3)\)).

(가) 정확한 값. \(P(X \le 8) = 0.8867\)이다.

(나) 수정 없는 정규근사. \(\text{SE} = 0.10247\)이므로

(다) 연속성 수정. \(np = 6\), \(\sqrt{npq} = 2.0494\)이므로

수정 없는 근사가 5퍼센트포인트 어긋나는 것을 수정 한 번으로 0.2퍼센트포인트로 줄였다. 어긋남의 방향에 주목할 만하다. 수정하지 않으면 \(P(X \le k)\)를 과소평가한다. 계단의 왼쪽 끝에서 재기 때문이다.

왜도는 위치와 척도에 불변이므로 \(X = n\hat p \sim \text{Binomial}(n, p)\)의 왜도를 구하면 된다. \(X = \sum_i X_i\)이고 \(X_i\)는 i.i.d. 베르누이다.

먼저 베르누이 하나의 3차 중심적률을 구한다. \(\mu = p\)이므로

독립인 확률변수의 합에서 3차 중심적률은 더해진다(누적률의 가법성). 따라서

이고

이다. \(\square\)

읽는 법. 분자 \(1-2p\)는 모집단의 비대칭을, 분모 \(\sqrt{npq}\)는 표본크기의 완화 효과를 담는다. \(p = 1/2\)이면 분자가 0이 되어 어떤 \(n\)에서도 왜도가 정확히 0이다. 이항분포가 \(p = 1/2\)에서 완벽히 대칭이라는 사실과 같은 말이다.

\(\gamma_1 = (1-2p)/\sqrt{npq} < 0.1\)을 \(n\)에 대해 풀면

이다.

\(p = 0.3\): \((0.4)^2 / (0.01 \times 0.21) = 0.16/0.0021 = 76.19\)이므로 \(n \ge 77\)이다. 검산하면 \(n = 77\)에서 \(\gamma_1 = 0.0995\)다.

\(p = 0.02\): \((0.96)^2 / (0.01 \times 0.0196) = 0.9216/0.000196 = 4702.0\)이므로 \(n \ge 4703\)이다.

60배 차이다. 같은 수준의 대칭성을 얻는 데 필요한 표본이 \(p\)에 따라 이렇게 달라진다. 이유는 식에 그대로 드러나 있다. \(p\)가 0에 가까우면 분자는 1에 가까워지고 분모의 \(pq\)는 0으로 가므로, 두 방향으로 동시에 나빠진다.

이 계산이 다음 페이지의 주제를 예고한다. "\(n\)이 크면 정규근사가 된다"는 말은 \(p\)를 고정하고 \(n \to \infty\)를 보낼 때만 참이다. 실무에서 \(n\)은 유한하고 \(p\)는 작을 수 있으므로, 판단의 기준이 \(n\) 하나가 아니라 \(np\)여야 한다.

구간 \([p - \text{SE},\ p + \text{SE}]\)의 길이는 \(2\text{SE} = 2\sqrt{pq/n}\)이고 격자 간격은 \(1/n\)이므로 개수는 대략

이다. \(\square\)

\(\sqrt n\)의 속도로만 늘어난다. \(n\)을 100배 키워야 격자점이 10배 촘촘해진다는 뜻이다. \(n = 1000\)에서도 표준오차 하나의 폭 안에 격자점이 29개뿐이다.

이것이 "이산성은 사라지지 않는다"는 말의 정량적 내용이다. 연속분포에서 나온 \(\bar X\)라면 이 개수가 처음부터 무한이다. \(\hat p\)에서는 유한하고, 그 유한한 개수가 느리게 늘어난다. 그래서 정확한 계산과 연속성 수정이 큰 \(n\)에서도 쓸모가 있다.

차수 확인. \(n\)이 \(k\)배가 되면 오차가 \(\sqrt k\)배로 줄어야 한다.

둘 다 \(n^{-1/2}\)다. 베리-에센 정리가 보장하는 차수이며, 연속성 수정은 이 차수를 바꾸지 못한다.

그러면 수정은 무엇을 없앴는가. 상수를 줄였다. 두 열의 비를 보면

로 일정하게 8.5배쯤이다. 에지워스 전개에서 \(n^{-1/2}\) 항은 두 조각이다. 하나는 격자의 계단에서 오는 톱니 모양 항이고 다른 하나는 왜도 항이다. 반 칸 밀기가 앞의 조각을 상쇄한다.

남은 오차를 왜도로 예측한다. 왜도 항은

이고 \(z\)에 대한 최댓값은 \(z = 0\)에서 \(\phi(0)\gamma_1/6 = 0.3989\,\gamma_1/6\)이다. 계산하면

거의 정확히 맞는다. 연속성 수정을 한 뒤 남는 오차가 전적으로 왜도에서 온다는 사실이 수치로 확인된다.

실무적 결론은 두 가지다. 첫째, 연속성 수정을 한 뒤에도 오차가 남고 그 크기는 \(\gamma_1/15\) 정도로 어림할 수 있다. 둘째, 그 오차를 더 줄이려면 왜도를 줄여야 하며 방법은 \(n\)을 키우거나 정확한 이항 계산을 쓰는 것뿐이다. 왜도가 큰 상황, 즉 \(p\)가 극단인 상황이 다음 페이지의 주제다.

델타 방법을 적용한다. \(g(p) = \arcsin\sqrt p\)이면

이므로

이다. \(p\)가 통째로 약분된다. \(\square\)

import numpy as np

rng = np.random.default_rng(0)
n = 100
print(f"{'p':>7}{'Var(phat)':>13}{'Var(arcsin)':>14}{'이론 1/(4n)':>14}")
for p in (0.05, 0.2, 0.5, 0.8, 0.95):
    ph = rng.binomial(n, p, 300_000) / n
    print(f"{p:>7}{ph.var():>13.6f}"
          f"{np.arcsin(np.sqrt(ph)).var():>14.6f}{1/(4*n):>14.6f}")

출력:

      p    Var(phat)   Var(arcsin)     이론 1/(4n)
   0.05     0.000474      0.002842      0.002500
    0.2     0.001604      0.002553      0.002500
    0.5     0.002505      0.002531      0.002500
    0.8     0.001595      0.002543      0.002500
   0.95     0.000476      0.002845      0.002500

\(\operatorname{Var}(\hat p)\)는 \(p\)에 따라 \(5\)배 넘게 변하는데(\(0.00047 \sim 0.00251\)) 변환 후에는 \(0.00253 \sim 0.00285\)로 거의 일정하다.

이것을 분산안정화 변환이라 한다. 분산이 평균에 의존하는 것이 문제가 되는 자리가 많다.

• 분산분석과 회귀는 등분산을 가정한다. 비율 자료를 그대로 넣으면 \(p\)가 극단인 집단의 분산이 작아 가중이 어긋난다.
• 관리도에서 관리한계를 그리려면 분산이 일정해야 편하다.
• 메타분석에서 여러 연구의 비율을 합칠 때 가중치가 \(p\)에 의존하지 않아 다루기 쉽다.

끝에서는 여전히 어긋난다. \(p = 0.05\)와 \(0.95\)에서 \(0.00284\)로 이론값보다 \(14\%\) 크다. 변환이 완벽하지 않은 것은 델타 방법이 1차 근사이기 때문이며, \(p\)가 극단이면 \(\hat p = 0\) 같은 경계 사건이 생겨 근사가 나빠진다.

대안으로 프리먼–투키 변환

이 있으며, 경계 근처에서 아크사인보다 안정적이다. 다만 오늘날에는 변환보다 이항 분포를 직접 모형화하는 쪽(로지스틱 회귀, 19장)이 표준이다. 변환은 컴퓨터가 없던 시절 정규 이론 도구를 쓰기 위한 우회로였고, 지금은 우회할 이유가 줄었다.

이웃한 두 확률의 비를 본다.

이 비가 \(1\)보다 크면 확률이 아직 오르고 있다는 뜻이다. 부등식을 풀면

이므로 \(k\)가 \((n+1)p\)를 넘는 순간부터 확률이 줄기 시작한다. 따라서 최빈값은 \(\lfloor (n+1)p \rfloor\)다. \(\square\)

import numpy as np
from scipy import stats

print(f"{'n':>5}{'p':>7}{'np':>8}{'최빈값':>8}{'floor((n+1)p)':>15}{'중앙값':>8}")
for n, p in [(20, 0.3), (20, 0.5), (25, 0.4), (7, 0.6)]:
    k = np.arange(n + 1)
    pmf = stats.binom.pmf(k, n, p)
    print(f"{n:>5}{p:>7}{n*p:>8.1f}{k[pmf.argmax()]:>8}"
          f"{int(np.floor((n+1)*p)):>15}{int(stats.binom.ppf(0.5, n, p)):>8}")

출력:

    n      p      np     최빈값  floor((n+1)p)     중앙값
   20    0.3     6.0       6              6       6
   20    0.5    10.0      10             10      10
   25    0.4    10.0      10             10      10
    7    0.6     4.2       4              4       4

셋이 대체로 붙어 있다. 평균 \(np\), 중앙값, 최빈값이 같거나 \(1\) 이내로 차이 난다. 일반적으로 이항분포에서는

임이 알려져 있다.

\((n+1)p\)이지 \(np\)가 아니라는 점이 요점이다. \(n=7\), \(p=0.6\)이면 \(np = 4.2\)이고 \((n+1)p = 4.8\)인데 둘 다 바닥이 \(4\)라 같은 답이 나온다. 그러나 \((n+1)p\)가 정수가 되는 경우에는 최빈값이 둘이다. \(n=9\), \(p=0.5\)면 \((n+1)p = 5\)이므로 \(k=4\)와 \(k=5\)의 확률이 정확히 같다.

왜 이 셋의 근접이 중요한가. 2장에서 본 대로 평균·중앙값·최빈값이 가까우면 분포가 거의 대칭이라는 신호다. 이항분포의 왜도가 \((1-2p)/\sqrt{npq}\)(연습문제 3)로 \(n\)이 커지면 \(0\)으로 가므로, 셋이 점점 붙는다.

반대로 \(p\)가 극단이면 벌어진다. \(n=20\), \(p=0.02\)면 평균 \(0.4\), 최빈값 \(0\), 중앙값 \(0\)이다. 평균이 최빈값과 중앙값에서 떨어져 나가며, 이것이 다음 페이지에서 다룰 왜도 문제의 다른 얼굴이다.

import numpy as np
from scipy import stats

n, p = 100, 0.3
print(f"{'c':>7}{'정확':>13}{'정규근사':>13}{'절대오차':>12}{'상대오차':>12}")
for c in (0.35, 0.40, 0.45, 0.50, 0.55):
    k = int(np.ceil(c * n))
    exact = stats.binom.sf(k - 1, n, p)
    z = (c - p - 0.5/n) / np.sqrt(p * (1 - p) / n)   # 연속성 수정
    approx = stats.norm.sf(z)
    print(f"{c:>7.2f}{exact:>13.3e}{approx:>13.3e}"
          f"{approx-exact:>+12.2e}{(approx-exact)/exact:>+12.2%}")

출력:

      c           정확         정규근사        절대오차        상대오차
   0.35    1.629e-01    1.631e-01   +1.96e-04      +0.12%
   0.40    2.099e-02    1.908e-02   -1.91e-03      -9.08%
   0.45    1.086e-03    7.777e-04   -3.08e-04     -28.38%
   0.50    2.206e-05    1.044e-05   -1.16e-05     -52.67%
   0.55    5.851e-08    4.488e-08   -1.36e-08     -23.30%

절대오차는 꼬리로 갈수록 작아진다. \(c=0.40\)에서 \(1.9\times10^{-3}\)이던 것이 \(c=0.55\)에서 \(1.4\times10^{-8}\)로 다섯 자릿수 줄어든다. 이것만 보면 "꼬리에서 근사가 더 좋다"고 착각하기 쉽다.

상대오차는 정반대로 간다. \(0.12\% \to 9\% \to 28\% \to 53\%\)로 폭증한다. \(c=0.50\)에서는 정규근사가 참값의 절반도 안 되는 값을 내놓는다(\(1.04\times10^{-5}\) 대 \(2.21\times10^{-5}\)). 확률 자체가 작아지므로 절대오차가 작아지는 것이 당연할 뿐이고, 정확도의 척도로는 상대오차가 맞다.

부호가 음수인 것도 읽어 둘 만하다. 정규근사가 꼬리 확률을 과소평가한다. \(p=0.3\)이라 분포가 오른쪽으로 치우쳐 있는데(왜도 \(+0.087\)) 정규는 대칭이라 오른쪽 꼬리를 얇게 잡기 때문이다.

어느 쪽을 봐야 하는가는 용도가 정한다.

\(p\)값에서 특히 문제가 된다. \(p = 0.001\)을 \(0.00105\)로 계산하면 절대오차는 \(5\times10^{-5}\)로 무시할 만하지만, 상대로는 \(5\%\) 틀린 것이다. 유의수준 근처에서 판정이 갈릴 수 있다.

극단 꼬리에서는 정규근사가 구조적으로 실패한다. 이항분포는 \(k > n\)에서 확률이 정확히 \(0\)인데 정규분포는 언제나 양수를 준다. \(c = 1.2\) 같은 값에서 정규근사는 \(10^{-20}\) 수준의 값을 내놓지만 참값은 \(0\)이므로 상대오차가 무한대다. 큰 편차 영역에서는 정규근사 대신 체르노프 경계나 정확 계산을 써야 한다.

import numpy as np
from scipy import stats

n, p = 50, 0.2
q = 1 - p
g1 = (1 - 2*p) / np.sqrt(n*p*q)          # 왜도
print(f"  n={n}, p={p}, 왜도 = {g1:.4f}\n")
print(f"{'k':>4}{'정확':>11}{'정규+연속성':>13}{'에지워스':>11}"
      f"{'수정오차':>11}{'에지워스오차':>13}")
for k in (5, 8, 10, 12, 15):
    exact = stats.binom.cdf(k, n, p)
    z = (k + 0.5 - n*p) / np.sqrt(n*p*q)
    cc = stats.norm.cdf(z)
    ed = cc - g1/6 * (z**2 - 1) * stats.norm.pdf(z)
    print(f"{k:>4}{exact:>11.6f}{cc:>13.6f}{ed:>11.6f}"
          f"{cc-exact:>+11.2e}{ed-exact:>+13.2e}")

출력:

  n=50, p=0.2, 왜도 = 0.2121

   k         정확       정규+연속성       에지워스       수정오차       에지워스오차
   5   0.048027     0.055806   0.049714  +7.78e-03    +1.69e-03
   8   0.307332     0.297942   0.306749  -9.39e-03    -5.82e-04
  10   0.583559     0.570158   0.583610  -1.34e-02    +5.08e-05
  12   0.813943     0.811620   0.813708  -2.32e-03    -2.35e-04
  15   0.969197     0.974085   0.968162  +4.89e-03    -1.03e-03

에지워스 보정이 오차를 크게 줄인다. 연속성 수정만 했을 때 \(10^{-2}\) 수준이던 오차가 대체로 \(10^{-3}\) 이하로 내려간다. \(k=10\)에서는 \(-1.34\times10^{-2}\)가 \(+5.1\times10^{-5}\)로 \(260\)배 줄었다.

다만 개선 폭이 일정하지는 않다. \(k=5\)와 \(k=15\)에서는 \(10^{-3}\) 수준이 남는데, 왜도 다음 항인 첨도가 아직 보정되지 않았기 때문이다. 에지워스 전개를 한 항 더 가져가면 이 부분도 줄어든다.

두 보정이 다른 것을 고친다.

연습문제 6이 "연속성 수정 뒤 남는 오차는 왜도 탓"이라고 했는데, 그 왜도를 직접 빼 주니 실제로 사라진다. 진단이 맞았음을 확인한 셈이다.

그런데 실무에서 에지워스를 쓰지는 않는다. 이항분포는 scipy.stats.binom.cdf로 정확히 계산할 수 있으므로 근사할 이유가 없다. 에지워스가 쓸모 있는 자리는 정확 계산이 불가능한 복잡한 통계량의 분포를 근사할 때이며, 부트스트랩의 정확도를 이론적으로 분석할 때(17장 BCa 구간의 근거) 핵심 도구가 된다.

여기서 얻을 교훈은 도구가 아니라 구조다. 정규근사의 오차는 무작위하게 흩어진 것이 아니라 왜도·첨도 같은 고차 적률로 설명되는 체계적인 것이다. 그래서 "얼마나 큰 \(n\)이면 되는가"라는 질문의 답이 모집단에 따라 달라지고, \(p\)가 극단일수록 커지는 것이다. \(\square\)