p̂₁ - p̂₂의 표본분포 (정규근사가 듣는 경우)¶

(1) 네 값 모두 닫힌 꼴이다. 이항의 누적률 \(\kappa_2 = npq\), \(\kappa_3 = npq(q-p)\), \(\kappa_4 = npq(1-6pq)\)를 \(\hat p = X/n\)의 척도로 옮기면 \(\kappa_j(\hat p) = \kappa_j(X)/n^j\)이므로

이다. 누적률은 독립인 확률변수의 합에서 더해지고 부호에서 \(\kappa_j(-Y) = (-1)^j\kappa_j(Y)\)이므로, 차에서는 홀수 차수만 빼고 짝수 차수는 더한다. 중심은 그대로 \(p_1-p_2 = 0.10\)이고

이다. \(p_1 = 0.5\)인 집단은 대칭이라 \(\kappa_3\)에 아무것도 보태지 않는다. 표준화하면

둘 다 \(0.01\)에 못 미친다. 왜도 \(-0.0099\)는 지수분포의 \(2\)에 견주면 200분의 일이고, 표준정규의 표본 20만 개로 재어도 왜도 추정의 몬테카를로 오차가 \(\sqrt{6/200000} = 0.0055\)이니 이 크기는 잡음과 구별되기도 어렵다. 그림에서는 아무 어긋남도 보이지 않을 것이라고 예측된다.

눈에 보일 것이 하나 남는데 모양이 아니라 격자다. \(n_1 = n_2 = 200\)이면 \(\hat p_1-\hat p_2\)가 가질 수 있는 값은 \(1/200 = 0.005\)의 배수뿐이다. 계급폭을 그 격자에 맞추지 않으면 막대가 들쭉날쭉해진다.

(2) 그려서 확인한다.

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

rng = np.random.default_rng(1)

n1 = n2 = 200
p1, p2 = 0.50, 0.40
mu = p1 - p2
sd = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)

B = 200_000
x1 = rng.binomial(n1, p1, B)
x2 = rng.binomial(n2, p2, B)
d = x1 / n1 - x2 / n2

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

# (좌) 차 자체의 분포. n1 = n2 이므로 가능한 값은 1/200 = 0.005 의 배수다.
# 계급폭을 그 격자에 맞추어야 막대가 들쭉날쭉해지지 않는다.
ax = axes[0]
ax.hist(d, bins=np.arange(d.min() - 0.0025, d.max() + 0.0075, 0.005),
        density=True, alpha=0.5, edgecolor="white",
        label=r"simulated $\hat p_1-\hat p_2$")
g = np.linspace(mu - 4.5 * sd, mu + 4.5 * sd, 400)
ax.plot(g, stats.norm(mu, sd).pdf(g), "--r", lw=2, label="normal PDF")
ax.axvline(mu, color="k", lw=1, ls=":")
ax.set_xlabel(r"$\hat p_1 - \hat p_2$")
ax.set_ylabel("Density")
ax.set_title(r"$n_1=n_2=200$,  $p_1=0.50$,  $p_2=0.40$")
ax.legend(fontsize=8)

# (우) 참 표준오차로 표준화한 통계량. 계급폭 0.101 = 0.005 / SE 는
# 격자 간격을 z 눈금으로 옮긴 값이다.
ax = axes[1]
ax.hist((d - mu) / sd, bins=np.arange(-4.5, 4.55, 0.101), density=True,
        alpha=0.5, edgecolor="white", label="standardized statistic")
g = np.linspace(-4.5, 4.5, 400)
ax.plot(g, stats.norm.pdf(g), "--r", lw=2, label=r"$N(0,1)$ PDF")
for c in (-1.96, 1.96):
    ax.axvline(c, color="gray", lw=1, ls=":")
ax.set_xlabel(r"$\left[(\hat p_1-\hat p_2)-(p_1-p_2)\right]/\mathrm{SE}$")
ax.set_xlim(-4.5, 4.5)
ax.set_title("standardized:  dotted lines at $\\pm 1.96$")
ax.legend(fontsize=8)

plt.tight_layout()
plt.show()

히스토그램과 빨간 곡선이 사실상 구분되지 않는다. 왼쪽 그림에서 중심이 \(0.10\)에 있고, 오른쪽 그림에서 표준화한 통계량이 \(N(0,1)\)과 겹친다. 오른쪽 계급폭 \(0.101\)은 우연한 값이 아니다. 가능한 \(\hat p_1 - \hat p_2\) 값의 간격 \(0.005\)를 표준오차 \(0.04950\)으로 나눈 값이며, 이 통계량이 실은 격자 위의 이산확률변수라는 사실을 계급폭이 그대로 반영한다.

예측대로 어긋남이 보이지 않는다. 히스토그램과 빨간 곡선이 사실상 구분되지 않는다. 왼쪽 그림에서 중심이 \(0.10\)에 있고, 오른쪽 그림에서 표준화한 통계량이 \(N(0,1)\)과 겹친다. 왜도 \(-0.0099\)와 초과첨도 \(-0.0048\)이 그리는 차이는 선 두께 안에 들어간다.

오른쪽 계급폭 \(0.101\)은 우연한 값이 아니다. 가능한 \(\hat p_1-\hat p_2\) 값의 간격 \(0.005\)를 표준오차 \(0.049497\)로 나눈 값이다.

곧 \(z\) 눈금으로 옮긴 격자 간격이다. 계급 하나가 격자점 하나를 담도록 맞춘 것이고, 이 수가 계급폭으로 쓰였다는 사실 자체가 이 통계량이 연속이 아니라 격자 위의 이산확률변수임을 말해 준다.

그래서 이 그림이 보이는 것과 가리는 것이 갈린다. 봉우리의 모양은 보이고 — 그리고 맞고 — 격자 간격은 계급폭에 숨어 있다. \(z\) 눈금으로 \(0.101\)이면 꼬리확률을 잴 때 임계값이 격자점 어디에 떨어지느냐에 따라 덩어리 하나만큼 답이 달라질 수 있다는 뜻이다. 그 크기가 얼마인지는 보기 3에서 수로 잰다.

(1) 둘은 정확하고 둘은 근사다.

평균과 표준편차는 정확하다. \(E[\hat p_i] = p_i\)와 \(\operatorname{Var}(\hat p_i) = p_iq_i/n_i\)는 이항분포의 정의에서 바로 나오는 등식이고, 두 표본이 독립이므로 차의 평균과 분산도 등식이다. 여기에는 정규근사가 한 번도 쓰이지 않는다. 그러므로 전수 계산은 \(0.1000\)과 \(0.04950\)을 자릿수가 허락하는 한 그대로 되돌려 주어야 한다.

포함률은 근사다. 근사가 들어간 곳이 둘이다. 하나는 모양(\(Z\)가 정규라는 것), 다른 하나는 분모(\(p_i\) 자리에 \(\hat p_i\)를 끼워 넣은 것)다.

평균 폭은 근사이고, 방향과 크기까지 손으로 나온다. 구간의 폭은 \(2z\widehat{\text{SE}}\)이고 \(\widehat{\text{SE}} = \sqrt{\hat V}\)인데, \(\hat V\)는 \(V\)의 불편추정량이 아니다.

이므로 \(E[\hat V] = \left(1-\frac1n\right)V\)로 \(1/n\)만큼 작게 겨눈다. 여기에 제곱근이 위로 오목하다는 옌센의 몫이 한 번 더 아래로 더해진다. \(\hat V\)의 변동계수를 \(\mathrm{cv}\)라 하면 델타법이 \(E[\sqrt{\hat V}] \approx \sqrt{E[\hat V]}\,(1-\mathrm{cv}^2/8)\)을 주므로

이다(\(\mathrm{cv} = 0.0150\)). 두 몫의 크기가 전혀 다르다. 앞의 \(-0.2503\%\)가 거의 전부이고 옌센은 \(-0.0028\%\)에 그친다. 따라서 평균 폭은

로, \(2z\,\text{SE} = 0.1940\)보다 조금 작아야 한다.

(2) 전수 계산.

\(n_1 = n_2 = 200\)이면 \((X_1, X_2)\)의 조합이 \(201 \times 201 = 40401\)개뿐이다. 모의실험을 돌릴 필요 없이 모든 조합에 이항 확률을 곱해 정확한 값을 얻을 수 있다.

import numpy as np
from scipy import stats

n1 = n2 = 200
p1, p2 = 0.50, 0.40
mu = p1 - p2
sd = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
z975 = stats.norm.ppf(0.975)

k1 = np.arange(n1 + 1)
k2 = np.arange(n2 + 1)

# W[i, j] = P(X1 = i) * P(X2 = j). 두 표본이 독립이므로 곱이 된다.
W = np.outer(stats.binom.pmf(k1, n1, p1), stats.binom.pmf(k2, n2, p2))

X1, X2 = np.meshgrid(k1, k2, indexing="ij")
ph1, ph2 = X1 / n1, X2 / n2
d = ph1 - ph2

# 비합동 표준오차를 쓴 왈드 구간이 참값 mu 를 담는 칸을 골라낸다.
se = np.sqrt(ph1 * (1 - ph1) / n1 + ph2 * (1 - ph2) / n2)
inside = (d - z975 * se <= mu) & (mu <= d + z975 * se)

print(f"이론 평균 {mu:.4f},  실제 평균 {(W * d).sum():.4f}")
print(f"이론 표준오차 {sd:.5f},  실제 표준편차 {np.sqrt((W * (d - mu) ** 2).sum()):.5f}")
print(f"명목 95% 구간의 정확 포함률 {W[inside].sum():.4f}")
print(f"구간의 평균 폭 {(W * 2 * z975 * se).sum():.4f}")

출력:

이론 평균 0.1000,  실제 평균 0.1000
이론 표준오차 0.04950,  실제 표준편차 0.04950
명목 95% 구간의 정확 포함률 0.9504
구간의 평균 폭 0.1935

평균과 표준편차는 근사가 아니라 정확히 이론값과 같다. 애초에 근사가 들어간 곳은 모양뿐이기 때문이다. 포함률 \(0.9504\)는 명목 \(0.95\)와 0.04퍼센트포인트밖에 차이 나지 않는다.

(1)의 가르기가 그대로 확인된다. 평균 \(0.1000\)과 표준편차 \(0.04950\)은 이론값과 찍어 낸 자리까지 같다. 애초에 근사가 없었으므로 당연한 결과이며, 이 두 줄은 전수 계산의 검산 노릇을 한다.

포함률 \(0.9504\)는 명목에서 \(0.0004\) 떨어져 있다. 모의실험이 아니라 \(201 \times 201 = 40401\)칸을 모두 센 값이므로 여기에 몬테카를로 오차는 없다. \(0.0004\)는 전부 근사의 몫이고, 기대도수가 \(80\) 이상인 이 설정에서 왈드 구간이 실용적으로 완벽하다는 뜻이다.

평균 폭이 예측한 \(0.1935\) 그대로다. \(2z\,\text{SE} = 0.1940\)보다 작고, 전수 계산으로 재면 \(E[\widehat{\text{SE}}] = 0.049372\)로 참값 \(0.049497\)의 \(0.99747\)배다. (1)에서 손으로 낸 \(0.99747\)과 소수 다섯째 자리까지 같다. \(\hat p\hat q\)의 \((1-1/n)\) 편향이 거의 전부이고 옌센의 몫은 백분의 일에 그친다.

구간이 평균적으로 살짝 짧은데도 포함률이 모자라지 않는 까닭은 짧아지는 표본과 중심이 빗나가는 표본이 서로 다르기 때문이다. \(\widehat{\text{SE}}\)는 \(\hat p_i\)가 \(0.5\)에서 멀어질 때 작아지는데, 그런 표본은 \(\hat p_i\)가 참값에서 크게 벗어난 표본이 아니라 그저 한쪽으로 치우친 표본이다. 두 어긋남이 같은 방향으로 겹치지 않아 상쇄가 일어난다. 다음 쪽에서 \(n\)을 줄이면 이 상쇄가 깨지고 왈드 구간의 포함률이 무너진다.

(1) 어긋남은 모양이 아니라 격자에서 온다. 보기 1에서 \(\hat p_1-\hat p_2\)의 왜도가 \(-0.0099\), 초과첨도가 \(-0.0048\)이었다. 모양으로 생길 수 있는 오차는 이 크기를 넘지 못한다. 남는 범인은 이 통계량이 \(0.005\) 간격의 격자 위에만 산다는 사실이다.

격자 위의 확률은 덩어리로 놓여 있으므로, 임계값이 어느 격자점의 바로 옆에 떨어지면 덩어리 하나가 통째로 안쪽이나 바깥쪽으로 쏠린다. 세 지점의 경계를 격자와 맞춰 본다.

격자점은 \(0.005\)의 배수다. \(c = 1.96\)의 \(0.19702\)는 \(0.195\)에서 칸 폭의 \(40\%\) 들어간 자리이고 \(c = 2.576\)의 \(0.22751\)도 \(0.225\)에서 \(50\%\) 들어간 한가운데다. 둘 다 격자점에서 멀다. 그런데 \(c = 1\)의 \(0.14950\)은 격자점 \(0.150\)의 왼쪽 \(0.00050\) 자리, 곧 칸 폭 \(0.005\)의 10% 지점이다. \(c = 1\)에서 가장 크게 어긋나야 한다.

크기를 어림하자. 격자점 \(0.150\)에 놓인 덩어리가 \(P(\hat p_1-\hat p_2 = 0.150) = 0.0243\)인데, 경계가 그 왼쪽에 있으므로 이 덩어리가 전부 꼬리로 들어간다. 덩어리를 칸 폭 \(0.005\)에 고르게 펴 놓은 "연속판"에서라면 경계 오른쪽으로 넘어갈 몫은 \((0.1525-0.1495)/0.005 = 60\%\)뿐이므로, 나머지 \(40\%\)가 초과분이다.

\(P(|Z|>1)\)이 표준정규보다 약 \(0.019\) 클 것으로 예측된다.

(2) 전수 계산.

포함률은 하나의 숫자이므로 안쪽과 꼬리의 오차가 서로를 가려 줄 수 있다. 꼬리확률을 여러 지점에서 따로 본다.

zstat = (d - mu) / sd

print("표준화한 통계량의 꼬리확률")
for c in (1.000, 1.960, 2.576):
    exact = W[np.abs(zstat) > c].sum()
    print(f"  P(|Z| > {c:.3f}) = {exact:.4f}   표준정규 {2 * (1 - stats.norm.cdf(c)):.4f}")

출력:

표준화한 통계량의 꼬리확률
  P(|Z| > 1.000) = 0.3372   표준정규 0.3173
  P(|Z| > 1.960) = 0.0487   표준정규 0.0500
  P(|Z| > 2.576) = 0.0099   표준정규 0.0100

실무에서 쓰는 두 지점(\(1.96\), \(2.576\))은 소수 셋째 자리까지 맞는다. 그런데 \(c = 1\)에서는 \(0.3372\) 대 \(0.3173\)으로 2퍼센트포인트나 어긋난다.

원인은 정규근사의 실패가 아니라 이산성이다. \(c = 1\)에 해당하는 경계 \(0.10 + 0.04950 = 0.14950\)은 격자점 \(0.150 = 30/200\)의 바로 왼쪽(\(0.00050\) 차이)에 있다. 그래서 \(P(\hat p_1 - \hat p_2 = 0.150) = 0.0243\)이라는 덩어리 하나가 통째로 꼬리에 들어간다. 이 덩어리를 폭 \(0.005\)에 고르게 펴 놓았다면 경계 오른쪽으로 넘어갈 몫은 60%뿐이므로 약 \(0.0243 \times 0.4 = 0.0097\)이 초과된다. 양쪽 꼬리를 합치면 \(0.0194\)이고, 관측된 차이 \(0.3372 - 0.3173 = 0.0199\)와 같은 크기다.

기대도수가 충분히 커도 이산성은 이렇게 남아 있다. 다음 페이지에서는 이 덩어리가 통계량 전체를 지배하게 된다.

예측 \(0.0194\)에 실제 \(0.3372 - 0.3173 = 0.0199\)다. 격자 덩어리 하나로 어긋남의 거의 전부가 설명된다. 나머지 두 지점은 소수 셋째 자리까지 맞는다.

실무에서 쓰는 두 지점이 맞은 것은 운이 좋아서다. 근사가 특별히 좋아서가 아니라 경계가 격자점에서 멀리 떨어진 자리에 놓였기 때문이다. \(n\)이나 \(p\)를 조금만 바꾸면 \(1.96\)의 경계도 격자점 옆으로 옮겨 가고, 그때는 \(c=1\)에서 본 것과 같은 크기의 어긋남이 거기서 나타난다.

이것이 포함률 하나만 보면 안 되는 이유다. 보기 2가 준 \(0.9504\)는 "안쪽 전체"를 재므로 양쪽 꼬리의 어긋남이 서로를 지우거나 묻힌다. 꼬리를 여러 지점에서 따로 보면 숨어 있던 \(2\)퍼센트포인트가 드러난다.

어긋남은 작아지지만 없어지지 않는다. 격자 간격은 \(1/n\)이고 표준오차는 \(1/\sqrt n\)에 비례하므로, \(z\) 눈금에서 잰 간격은 \(1/\sqrt n\)으로 줄어든다. \(n\)을 네 배로 키워야 절반이 된다. 0으로 가기는 하되 느리게 가며, 연속성 수정이 바로 이 몫을 손으로 메우는 장치다. 다음 쪽에서는 이 덩어리가 작아지기는커녕 통계량 전체를 지배한다.

(1) 비 하나가 답을 정한다. 구간의 반폭은 \(z_{0.975}\widehat{\text{SE}}\)이므로, 쓰는 표준오차가 참값의 \(r\)배이면 구간은 참 표준오차 눈금으로 \(\pm 1.96r\)만큼 뻗는다. 표본이 커서 \(\hat p_1-\hat p_2\)가 거의 정규라면

이다. 비합동은 \(r \approx 1\)이라 \(0.95\)를 준다. 합동은 \(p_1 \ne p_2\)일 때 \(r\)이 \(1\)에서 벗어난다. 두 표본크기가 같으면 \(\hat p \to (p_1+p_2)/2 \equiv \bar p\)이므로

이다. \(p_1 = \bar p + \delta/2\), \(p_2 = \bar p - \delta/2\)(\(\delta = p_1-p_2\))를 넣으면

가 정확히 성립하므로

임을 알 수 있다. 합동 표준오차는 언제나 참값 이상이고, 두 비율이 멀어질수록 커진다. 네 조합에 넣으면 이렇다.

\((0.9, 0.1)\)에서 \(r = 5/3\)이다. 합동 표준오차 \(\sqrt{0.25 \times 2/200} = 0.05\)는 두 비율이 둘 다 \(0.5\)인 양 셈한 값이고, 실제로는 \(0.9\)와 \(0.1\)이라 각 비율의 분산 \(pq = 0.09\)가 \(0.25\)의 삼분의 일 남짓이다. 구간이 참값 눈금으로 \(\pm 3.27\)까지 뻗으니 거의 아무것도 놓치지 않는다.

검정 쪽은 사정이 다르다. \(H_0\)가 참이면 \(p_1 = p_2\)이므로 위 식에서 \((p_1-p_2)^2 = 0\), 곧 \(r = 1\)이다. 합동과 비합동이 겨누는 값이 같아진다. 그래서 네 줄이 모두 \(0.05\) 둘레에 있어야 한다.

(2) 전수 계산.

\(H_0: p_1 = p_2\) 검정의 실제 1종오류율과, 신뢰구간에 합동 표준오차를 잘못 썼을 때의 포함률을 함께 계산한다.

ppool = (X1 + X2) / (n1 + n2)
se_pool = np.sqrt(ppool * (1 - ppool) * (1 / n1 + 1 / n2))

# 표준오차가 0 인 칸(두 표본이 모두 전원 성공 또는 전원 실패)은
# 차도 0 이므로 기각하지 않는 것으로 둔다.
z_pool = np.divide(d, se_pool, out=np.zeros_like(d), where=se_pool > 0)
z_unpool = np.divide(d, se, out=np.zeros_like(d), where=se > 0)

print("검정: 명목 유의수준 5% 의 실제 1종오류율")
for p in (0.05, 0.20, 0.40, 0.50):
    W0 = np.outer(stats.binom.pmf(k1, n1, p), stats.binom.pmf(k2, n2, p))
    print(f"  p1 = p2 = {p:.2f}:  합동 {W0[np.abs(z_pool) > z975].sum():.4f}"
          f"   비합동 {W0[np.abs(z_unpool) > z975].sum():.4f}")

print("신뢰구간: 명목 95% 의 실제 포함률")
for a, b in [(0.50, 0.40), (0.60, 0.40), (0.80, 0.20), (0.90, 0.10)]:
    Wab = np.outer(stats.binom.pmf(k1, n1, a), stats.binom.pmf(k2, n2, b))
    t = a - b
    ok = Wab[(d - z975 * se <= t) & (t <= d + z975 * se)].sum()
    bad = Wab[(d - z975 * se_pool <= t) & (t <= d + z975 * se_pool)].sum()
    print(f"  p1 = {a:.2f}, p2 = {b:.2f}:  비합동 {ok:.4f}   합동 {bad:.4f}")

출력:

검정: 명목 유의수준 5% 의 실제 1종오류율
  p1 = p2 = 0.05:  합동 0.0481   비합동 0.0489
  p1 = p2 = 0.20:  합동 0.0508   비합동 0.0509
  p1 = p2 = 0.40:  합동 0.0482   비합동 0.0497
  p1 = p2 = 0.50:  합동 0.0510   비합동 0.0510
신뢰구간: 명목 95% 의 실제 포함률
  p1 = 0.50, p2 = 0.40:  비합동 0.9504   합동 0.9513
  p1 = 0.60, p2 = 0.40:  비합동 0.9474   합동 0.9536
  p1 = 0.80, p2 = 0.20:  비합동 0.9460   합동 0.9853
  p1 = 0.90, p2 = 0.10:  비합동 0.9494   합동 0.9987

검정에서는 두 방식이 사실상 같다. \(H_0\)가 참이면 \(\hat p_1\)과 \(\hat p_2\)가 모두 공통값 \(p\)에 가까우므로 합동과 비합동이 같은 값으로 수렴하기 때문이다. 합동을 권하는 이유는 그 수렴이 더 빠르다는 것뿐이다.

신뢰구간에서는 사정이 다르다. 비합동 구간은 어느 조합에서든 \(0.95\) 근처를 지키지만, 합동 구간은 두 비율이 멀어질수록 \(0.9853\), \(0.9987\)로 치솟는다. 구간이 명목보다 훨씬 넓다는 뜻이다. 필요 이상으로 넓은 구간은 안전해 보이지만 아무것도 말해 주지 않는 구간이기도 하다.

예측이 네 자리까지 맞는다.

가장 크게 어긋난 줄이 \((0.60,0.40)\)의 \(0.9545\) 대 \(0.9536\)으로 \(0.0009\)이고, 나머지 세 줄은 \(0.0004\) 안이다. 표준오차의 비 하나로 포함률이 정해진다는 (1)의 주장이 그대로 섰다.

검정 열도 예측대로다. 네 줄이 \(0.0481\)에서 \(0.0510\) 사이이고 합동과 비합동의 차이가 최대 \(0.0015\)다. \(H_0\)가 참이면 \(r = 1\)이므로 두 방식이 같은 것을 겨누고, 남는 차이는 유한표본에서 \(\hat p\)와 \((\hat p_1, \hat p_2)\)가 서로 다르게 흔들리는 몫과 이산성뿐이다. 합동을 권하는 까닭도 거기에 있다. 두 표본을 합쳐 \(400\)개로 \(p\)를 추정하므로 같은 값을 더 안정적으로 겨눈다.

그러니 둘의 자리는 이렇게 갈린다. 합동은 \(p_1 = p_2\)를 가정하고 계산한 흔들림이다. 그 가정이 전제인 검정에서는 정당하고 더 낫다. 신뢰구간에는 그 가정이 없으므로, 합동을 쓰는 것은 묻지 않기로 한 것을 답으로 미리 써넣는 일이다. 대가는 \((0.9,0.1)\)에서 포함률 \(0.9987\), 곧 명목 \(95\%\)를 내걸고 사실상 \(99.9\%\) 구간을 보고하는 것이며, 그만큼 넓어 아무것도 배제하지 못한다.

분산을 더한다.

기대도수는

으로 네 값이 모두 5보다 훨씬 크다. 조건이 넉넉히 충족된다.

눈여겨볼 점. \(p(1-p)\)는 \(p = 0.5\)에서 최대 \(0.25\)이고 \(p = 0.4\)에서 \(0.24\)로 거의 차이가 없다. \(p\)가 \(0.3\)과 \(0.7\) 사이이면 \(0.21\)에서 \(0.25\) 사이에 머물므로 표준오차가 \(p\)에 거의 의존하지 않는다. 표본크기 계산에서 \(p\)를 모를 때 \(0.5\)를 쓰는 관행이 여기서 나온다.

\(X_i \sim \text{Binomial}(n_i, p_i)\)이므로 \(E[X_i] = n_ip_i\), \(\text{Var}(X_i) = n_ip_iq_i\)이고

이다. 기댓값은 선형이고 독립이면 공분산이 0이므로

이다. 어디에도 근사가 없으며 \(n_i = 1\)이어도 성립한다. \(\square\)

근사는 모양에만 쓰인다. \(\hat p_1 - \hat p_2\)의 참 분포는 두 이항분포의 차를 척도변환한 것으로, \(n_1 = n_2 = n\)이면 \(\{-1, -1+1/n, \dots, 1\}\) 위의 이산분포다. 이 이산분포를 연속분포로 바꾸어 읽는 것이 근사의 전부다. 보기 2에서 평균과 표준편차는 소수점까지 맞고 포함률만 \(0.9504\)로 어긋난 이유가 이것이며, 어긋나는 것은 언제나 확률을 읽는 방식이다.

검정. 검정은 "\(H_0\)가 참이라면 이만한 차가 얼마나 흔한가"를 묻는다. 확률을 계산하는 무대가 \(H_0\) 아래이므로 표준오차도 \(H_0\) 아래에서 재야 한다. \(H_0: p_1 = p_2 = p\)라면 공통값 \(p\)의 최량 추정값은 두 표본을 합친 \(\hat p = (X_1+X_2)/(n_1+n_2)\)이며, 자료를 두 배 쓰므로 표준오차가 더 안정적이다.

신뢰구간. 신뢰구간은 \(p_1 - p_2\)가 무엇이든 \(95\%\)를 담아야 하는 물건이므로 어떤 값도 가정할 수 없다. 특히 \(p_1 - p_2 = 0\)은 신뢰구간이 판정해야 할 대상이니, 그것을 재료로 쓰면 순환논법이 된다.

무엇이 잘못되는가. \(n_1 = n_2\)이면 합동 표준오차는 \(\bar p = (p_1+p_2)/2\)를 써서 \(\sqrt{2\bar p\bar q/n}\)에 수렴한다. \(p(1-p)\)가 오목함수이므로 옌센 부등식에 따라

이고, 등호는 \(p_1 = p_2\)에서만 성립한다. 합동 표준오차는 언제나 참 표준오차보다 크거나 같다. 따라서 합동으로 만든 구간은 항상 지나치게 넓고 두 비율이 멀어질수록 더 넓어진다(\(p_1=0.80\), \(p_2=0.20\)에서 \(\bar p\bar q = 0.25\) 대 \(0.16\), 포함률 \(0.9853\)).

넓은 구간은 무해하지 않다. 대응하는 검정의 검정력이 낮아지고, 보고된 불확실성이 실제보다 커져 자료가 준 정보를 버리게 된다.

\(n_1 = n_2 = 200\)이면 \(\hat p_1 - \hat p_2 = (X_1 - X_2)/200\)이므로 가능한 값이 \(0.005\) 간격의 격자 위에 있고, 표준오차 \(0.04950\)으로 나누면 격자 간격이 \(z\) 눈금으로 \(0.005/0.04950 = 0.1010\)이 된다. 경계 \(0.14950\)은 격자점 \(0.150\)의 왼쪽, 간격의 10%밖에 안 되는 거리에 놓인다. 그래서 확률 \(0.0243\)짜리 덩어리가 통째로 꼬리로 넘어가고, 보기 3에서 본 \(0.0194\)의 초과가 생긴다.

줄이는 방법. 첫째, 연속성 수정으로 경계를 격자 반칸만큼 밀어 덩어리를 나눈다. 다만 평균적으로 지나치게 보수적이라는 비판을 받는다. 둘째, 표본을 키운다. \(z\) 눈금의 격자 간격은 \(1/\sqrt n\)에 비례해 줄어들므로 \(n = 200\)의 \(0.101\)이 \(n = 2000\)에서는 \(0.032\)가 된다. 셋째, 이항 확률을 그대로 쓴다. 보기 2와 4가 한 일이다.

오차는 근사의 질이 아니라 임계값의 위치에 달렸다. 임계값이 격자점 중간에 떨어지면 오차가 거의 0이 되고, 격자점 옆에 떨어지면 최대가 된다. 그래서 \(c = 1.96\)에서는 잘 맞고 \(c = 1\)에서는 어긋난다. 이 현상이 극단으로 치달으면 다음 페이지의 톱니 그래프가 된다.

보기 4의 코드에서 \(p\)만 바꾸면 된다. 전수 계산이므로 모의오차가 없다.

첫 줄이 보기 4의 출력과 일치하는 것이 검산이다.

문턱값이 그대로 확인된다. \(np = 10\)에서는 \(0.048\)로 명목 \(0.05\)와 사실상 같다. \(5\)로 내려가면 \(0.046\), \(2\)에서 \(0.041\), \(1\)에서 \(0.014\)다. 모양만 눈으로 볼 때는 5로 족하지만 오류율을 명목에 맞추려면 10이 필요하다는 4.1절의 약속이 숫자로 나타난 것이다.

방향은 보수적이다. 명목보다 덜 기각하므로 1종오류만 보면 안전해 보인다. 그러나 안전한 것이 아니라 자기가 몇 %짜리 검정인지 모르는 것이다. 실제 수준이 \(0.014\)라면 기각역이 명목보다 훨씬 좁다는 뜻이고, 참으로 두 비율이 다를 때 그것을 잡아낼 힘도 그만큼 줄어든다.

원인은 정규근사의 실패라기보다 이산성이다. \(np\)가 작으면 \(\hat p_1 - \hat p_2\)가 실질적으로 취하는 값이 0 언저리의 몇 개로 줄어들어 달성 가능한 확률이 띄엄띄엄해진다. 그 띄엄띄엄한 값 가운데 명목 \(0.05\)에 맞아떨어지는 것이 있을 이유가 없다.

표의 두 열이 거의 같다는 점도 눈여겨볼 만하다. 귀무가설이 참인 자리에서는 합동과 비합동이 같은 값으로 모이므로, 어긋남은 두 표준오차의 선택이 아니라 이산성에서 온다.

격자. 통분하면

이고 \(3k_1 - 4k_2\)는 정수를 모두 만들 수 있으므로(\(3 \times (-1) - 4 \times (-1) = 1\)) 격자 간격이 \(1/600 = 0.001667\)이다. \(200\)과 \(150\)의 최소공배수가 \(600\)이기 때문이며, 등크기일 때의 \(1/200\)보다 세 배 촘촘하다.

표준오차는 \(\sqrt{0.25/200 + 0.24/150} = 0.05339\)이므로 \(z\) 눈금으로 옮긴 격자 간격은 \(0.001667/0.05339 = 0.0312\)다. 전수 계산으로 꼬리확률을 재면 다음과 같다.

표본을 늘리지 않고도 이산성이 줄었다. 오히려 \(n_2\)를 \(50\)명 줄였는데 어긋남이 \(0.0199\)에서 \(0.0068\)로 작아졌다. 격자가 촘촘해지면 어느 한 덩어리도 경계를 통째로 넘나들지 못하기 때문이다. 표본을 열 배로 키워 얻는 격자 간격 \(0.032\)를 \(n_2\)를 줄이는 것만으로 얻은 셈이다.

그렇다고 불균형 설계를 권하는 것은 아니다. 이산성은 근사의 한 축일 뿐이다. \(n_2\)를 줄이면 표준오차가 \(0.04950\)에서 \(0.05339\)로 커지고 기대도수도 함께 줄어드니, 표본분포의 폭은 넓어지고 조건은 빠듯해진다. 여기서 얻을 교훈은 설계 지침이 아니라 원인의 정체다. 이산성에서 오는 오차는 근사의 질이 아니라 격자의 성김과 임계값의 위치가 정한다.

차에 대한 연속성 수정은 \(|\hat p_1-\hat p_2|\)에서 \(\frac12(1/n_1+1/n_2)\)를 빼는 형태다.

import numpy as np

rng = np.random.default_rng(0)
z = 1.959964
n1 = n2 = 25
REP = 200_000
for p in (0.3, 0.5):
    a = rng.binomial(n1, p, REP) / n1
    b = rng.binomial(n2, p, REP) / n2
    pool = (a*n1 + b*n2) / (n1 + n2)
    se = np.sqrt(pool*(1-pool)*(1/n1 + 1/n2))
    ok = se > 0
    Z = np.zeros(REP); Z[ok] = (a[ok] - b[ok]) / se[ok]
    cc = np.zeros(REP)
    d = np.abs(a - b) - 0.5*(1/n1 + 1/n2)
    cc[ok] = np.maximum(d[ok], 0) / se[ok]
    print(f"  p={p}: 수정없음 {np.mean(np.abs(Z)>z):.4f}   "
          f"연속성수정 {np.mean(cc>z):.4f}")

출력:

  p=0.3: 수정없음 0.0536   연속성수정 0.0228
  p=0.5: 수정없음 0.0655   연속성수정 0.0335

연속성 수정이 지나치게 보수적이다. \(p=0.3\)에서 \(0.0536\)이던 오류율이 \(0.0228\)로 명목의 절반 이하로 떨어진다. 수정하지 않은 쪽이 오히려 \(0.05\)에 가깝다.

하나짜리 비율과 사정이 다른 이유는 수정량이 격자 한 칸을 통째로 밀기 때문이다. \(\hat p\) 하나에서는 격자 간격이 \(1/n\)이고 수정량 \(\frac{1}{2n}\)이 정확히 격자의 절반이라 잘 맞았다. 그런데 차의 격자는 등크기일 때 좁아지지 않는다. 연습문제 10에서 보듯 격자는 두 분모의 최소공배수가 정하므로, \(n_1=n_2=25\)이면 간격이 \(1/25 = 0.04\)로 한 비율일 때와 똑같다. 반면 수정량은 \(\frac12(1/n_1+1/n_2) = 0.04\)로 반 칸이 아니라 한 칸이다. 밀어야 할 거리의 두 배를 미는 셈이니 과한 보정이다.

실무 권고. 두 비율의 차에는 연속성 수정을 쓰지 않는 것이 일반적이다. 정확한 판정이 필요하면 수정 대신 피셔 정확검정이나 바너드 검정으로 가는 것이 낫다(다음 쪽). scipy.stats.chi2_contingency가 \(2\times2\)에서 예이츠 보정을 기본으로 걸어 두는 것도 같은 이유로 논란이 있으며, 많은 통계학자가 correction=False를 권한다.

일반 교훈. 한 상황에서 잘 듣는 보정이 비슷해 보이는 다른 상황에서도 잘 들으리라는 보장이 없다. 보정은 그것이 유도된 조건 안에서만 타당하다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
n1 = n2 = 50
p = 0.1
REP = 300_000
k1 = rng.binomial(n1, p, REP)
k2 = rng.binomial(n2, p, REP)

d = k1/n1 - k2/n2
ok = (k1 > 0) & (k1 < n1) & (k2 > 0) & (k2 < n2)
lor = np.log((k1[ok]*(n2-k2[ok])) / ((n1-k1[ok])*k2[ok]))

print(f"  차      : 왜도 {stats.skew(d):+.4f}   첨도 {stats.kurtosis(d, fisher=False):.4f}")
print(f"  로그오즈비: 왜도 {stats.skew(lor):+.4f}   첨도 {stats.kurtosis(lor, fisher=False):.4f}")
print(f"  (로그오즈비 계산 가능한 비율: {ok.mean():.4f})")

출력:

  차      : 왜도 +0.0024   첨도 3.0380
  로그오즈비: 왜도 +0.0057   첨도 3.3970
  (로그오즈비 계산 가능한 비율: 0.9897)

\(p_1 = p_2\)이면 둘 다 거의 대칭이다. 왜도가 \(0.006\) 이하로, 두 집단이 대칭적으로 들어가므로 차의 왜도가 상쇄된다. 이 쪽 연습문제 2에서 본 성질이다.

그런데 첨도는 로그 오즈비 쪽이 더 크다(\(3.40\) 대 \(3.04\)). 분모에 들어가는 \(k_2\)가 작을 때 로그 오즈비가 크게 튀기 때문이며, \(p=0.1\)이라 \(k_2\)가 한 자릿수인 경우가 흔하다. 로그 변환이 왜도는 잡아도 꼬리까지 얇게 만들지는 않는다.

그리고 로그 오즈비에는 계산 자체의 문제가 있다. \(k_1\)이나 \(k_2\)가 \(0\)이면 오즈비가 \(0\) 또는 \(\infty\)가 되어 로그를 취할 수 없다. 위에서 계산 가능한 비율이 \(0.9897\)로 \(100\%\)가 아닌 것이 그 때문이며, \(p\)가 더 작거나 \(n\)이 작으면 이 비율이 크게 떨어진다.

\(p_1 \ne p_2\)일 때 차이가 드러난다. 차는 참값이 \(0\)에서 멀어질수록 경계 효과로 분포가 치우치지만, 로그 오즈비는 그렇지 않다. 그래서 연구 간 통합이나 회귀에서는 로그 오즈비를 쓰고, 단일 연구의 보고에는 차를 쓰는 관행이 자리 잡았다.

\(k=0\) 문제의 표준 대처는 \(0.5\)씩 더하는 것이다(하네스–헐 보정). 5.5절 희귀사건 문서 연습문제 10에서 쓴 방법이며, 편향을 조금 넣는 대신 계산 가능성을 얻는다.

import numpy as np

rng = np.random.default_rng(2)
z = 1.959964
n1 = n2 = 100
REP = 200_000

def rates(p1, p2):
    k1 = rng.binomial(n1, p1, REP); k2 = rng.binomial(n2, p2, REP)
    a, b = k1/n1, k2/n2
    pool = (k1 + k2) / (n1 + n2)
    se_p = np.sqrt(pool*(1-pool)*(1/n1 + 1/n2))
    se_u = np.sqrt(a*(1-a)/n1 + b*(1-b)/n2)
    rp = np.mean(np.abs(a-b) > z*np.where(se_p > 0, se_p, np.inf))
    ru = np.mean(np.abs(a-b) > z*np.where(se_u > 0, se_u, np.inf))
    return rp, ru

print(f"{'p1,p2':>12}{'합동':>10}{'비합동':>10}")
for p1, p2 in [(0.3,0.3), (0.5,0.5), (0.3,0.45), (0.1,0.2)]:
    rp, ru = rates(p1, p2)
    print(f"{f'{p1},{p2}':>12}{rp:>10.4f}{ru:>10.4f}")

출력:

       p1,p2        합동       비합동
     0.3,0.3    0.0508    0.0527
     0.5,0.5    0.0560    0.0560
    0.3,0.45    0.5947    0.5982
     0.1,0.2    0.5184    0.5184

차이가 생각보다 작다. \((0.3,0.3)\)에서 합동 \(0.0508\), 비합동 \(0.0527\)로 합동이 조금 낫지만, \((0.5,0.5)\)에서는 두 값이 똑같이 \(0.0560\)이다. \(n_1=n_2\)이고 \(p=0.5\)이면 두 표준오차가 거의 같아져 판정이 갈리는 표가 사실상 없기 때문이다.

이유는 귀무가설을 실제로 쓰기 때문이다. \(H_0: p_1=p_2\)가 참이면 두 표본이 같은 모집단에서 온 것이므로, 둘을 합쳐 \(p\)를 추정하는 편이 자유도를 더 쓰는 것이라 정밀하다. 5.5절 윌슨 구간에서 본 "추정값 대신 가설값으로 표준오차를 계산하라"는 점수 검정의 원리와 같다.

대립가설 아래에서도 거의 같다. 검정력이 \(0.5947\) 대 \(0.5982\)이고, \((0.1,0.2)\)에서는 완전히 같다. 실무적으로는 어느 쪽을 써도 결과가 거의 달라지지 않는다.

그래서 관행이 갈린다.

이 불일치가 혼란을 부른다. 검정은 유의한데 신뢰구간은 \(0\)을 포함하는(또는 그 반대의) 일이 드물게 생긴다. 두 절차가 다른 표준오차를 쓰기 때문이며, 논리적 모순이 아니라 설계의 차이다. 보고할 때는 어느 쪽을 썼는지 밝히는 것이 좋다.

격자는 두 분모의 최소공배수가 정한다. \(\hat p_1 - \hat p_2 = k_1/n_1 - k_2/n_2\)이므로 가능한 값은

의 형태다. \(n_1=n_2=200\)이면 분모가 \(200\)으로 약분되어 간격 \(1/200 = 0.005\)의 균일 격자다.

from math import gcd

# 부동소수점으로 세면 수학적으로 같은 값이 미세하게 달라져 과대계수된다.
# 정수 산술로 (k1*n2 - k2*n1)/(n1*n2) 를 약분해 센다.
for n1, n2 in [(200,200), (200,150), (200,151)]:
    g = gcd(n1, n2)
    L = n1 * n2 // g
    vals = {(k1*n2 - k2*n1) // g
            for k1 in range(n1+1) for k2 in range(n2+1)}
    print(f"  n1={n1}, n2={n2}: gcd={g:>3}  간격 1/{L} = {1/L:.6f}  "
          f"서로 다른 값 {len(vals):,}개")

출력:

  n1=200, n2=200: gcd=200  간격 1/200 = 0.005000  서로 다른 값 401개
  n1=200, n2=150: gcd= 50  간격 1/600 = 0.001667  서로 다른 값 1,195개
  n1=200, n2=151: gcd=  1  간격 1/30200 = 0.000033  서로 다른 값 30,551개

(부동소수점 np.unique로 세면 수학적으로 같은 값이 미세하게 달라져 \(200,200\)에서도 \(1{,}293\)개가 나온다. 이산 구조를 다룰 때는 정수 산술로 세야 한다.)

\(n_2\)를 \(150\)에서 \(151\)로 하나만 바꾸면 격자가 \(50\)배 조밀해진다. \(\gcd(200,150)=50\)인데 \(\gcd(200,151)=1\)이기 때문이다. 표본크기의 정수론적 성질이 격자를 좌우한다.

그러나 격자가 조밀하다고 포함률이 좋아지지는 않는다. 값의 개수는 늘어나도 확률질량은 여전히 몇 개 값에 몰려 있다. \(30{,}551\)개 값 중 대부분은 확률이 사실상 \(0\)이고, 실제로 일어나는 것은 \(\hat p_1\)과 \(\hat p_2\)가 각자의 평균 근처일 때의 조합뿐이다.

포함률의 톱니는 사라지지 않는다. 5.5절 포함률 문서에서 본 대로 톱니는 "어떤 \(k\)의 구간이 참값을 담느냐 마느냐"가 갑자기 바뀌는 데서 생기는데, 그 \(k\)들이 가진 확률질량이 문제이지 격자의 조밀함이 문제가 아니다.

실무적 함의는 소박하다. 표본크기를 정할 때 \(\gcd\)를 신경 쓸 필요는 없다. 다만 "\(n\)을 조금 늘렸는데 결과가 이상하게 달라졌다"는 경험의 배경에 이런 이산 구조가 있다는 것은 알아 둘 만하다. 근본적인 해결은 격자를 조작하는 것이 아니라 이산성을 제대로 다루는 구간(윌슨, 뉴콤)을 쓰는 것이다.