X̄의 표본분포 (Bernoulli)¶

(1) 이론값. \(\hat p\)는 \(0\)과 \(1\)의 표본평균이므로 \(\mu = p\), \(\sigma^2 = p(1-p)\)를 그대로 넣으면 된다.

왜도는 \(\gamma_1/\sqrt{n}\) 공식에 베르누이의 왜도 \((1-2p)/\sqrt{p(1-p)}\)를 넣어 얻는다.

왜도가 넷 다 \(\lvert \cdot \rvert \le 0.09\)로 아주 작다. \(n = 100\)에서 치우침은 사실상 없다는 예측이다. 그렇다면 정규곡선이 어긋날 자리도 없어야 할 텐데, 꼭 그렇지는 않다. 남은 것이 하나 있다.

\(n\hat p \sim \text{Binomial}(n, p)\)이므로 \(\hat p\)의 정확한 분포는 \(0,\, 0.01,\, 0.02,\, \ldots,\, 1\)에만 확률이 얹힌 이산분포다. 정규곡선은 연속이다. 이 두 가지를 누적분포로 견주면 어긋남의 크기가 곧바로 나오는데, 격자점 \(k/n\) 바로 아래에서 정규누적분포는 그 점의 확률을 절반쯤만 세므로 어긋남이 대략

이 된다. \(p = 0.5\), \(n = 100\)이면 \(0.3989/(2 \times 5) = 0.0399\)다. 이것은 치우침과 아무 상관이 없고 \(1/\sqrt{n}\)으로만 줄어든다. 격자를 반 칸 밀어 주는 연속성 보정을 넣으면 사라진다.

(2) 모의실험.

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

np.random.seed(1)

n_population = 10_000
n_sample = 100
n_sim = 1_000
p_values = [0.4, 0.5, 0.6, 0.7]

fig, axes = plt.subplots(1, len(p_values), figsize=(14, 3.5))

# 유한모집단 수정. 1만 개에서 100 개를 비복원으로 뽑으므로 폭이 조금 줄어든다.
fpc = np.sqrt((n_population - n_sample) / (n_population - 1))
print(f"유한모집단 수정 = {fpc:.5f}")
print("  p   실현비율   이론 SE   겨냥 SE   모의 평균   모의 표준편차   MC오차(평균/표준편차)")

# p를 0.4에서 0.7까지 바꿔 가며 네 패널을 그린다.
# 모집단은 0과 1뿐인 가장 비정규적인 분포인데도
# 표본비율의 표집분포는 어느 p에서나 종 모양이 된다.
for ax, p in zip(axes, p_values):
    population = stats.binom(n=1, p=p).rvs(n_population, random_state=1)

    # 0/1 자료의 평균이 곧 비율이므로 p-hat 은 표본평균의 한 경우다.
    p_hat_sims = np.array([
        np.random.choice(population, size=n_sample, replace=False).mean()
        for _ in range(n_sim)
    ])

    # 이론값과 모의값을 나란히 적는다. 비복원이므로 겨냥하는 것은 명목 p 가 아니라
    # 실현된 1만 개의 비율이고, 거기에 유한모집단 수정이 붙는다.
    p_real = population.mean()
    se = np.sqrt(p * (1 - p) / n_sample)
    se_target = np.sqrt(p_real * (1 - p_real) / n_sample) * fpc
    print(f"{p:.1f}    {p_real:.4f}    {se:.4f}    {se_target:.4f}    {p_hat_sims.mean():.4f}      "
          f"{p_hat_sims.std(ddof=1):.4f}        {se_target / np.sqrt(n_sim):.4f} / {se_target / np.sqrt(2 * n_sim):.4f}")

    # 모의실험으로 얻은 값들의 히스토그램.
    _, bins, _ = ax.hist(p_hat_sims, density=True, bins=15,
                         alpha=0.5, edgecolor="white",
                         label=r"simulated $\hat{p}$")

    # 정규근사를 겹쳐 그린다.
    # 베르누이의 분산이 p(1-p) 이므로 표준오차는 sqrt(p(1-p)/n) 이다.
    # 이 값은 p = 0.5 에서 최대가 되고 0이나 1에 가까울수록 작아진다.
    # 네 패널의 폭이 조금씩 다른 이유가 그것이다.
    se = np.sqrt(p * (1 - p) / n_sample)
    x_grid = np.linspace(bins[0], bins[-1], 200)
    pdf = stats.norm(loc=p, scale=se).pdf(x_grid)
    ax.plot(x_grid, pdf, "--r", lw=2, alpha=0.7, label="Normal approx.")
    ax.set_title(f"p = {p}")
    ax.set_xlabel(r"$\hat{p}$")

axes[0].set_ylabel("Density")
axes[-1].legend(fontsize=8)
plt.tight_layout()
plt.show()

# 근사의 오차를 정확한 이항분포와 견주어 잰다. 모의실험이 아니라 닫힌 꼴 계산이다.
print("\n  p   최대 |F_binom - F_normal|   보정 후   왜도     최빈값 확률/2   P(p-hat >= 0.5) 정확 / 근사 / 보정")
for p in p_values:
    k = np.arange(n_sample + 1)
    sd = np.sqrt(p * (1 - p) / n_sample)
    f_binom = stats.binom.cdf(k, n_sample, p)
    f_normal = stats.norm.cdf(k / n_sample, loc=p, scale=sd)
    f_corrected = stats.norm.cdf((k + 0.5) / n_sample, loc=p, scale=sd)
    print(f"{p:.1f}        {np.abs(f_binom - f_normal).max():.5f}            {np.abs(f_binom - f_corrected).max():.5f}"
          f"   {(1 - 2 * p) / np.sqrt(n_sample * p * (1 - p)):+.5f}   {stats.binom.pmf(k, n_sample, p).max() / 2:.5f}"
          f"     {stats.binom.sf(49, n_sample, p):.6f} / {stats.norm.sf(0.5, loc=p, scale=sd):.6f} / {stats.norm.sf(0.495, loc=p, scale=sd):.6f}")

출력:

유한모집단 수정 = 0.99504
  p   실현비율   이론 SE   겨냥 SE   모의 평균   모의 표준편차   MC오차(평균/표준편차)
0.4    0.3942    0.0490    0.0486    0.3909      0.0479        0.0015 / 0.0011
0.5    0.5010    0.0500    0.0498    0.5016      0.0503        0.0016 / 0.0011
0.6    0.6058    0.0490    0.0486    0.6038      0.0494        0.0015 / 0.0011
0.7    0.7049    0.0458    0.0454    0.7058      0.0460        0.0014 / 0.0010

  p   최대 |F_binom - F_normal|   보정 후   왜도     최빈값 확률/2   P(p-hat >= 0.5) 정확 / 근사 / 보정
0.4        0.04329            0.00272   +0.04082   0.04061     0.027099 / 0.020613 / 0.026240
0.5        0.03979            0.00027   +0.00000   0.03979     0.539795 / 0.500000 / 0.539828
0.6        0.03792            0.00272   -0.04082   0.04061     0.983238 / 0.979387 / 0.983956
0.7        0.03766            0.00578   -0.08729   0.04339     0.999991 / 0.999994 / 0.999996

중심과 폭. 비복원으로 뽑으므로 겨냥값은 명목 \(p\)가 아니라 실현된 1만 개의 비율이다. \(p = 0.4\)를 주문했는데 실제로 만들어진 모집단의 비율은 \(0.3942\)이고, 모의 평균 \(0.3909\)는 거기서 몬테카를로 오차 \(0.0015\)의 \(2.2\)배만큼 떨어져 있다. 나머지 셋은 각각 \(0.4\)배, \(1.3\)배, \(0.6\)배다. 폭도 네 경우 모두 겨냥 SE에서 표준편차의 몬테카를로 오차 \(0.0011\)의 \(0.6\)배, \(0.5\)배, \(0.7\)배, \(0.6\)배 안에 든다. 이론이 준 네 숫자가 그대로 재현된다.

되풀이가 1천 번뿐이라 오차가 앞 쪽들보다 세 배 크다. \(p = 0.4\)에서 평균이 \(2.2\) 오차만큼 비껴난 것은 그 때문이고, 되풀이를 늘리면 줄어든다.

어긋남의 원인 가리기 — 이것이 (2)의 본론이다. 아래 표의 셋째 열에 적힌 왜도는 넷 다 \(0.09\) 아래다. 그런데 정확한 이항 누적분포와 정규 누적분포의 차는 \(0.038\)–\(0.043\)으로, 왜도로 설명될 크기가 전혀 아니다. 범인은 다음 열이 가리킨다. 최빈값 확률의 절반이 \(0.0406\), \(0.0398\), \(0.0406\), \(0.0434\)로 어긋남과 거의 같은 값이다. \(p = 0.5\)에서는 \(0.03979\)로 소수 다섯째 자리까지 일치한다.

곧 어긋남은 치우침이 아니라 이산성이다. 연속성 보정을 넣으면 \(0.04\)대였던 차가 \(0.00027\)–\(0.00578\)로 한 자리에서 두 자리 줄어든다. 남은 것이 비로소 치우침 몫이며, 왜도가 가장 큰 \(p = 0.7\)에서 \(0.00578\)로 가장 크다는 순서까지 맞는다.

마지막 열이 이 구별을 한눈에 보여 준다. \(p = 0.5\)에서 \(P(\hat p \ge 0.5)\)는 정확히 \(0.539795\)다. \(\hat p = 0.5\)라는 점 하나에만 확률 \(0.0796\)이 얹혀 있고 그것을 전부 세기 때문이다. 정규근사는 \(0.5\)를 주어 \(0.0398\)만큼 틀리고, 보정을 넣으면 \(0.539828\)이 되어 넷째 자리까지 맞는다. 분포가 종 모양으로 보이는 것과 확률을 바르게 주는 것은 다른 문제다. 그림의 네 패널에서 히스토그램과 붉은 곡선이 잘 겹쳐 보이는 것은 15 개 구간으로 묶어 격자를 지워 놓았기 때문이기도 하다.

\(\hat{p} = \frac{1}{n}\sum_{i=1}^n X_i\)이고 \(X_i\)들이 독립이므로:

\(\square\)

정규근사를 사용하면 95% 신뢰구간은:

추정된 표준오차는:

\(z_{0.025} = 1.96\)이므로:

95% 신뢰구간은 약 \((0.481, 0.579)\)이다. 이 구간이 0.5를 포함하므로 95% 수준에서 이 후보가 과반의 지지를 받는다고 결론지을 수 없다. \(\square\)

\(p \in [0, 1]\)에서 \(g(p) = p(1-p) = p - p^2\)이라 하자.

\(g''(p) = -2 < 0\)이므로 최대점이다. 최댓값은:

따라서 표준오차는 다음을 만족한다:

이 상한은 표본크기를 계획할 때 유용하다. 미지의 \(p\)가 무엇이든 표준오차는 결코 \(1/(2\sqrt{n})\)을 넘지 않는다. 예를 들어 \(\text{SE} \le 0.03\)을 보장하려면 \(n \ge 1/(4 \times 0.03^2) \approx 278\)이 필요하다. \(\square\)

오차한계는 \(E = z_{0.025} \cdot \text{SE}(\hat{p}) = 1.96 \sqrt{p(1-p)/n}\)이다.

최악의 경우 \(p(1-p) \le 1/4\)를 사용하면:

\(E \le 0.02\)로 두면:

표본크기가 최소 \(n = 2401\)이면 오차한계 0.02 이하가 보장된다. \(\square\)

계산하면:

\(np = 1 < 5\)이므로 느슨한 기준조차 충족되지 않으며 정규근사를 신뢰할 수 없다. \(n\hat{p} = \sum X_i\)의 분포는 \(\text{Binomial}(100, 0.01)\)로 오른쪽으로 심하게 치우쳐 있고 0 근처에 몰려 있다.

대안으로는 다음이 있다:

• 정확한 binomial 방법: 신뢰구간과 검정에 정확한 binomial 분포를 사용한다(예: Clopper–Pearson 구간).
• 포아송 근사: \(n\)이 크고 \(p\)가 작으므로 \(\sum X_i \approx \text{Poisson}(\lambda = np = 1)\)이며 다루기가 더 간단한 경우가 많다.
• Wilson 구간: \(p\)가 0이나 1에 가까울 때 Wald(정규 기반) 구간보다 잘 작동하도록 수정된 신뢰구간이다.

일반적으로 관심 사건이 드물면 \(n\)을 크게 늘리거나 정규근사에 의존하지 않는 방법을 써야 한다. \(\square\)

\(\hat p \ge 0.5\)는 \(\sum X_i \ge 10\)과 같다.

(가) 정확한 값. \(\sum X_i \sim \text{Binomial}(20, 0.3)\)이므로

(나) 수정 없는 정규근사. \(\operatorname{SE}(\hat p) = \sqrt{0.3\times0.7/20} = 0.1025\)이므로

참값의 절반밖에 안 된다.

(다) 연속성 수정. 이산확률변수 \(X \ge 10\)을 연속 척도에서 \(X \ge 9.5\)로 바꾼다. \(E[X]=6\), \(\operatorname{SD}(X) = \sqrt{4.2} = 2.049\)이므로

참값 0.0480에 훨씬 가깝다.

정리.

\(np = 6\)이니 느슨한 기준 \(np \ge 5\)는 턱걸이로 넘긴 셈이다. 그런데도 오차가 이렇게 큰 것은 묻는 값이 분포의 가운데가 아니라 꼬리확률이기 때문이다. 느슨한 기준이 약속하는 것은 가운데의 모양뿐이고, \(n = 20\)에서 격자 간격이 여전히 성긴 탓에 꼬리에서는 한 칸의 차이가 그대로 드러난다. 그래도 연속성 수정 하나로 오차가 다섯 배 줄어든다. 이산분포를 연속분포로 근사할 때는 언제나 \(\pm0.5\) 보정을 넣어야 한다. 특히 꼬리 확률에서 차이가 크다.

델타 방법에 따르면 \(\operatorname{Var}\{g(\hat p)\} \approx \{g'(p)\}^2\operatorname{Var}(\hat p)\)이다.

\(g(p) = \arcsin\sqrt p\)를 미분하면 \(u = \sqrt p\)로 두어

이다. 따라서

로 \(p\)가 약분되어 사라진다. \(\square\)

왜 이런 변환을 찾는가. 일반적으로 \(\operatorname{Var}(\hat\theta) = \sigma^2(\theta)\)일 때

로 두면 분산이 상수가 된다. 여기서는 \(\sigma(p) = \sqrt{p(1-p)/n}\)이므로 적분이 \(\arcsin\sqrt p\)를 준다. 이런 변환을 분산안정화 변환이라 하며, 포아송의 \(\sqrt X\), 상관계수의 \(\operatorname{arctanh}\)(피셔 \(z\))도 같은 계산에서 나온다.

쓸모. 분산이 모수에 의존하지 않으면 (가) 등분산성이 필요한 분산분석·회귀에 비율 자료를 넣을 수 있고, (나) 신뢰구간의 폭이 일정해지며, (다) 여러 비율을 합치는 메타분석에서 가중치를 표본크기만으로 정할 수 있다.

한계. \(p\)가 0이나 1에 아주 가까우면 근사가 나빠진다. 델타 방법 자체가 \(\hat p\)가 \(p\) 근처에 머문다는 가정에 기대는데, 경계에서는 그렇지 않기 때문이다. 요즘은 변환 대신 이항 일반화선형모형을 직접 적합하는 편이 일반적이다.

전체 지지율은 \(p = 0.5\)다.

(가) 단순무작위추출.

표준오차가 0.025다.

(나) 층화추출. 층별 추정량을 \(\hat p_1, \hat p_2\)라 하면 \(\hat p_{\text{st}} = 0.5\hat p_1 + 0.5\hat p_2\)이고 두 층이 독립이므로

표준오차가 0.0245다.

차이의 정체. 단순무작위의 분산을 쪼개면

이다. 층화는 층 간 변동을 제거한다. 층을 고정된 크기로 뽑으므로 "남성이 우연히 많이 뽑히는" 변동이 아예 생기지 않기 때문이다.

이 예에서 이득이 4%로 작은 것은 층 간 차이(0.4 대 0.6)가 층 내 변동에 비해 작기 때문이다. 층 간 차이가 클수록 층화의 이득이 커진다. 반대로 층별 지지율이 모두 같으면 이득이 전혀 없다.

실무에서는 층 크기를 비례배분하지 않고 층별 분산에 비례해 배분하면(네이만 배분) 더 줄일 수 있다. 여기서는 두 층의 분산이 같아 비례배분이 이미 최적이다.

최대가능도. \(k = \sum x_i\)일 때 가능도는

이고

이다. \(\ell''(p) = -k/p^2 - (n-k)/(1-p)^2 < 0\)이므로 최대다.

최소분산불편성. 두 가지 길이 있다.

(길 1) 크라메르-라오. 관측값 하나의 피셔 정보량은

이므로 하한이 \(\frac{1}{nI(p)} = \frac{p(1-p)}{n}\)이다. \(\operatorname{Var}(\hat p) = p(1-p)/n\)이 정확히 이 값이므로 하한을 달성한다.

(길 2) 레만-셰페. 가능도를

로 쓰면 자료가 \(k = \sum x_i\)를 통해서만 들어오므로 \(\sum X_i\)가 충분통계량이고, 베르누이는 완비지수족이므로 완비이기도 하다. \(E[\bar X] = p\)로 불편이므로 \(\bar X\)가 유일한 최소분산불편추정량이다. \(\square\)

덧붙임. \(\hat p\)는 최적이지만 평균제곱오차 기준의 최적은 아니다. 예컨대 \(\tilde p = (k+1)/(n+2)\)(라플라스의 계승 규칙)는 편향되어 있지만 \(p\)가 0이나 1 근처일 때 MSE가 더 작고, 무엇보다 \(k=0\)일 때도 0이 아닌 값을 준다. 베이즈 관점에서는 \(\text{Beta}(1,1)\) 사전분포의 사후평균이다. 불편성이 언제나 최선의 기준은 아니라는 점을 다시 보여 준다.

\(\hat p = 10/10^6 = 10^{-5}\)이고

이다. 절댓값으로 보면 대단히 작아 보인다. 그러나 변동계수는

으로 32% 다. 추정값이 참값의 두 배이거나 절반이어도 이상하지 않다.

이유. 사건 수 \(k = n\hat p\)가 근사적으로 \(\text{Poisson}(np)\)을 따르고, 포아송의 변동계수는

이다. 여기서는 \(np = 10\)이므로 \(1/\sqrt{10} = 0.316\)이다. 정밀도를 정하는 것은 \(n\)이 아니라 \(np\), 즉 실제로 관측된 사건의 개수다.

상대오차를 10%로 줄이려면 \(np = 100\), 즉 표본을 1000만 명으로 늘려야 한다. 100명을 찾아야 한다는 뜻이다.

실무적 함의. 희귀 사건 연구에서 전체 표본크기를 자랑하는 것은 의미가 없다. 신뢰구간도 정규근사가 아니라 포아송 정확 구간을 써야 한다. \(k=10\)에 대한 95% 정확 구간은 대략 \((4.8,\ 18.4)\)건이므로 유병률 구간이 \((4.8, 18.4)\times10^{-6}\)이고 비대칭이다.

같은 이유로 희귀 사건 연구에서는 환자-대조군 설계가 쓰인다. 사건이 일어난 사람을 먼저 모으고 대조군을 맞추면, 전체 인구를 훑는 것보다 훨씬 적은 표본으로 같은 정밀도를 얻는다.