p의 신뢰구간¶

(1) 구간. \(\hat p = 120/200 = 0.6\)이고 \(z_{0.025} = 1.95996\)이다.

왜 산포를 따로 추정하지 않는가. 베르누이 분포에는 모수가 하나뿐이다. \(X \sim \text{Bernoulli}(p)\)이면

로 분산이 평균에 딸려 온다. 그래서 \(\hat p\) 하나가 구간의 중심과 폭을 모두 정한다. 정규모집단에서는 \(\mu\)와 \(\sigma^2\)이 따로 움직이는 두 모수라 \(\bar x\) 말고 \(s\)를 하나 더 추정해야 했고, 그 추정의 대가로 \(t\) 임계값을 치렀다. 비율에는 치를 대가가 없어 보인다.

그러나 공짜가 아니다. 분산을 따로 추정하지 않는 대신, 표준오차 안의 \(p\) 자리에 추정값 \(\hat p\)을 꽂았다. 그 꽂기가 \(\hat p\)이 극단일 때 구간을 무너뜨린다. 여기서는 그 일이 일어나지 않는다는 것이 타당성 조건이다.

보수적 기준 \(10\)을 넉넉히 만족한다.

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

이다. 모의실험이 필요 없다. 아래 코드가 이 합을 그대로 계산한다.

확인.

import numpy as np
import scipy.stats as stats

n = 200          # 표본크기
x = 120          # 성공 횟수
confidence_level = 0.95

p_hat = x / n
z_critical = stats.norm.ppf(1 - (1 - confidence_level) / 2)
# 평균의 구간과 달리 여기서는 산포를 따로 추정하지 않는다.
# Bernoulli 분포에서는 분산 p(1-p)가 평균 p에 딸려 오기 때문이다.
# 그래서 p_hat 하나로 중심과 너비가 모두 정해진다.
standard_error = np.sqrt((p_hat * (1 - p_hat)) / n)
margin_of_error = z_critical * standard_error
confidence_interval = (p_hat - margin_of_error, p_hat + margin_of_error)

print(f"{p_hat = }, {standard_error = :.7f}, {margin_of_error = :.7f}")
print(f"{confidence_interval = }")


# (2) 포함률은 유한합이다. 열거해서 더하면 끝이다.
def exact_coverage(n, p, z):
    k = np.arange(n + 1)
    ph = k / n
    pmf = stats.binom.pmf(k, n, p)
    half = z * np.sqrt(ph * (1 - ph) / n)
    return pmf[(ph - half <= p) & (p <= ph + half)].sum()


print("\n참값 p = 0.6 에서 명목 95% 왈드 구간의 정확한 포함률")
for m in (198, 199, 200, 201, 202, 205):
    print(f"  n = {m}:  {exact_coverage(m, 0.6, z_critical):.5f}")

출력:

p_hat = 0.6, standard_error = 0.0346410, margin_of_error = 0.0678951
confidence_interval = (0.5321048559554297, 0.6678951440445703)

참값 p = 0.6 에서 명목 95% 왈드 구간의 정확한 포함률
  n = 198:  0.94077
  n = 199:  0.94920
  n = 200:  0.94895
  n = 201:  0.94756
  n = 202:  0.94771
  n = 205:  0.95386

(1)과 맞는다. \(\hat p = 0.6\), \(\text{SE} = 0.0346410\), \(\text{MOE} = 0.0678951\), 구간 \((0.5321049,\ 0.6678951)\) 이 손 계산과 같다.

(2)의 답은 \(0.95\)가 아니라 \(0.94895\)다. 타당성 조건을 넉넉히 만족하고 \(p = 0.6\)이 경계에서 멀며 \(n\)이 \(200\)인데도, 명목 \(95\%\) 구간이 실제로는 \(94.90\%\)를 담는다. 모자라는 \(0.00105\)는 작지만 잡음이 아니다. 이 값은 모의실험이 아니라 유한합으로 구한 것이라 몬테카를로 오차가 아예 없다.

더 눈에 띄는 것은 \(n\)에 대한 움직임이다. \(n\)을 \(198\)에서 \(205\)까지 하나씩 올리는 동안 포함률이

으로 오르내린다. \(n = 198\)에서 \(199\)로 하나 늘렸을 뿐인데 \(0.0084\) 올랐고, 다시 \(200\)에서 \(201\)로 가며 내렸다. 포함률은 \(n\)의 증가함수가 아니다.

까닭은 위 합의 지시함수에 있다. \(n\)이 하나 늘면 격자 \(k/n\) 전체가 조금씩 움직이고, 그러다 어떤 \(k\)의 구간이 \(p = 0.6\)을 담던 상태에서 담지 못하는 상태로 갑자기 넘어간다. 그 \(k\)의 점질량이 통째로 합에서 빠지거나 들어온다. \(\hat p\)이 이산이라서 생기는 일이며, \(n\)을 키운다고 사라지지 않고 톱니가 잘아질 뿐이다.

이 쪽에서 다루는 \(p = 0.6\), \(n = 200\)은 왈드 구간에 가장 유리한 자리다. \(p\)가 극단으로 가면 사정이 훨씬 나빠지며, 그 모습은 아래 그림과 비율 신뢰구간의 실제 포함률에서 본다.

95% 신뢰수준에서 \(z_{\alpha/2} \approx 1.96\)이다.

그 후보를 지지하는 유권자의 참 비율이 0.5321과 0.6679 사이에 있다고 95% 신뢰한다.

최악의 표준오차는 (\(\hat{p}(1-\hat{p})\)를 최대화하는) \(\hat{p} = 0.5\)에서 나온다.

import scipy.stats as stats
import numpy as np

confidence_level = 0.95
alpha = 1 - confidence_level
z_star = stats.norm().ppf(1 - alpha / 2)
margin_of_error_max = 0.02
p_max = 0.5

# p를 모르는 채로 표본크기를 정해야 하므로 최악을 가정한다.
# p(1-p)는 p=0.5에서 최대(0.25)이니, 이 n이면 참 p가 무엇이든 안전하다.
n = 1
while True:
    n += 1
    me = z_star * np.sqrt(p_max * (1 - p_max) / n)
    if me <= margin_of_error_max:
        break
print(f"{n = }")

출력:

n = 2401

공식으로 계산하면 \(n = (z_{0.025}/2E)^2 = 2400.91\)이고 올림하여 2401이다(\(z = 1.96\)으로 반올림해 쓰면 정확히 2401이 된다). 참 비율이 0.5에서 멀면 이만큼 필요하지 않다. 사전 정보로 \(p \approx 0.1\)을 쓸 수 있다면 같은 오차한계에 865개면 된다.

import numpy as np
from scipy import stats

confidence_level = 0.99
alpha = 1 - confidence_level
p_hat = 20 / 50
n = 50

z_star = stats.norm().ppf(1 - alpha / 2)
margin_of_error = z_star * np.sqrt(p_hat * (1 - p_hat) / n)
print(f"{p_hat} ± {margin_of_error:.3f}")

출력:

0.4 ± 0.178

99% 신뢰구간은 대략 \((0.222, 0.578)\)이다.

여기서 한 가지 짚어 둘 것이 있다. 500곡 중 50곡을 뽑았으므로 \(n/N = 0.1\)이고, i.i.d. 근사가 아슬아슬한 경계에 있다. 유한모집단 수정을 넣으면 \(\sqrt{(500-50)/499} = 0.950\)이 곱해져 오차한계가 0.178에서 0.169로 줄어든다.

\(\hat p = 0.74\). \(\mathrm{SE} = \sqrt{0.74 \cdot 0.26/100} \approx 0.0439\). 오차한계 = \(1.96 \cdot 0.0439 \approx 0.086\).

신뢰구간: \((0.654, 0.826)\).

\(\hat p = 0.60\). \(z_{0.05} = 1.645\). 오차한계 = \(1.645 \sqrt{0.60 \cdot 0.40/200} \approx 0.057\).

신뢰구간: \((0.543, 0.657)\).

(a) \(\hat p = 0.56\). 오차한계 = \(1.96 \cdot 0.0222 = 0.0435\). Wald 신뢰구간: \((0.517, 0.604)\).

(b) Wilson 신뢰구간: 중심 \(= (0.56 + 1.96^2/1000)/(1 + 1.96^2/500) \approx 0.560\). 신뢰구간 \(\approx (0.516, 0.603)\). Wald와 거의 같다 — \(n\)이 크고 \(\hat p\)가 중간이면 Wilson과 Wald가 일치한다.

(c) 신뢰구간 전체가 0.5보다 위에 있다: 5% 수준에서 과반이 이 안건을 지지한다는 증거이다.

Wald: 오차한계 \(= 1.96 \sqrt{0.1 \cdot 0.9/20} \approx 0.131\). Wald 신뢰구간: \((0.1 - 0.131, 0.1 + 0.131) = (-0.031, 0.231)\) — 0 아래로 뻗는다.

Wilson: 중심 \(\approx (0.1 + 1.96^2/40)/(1 + 1.96^2/20) \approx 0.196/1.192 \approx 0.164\). 반너비: \(1.96\sqrt{0.1 \cdot 0.9/20 + 1.96^2/1600}/1.192 \approx 0.137\). 신뢰구간 \(\approx (0.028, 0.301)\) — \([0, 1]\) 안에 머문다.

경계 근처에서는 Wilson이 훨씬 낫다. 이항 신뢰구간에는 Wilson을 기본으로 삼으라.

\(\mathrm{ME} = 1.96\sqrt{\hat p(1-\hat p)/n} \le 0.03\).

최악의 경우는 \(\hat p = 1/2\)에서 \(\hat p(1-\hat p) = 0.25\)일 때이다.

\(1.96 \sqrt{0.25/n} \le 0.03 \Rightarrow n \ge (1.96)^2 \cdot 0.25/(0.03)^2 = 0.9604/0.0009 \approx 1068\).

참 \(p\)가 무엇이든 \(n = 1068\)이면 95% 신뢰수준에서 오차한계 \(0.03\) 이하가 된다. 여론조사의 "\(n \approx 1000\)" 규칙이 여기서 나온다.

\(p\)가 0.1이나 0.9 근처라고 짐작된다면 \(p(1-p) = 0.09\)이므로 \(n \ge 385\)면 된다. \(p\)의 대략적인 값을 알면 필요한 \(n\)을 약 64% 줄일 수 있다.

연속성 보정은 이항분포에 정규근사를 쓸 때 반너비에 \(1/(2n)\)을 더한다. 연습문제 1에서는 보정한 오차한계가 \(\approx 0.086 + 1/200 = 0.091\)이다.

보정한 신뢰구간: \((0.649, 0.831)\), 보정하지 않은 것: \((0.654, 0.826)\). 약간 더 넓다.

언제 쓰는가:

• 연속성이 문제가 되는 작거나 중간 크기의 \(n\).
• 꼬리 확률이 중요할 때(단측 검정).
• 이항분포의 이산성이 무시할 수 없을 때.

현대의 실무: 보정 대신 정확한 방법(Clopper-Pearson)이나 Wilson 구간을 쓴다. 연속성 보정은 컴퓨터 이전 시대 근사의 유산이다.

기본식. 오차한계 \(E=z\sqrt{p(1-p)/n}\)을 \(n\)에 대해 풀면

\(p\)를 모를 때. \(p(1-p)\)는 \(p=0.5\)에서 최대 0.25다. 가장 보수적인 선택이다.

import numpy as np
from scipy import stats

E = 0.03
print(f"{'신뢰수준':>8s} {'z':>7s} {'p=0.5':>8s} {'p=0.2':>8s} {'p=0.05':>8s}")
for lv in [0.90, 0.95, 0.99]:
    z = stats.norm.ppf(0.5 + lv / 2)
    ns = [int(np.ceil(z**2 * p * (1 - p) / E**2)) for p in [0.5, 0.2, 0.05]]
    print(f"{lv:8.0%} {z:7.4f} {ns[0]:8d} {ns[1]:8d} {ns[2]:8d}")

    신뢰수준       z    p=0.5    p=0.2   p=0.05
     90%  1.6449      752      481      143
     95%  1.9600     1068      683      203
     99%  2.5758     1844     1180      351

읽기.

• 95%에서 \(p\)를 모르면 \(n=1068\). 여론조사에서 "표본 1000명, 오차 \(\pm3.1\%\)포인트"라는 문구가 여기서 나온다.
• 사전정보가 큰 차이를 만든다. \(p\approx0.2\)임을 알면 683명으로 충분하다. 36% 절약이다.
• \(p=0.05\)처럼 극단이면 203명이면 되지만, 이때는 정규근사 자체가 나빠 윌슨 기준으로 다시 계산해야 한다.

\(E\)에 대한 민감도. \(n\propto E^{-2}\)이므로

오차를 절반으로 줄이려면 표본이 네 배다. 1%포인트를 노리면 1만 명이 필요하고, 이 지점에서 표집오차보다 무응답 편향이 훨씬 커지므로 투자할 가치가 없다.

주의. 이 계산은 단순임의추출을 가정한다. 실제 조사는 층화·군집 설계를 쓰므로 다음 문제에서 볼 설계효과를 곱해야 한다.

설계효과. 크기 \(m\)인 군집을 뽑을 때

이고, 유효표본크기는 \(n_{\text{eff}}=n/\text{DEFF}\)다.

import numpy as np
from scipy import stats

n, m, rho, phat = 1000, 25, 0.05, 0.42
deff = 1 + (m - 1) * rho
n_eff = n / deff
z = stats.norm.ppf(0.975)

se0 = np.sqrt(phat * (1 - phat) / n)
se1 = np.sqrt(phat * (1 - phat) / n_eff)
print(f"DEFF = 1 + ({m}-1)×{rho} = {deff:.2f}")
print(f"유효표본크기 {n_eff:.1f}명 (실제 {n}명)")
print(f"군집 무시  SE {se0:.5f}  구간 ({phat - z * se0:.4f}, {phat + z * se0:.4f})  "
      f"오차한계 {z * se0:.4f}")
print(f"설계 반영  SE {se1:.5f}  구간 ({phat - z * se1:.4f}, {phat + z * se1:.4f})  "
      f"오차한계 {z * se1:.4f}")
print(f"오차한계가 {se1 / se0:.2f}배로 커진다")

DEFF = 1 + (25-1)×0.05 = 2.20
유효표본크기 454.5명 (실제 1000명)
군집 무시  SE 0.01561  구간 (0.3894, 0.4506)  오차한계 0.0306
설계 반영  SE 0.02315  구간 (0.3746, 0.4654)  오차한계 0.0454
오차한계가 1.48배로 커진다

1000명을 조사했지만 정보량은 455명어치다. 같은 학교 학생들이 서로 비슷하기 때문에, 한 학교에서 25명을 뽑아도 25명분의 독립정보를 얻지 못한다.

\(\rho\)가 작아도 \(m\)이 크면 심각하다.

\((m-1)\rho\)가 관건이다. 급내상관이 0.01처럼 작아 보여도 군집이 100명이면 유효표본이 절반이 된다.

설계 함의. 같은 예산에서 군집을 많이, 각 군집은 작게 뽑는 것이 유리하다. 학교 40곳×25명보다 학교 100곳×10명이 낫다(DEFF가 2.20에서 1.45로 준다). 다만 학교를 추가하는 비용이 학생을 추가하는 비용보다 크므로, 최적 \(m\)은 비용비와 \(\rho\)로 결정된다.

보고. 조사 결과를 낼 때 설계효과나 유효표본크기를 함께 적어야 한다. "\(n=1000\), 오차 \(\pm3.1\%\)p"는 군집설계에서 거짓이다.

층화추정량.

여기서 \(W_h\)는 모집단에서의 층 비중이다.

import numpy as np
from scipy import stats

W = np.array([0.5, 0.3, 0.2])
nh = np.array([300, 200, 200])
ph = np.array([0.40, 0.55, 0.30])
z = stats.norm.ppf(0.975)

p_st = (W * ph).sum()
var_st = (W**2 * ph * (1 - ph) / nh).sum()
se_st = np.sqrt(var_st)

n = nh.sum()
p_pool = (nh * ph).sum() / n                 # 층화를 무시하고 단순 합산
se_pool = np.sqrt(p_pool * (1 - p_pool) / n)

print(f"층화 추정 {p_st:.4f}  SE {se_st:.5f}  "
      f"구간 ({p_st - z * se_st:.4f}, {p_st + z * se_st:.4f})")
print(f"단순 합산 {p_pool:.4f}  SE {se_pool:.5f}  "
      f"구간 ({p_pool - z * se_pool:.4f}, {p_pool + z * se_pool:.4f})")
print(f"설계효과 DEFF = {var_st / (p_st * (1 - p_st) / n):.4f}")

nh_prop = np.array([350, 210, 140])          # 비례배분이었다면
var_prop = (W**2 * ph * (1 - ph) / nh_prop).sum()
print(f"비례배분일 때 SE {np.sqrt(var_prop):.5f}  "
      f"DEFF = {var_prop / (p_st * (1 - p_st) / n):.4f}")

층화 추정 0.4250  SE 0.01880  구간 (0.3882, 0.4618)
단순 합산 0.4143  SE 0.01862  구간 (0.3778, 0.4508)
설계효과 DEFF = 1.0122
비례배분일 때 SE 0.01837  DEFF = 0.9668

두 가지가 다르다.

1. 추정값이 다르다. 0.4250 대 0.4143. 표본배분 \((300,200,200)\)이 모집단 비중 \((0.5,0.3,0.2)\)과 다르기 때문이다. 비례배분이 아니면 단순 합산은 편향된다. 1지역이 모집단의 50%인데 표본에서는 42.9%뿐이라 과소대표되고, 비율이 가장 낮은 3지역이 20% 대신 28.6%로 과대대표되어 합산값이 아래로 끌린다.

2. 표준오차도 다르다. 이 배분에서는 DEFF가 1.012로 1보다 크다. 층화 자체는 분산을 줄이지만, 배분이 비례에서 벗어난 손실이 그 이득을 조금 넘어섰다. 같은 700명을 비례배분 \((350,210,140)\)으로 뽑았다면 DEFF가 0.967로 내려가 단순임의추출보다 유리해진다.

층화의 두 얼굴.

• 편향 제거. 배분이 비례가 아닐 때 가중이 필수다. 이것은 언제나 이득이다.
• 분산. 비례배분이면 DEFF \(\le1\)이 보장되지만, 비례에서 벗어나면 1을 넘을 수 있다. 층별 추정을 위해 작은 층을 과대표집하면 전체 추정의 정밀도를 내주는 셈이다.

군집과 반대 방향이다. 앞 문제에서 군집은 DEFF를 1보다 크게 만들었다. 실제 조사는 층화와 군집을 함께 쓰므로 두 효과가 상쇄되며, 총 DEFF는 보통 1.5~3 사이다.

최적 배분. 분산을 최소화하는 네이만 배분은 \(n_h\propto W_h\sigma_h\)다. 비율이라면 \(\sigma_h=\sqrt{p_h(1-p_h)}\)이므로 \(p_h\)가 0.5에 가까운 층에 더 많이 배분한다. 다만 층별 추정도 필요하다면 최소 표본수를 보장해야 한다.

착안. \(\hat p\)는 \([0,1]\)에 갇혀 있어 경계 근처에서 정규근사가 나쁘다. 로짓

은 실수 전체를 값으로 가지므로 정규근사가 낫다. 델타법으로

구간을 로짓 척도에서 만들고 역로짓으로 되돌린다.

import numpy as np
from scipy import stats

n, k = 30, 3
z = stats.norm.ppf(0.975)
ph = k / n

# 로짓 구간
eta = np.log(ph / (1 - ph))
se_eta = np.sqrt(1 / k + 1 / (n - k))
lo_l, hi_l = [1 / (1 + np.exp(-(eta + s * z * se_eta))) for s in (-1, 1)]

# 왈드
se = np.sqrt(ph * (1 - ph) / n)
lo_w, hi_w = ph - z * se, ph + z * se

# 윌슨
d = 1 + z**2 / n
c = (ph + z**2 / (2 * n)) / d
h = z / d * np.sqrt(ph * (1 - ph) / n + z**2 / (4 * n**2))

# 클로퍼-피어슨
lo_c = stats.beta.ppf(0.025, k, n - k + 1)
hi_c = stats.beta.ppf(0.975, k + 1, n - k)

for name, lo, hi in [("왈드", lo_w, hi_w), ("로짓", lo_l, hi_l),
                     ("윌슨", c - h, c + h), ("클로퍼-피어슨", lo_c, hi_c)]:
    print(f"{name:12s} ({lo:.4f}, {hi:.4f})  폭 {hi - lo:.4f}")

왈드           (-0.0074, 0.2074)  폭 0.2147
로짓           (0.0326, 0.2681)  폭 0.2355
윌슨           (0.0346, 0.2562)  폭 0.2216
클로퍼-피어슨      (0.0211, 0.2653)  폭 0.2442

읽기.

• 왈드는 하한이 음수다. 확률이 음수일 수 없으므로 명백히 잘못이다.
• 로짓 구간은 \((0,1)\) 안에 있다. 역로짓이 항상 \((0,1)\)로 사상하기 때문에 구조적으로 보장된다.
• 로짓과 윌슨은 하한이 거의 같다(0.033 대 0.035). 상한은 로짓이 0.268로 윌슨의 0.256보다 높아, 로짓 구간이 조금 더 넓고 오른쪽으로 더 뻗는다.
• 클로퍼-피어슨이 가장 넓다. 보수적인 정확 구간이다.

로짓의 치명적 결함. \(k=0\)이나 \(k=n\)이면 \(\hat\eta=\mp\infty\)이고 분산도 무한대라 구간을 만들 수 없다. 흔한 대처는 \(k\)와 \(n-k\)에 0.5씩 더하는 것인데(하일랜드 수정), 임시방편이다.

어느 것을 쓰는가. 단독 비율이라면 윌슨이 가장 낫다. 로짓이 유용한 곳은 회귀모형이다. 로지스틱회귀의 계수 구간을 만들고 되돌리면 자동으로 로짓 구간이 되며, 여러 변수를 함께 다룰 수 있다.