베르누이분포¶

\(X \sim \text{Bernoulli}(0.02)\)이므로

이다.

표준편차 0.14가 평균 0.02보다 일곱 배 크지만 전혀 이상하지 않다. \(X\)는 0 아니면 1만 가지는데 평균은 0.02로 둘 사이 어디에도 없는 값이다. 한 번의 관측으로는 평균 근처의 값이 아예 나올 수 없는 분포이며, 이런 상황에서 변동계수 \(\text{SD}/E = \sqrt{(1-p)/p} = 7\)이 큰 것은 당연하다.

실무적 함의는 분명하다. 클릭률처럼 \(p\)가 작은 양을 추정하려면 표본이 아주 많이 필요하다. 방문자 \(n\)명의 클릭률 추정값 \(\hat p\)의 표준오차가 \(\sqrt{p(1-p)/n} = 0.14/\sqrt n\)이므로, 오차를 0.002 안에 넣으려면 \(n \approx 4{,}900\)명이 필요하다.

(1) 해석적으로. 거울상은 이름 바꾸기다. 본문 "무엇을 1이라 부를지는 우리가 정한다"에서 말한 대로, \(X \sim \text{Bernoulli}(p)\)이면 성공과 실패를 맞바꾼 \(Y = 1 - X\)는

이므로 \(Y \sim \text{Bernoulli}(1-p)\)다. 곧 \(\text{Bernoulli}(0.2)\)와 \(\text{Bernoulli}(0.8)\)은 같은 분포를 가로축만 뒤집어 본 것이다. 두 칸이 거울상인 데에는 그 이상의 뜻이 없다.

거울상이면 평균은 \(p \leftrightarrow 1-p\)로 바뀌고 분산 \(p(1-p)\)는 그대로이며, 왜도는 부호만 뒤집힌다. 성질 표의 식에 넣으면

이다. 세 경우를 계산한다.

첫째 줄과 셋째 줄이 분산은 같고 왜도는 부호만 다르다. 가운데 줄은 자기 자신의 거울상이므로(\(1 - 0.5 = 0.5\)) 왜도가 \(0\)일 수밖에 없다. 대칭인 분포의 왜도가 \(0\)이라는 것을 식을 보지 않고도 알 수 있는 경우다.

(2) 해석적으로. 최빈값은 두 확률 \(1-p\)와 \(p\) 가운데 큰 쪽이 놓인 자리다. 비로 쓰면

이고 이것이 \(1\)보다 크면 최빈값이 \(1\), 작으면 \(0\)이다. 비가 정확히 \(1\)이 되는 자리는

하나뿐이고, 그때는 \(P(X=0) = P(X=1) = 1/2\)이라 최빈값이 둘이다. 성질 표가 \(p = 1/2\)을 적지 않은 것은 빠뜨린 것이 아니라 그 자리에서 최빈값이 하나로 정해지지 않기 때문이다. 이 책에서 이산분포의 최빈값을 구할 때마다 되풀이되는 구조이고, 이항분포의 \((n+1)p \in \mathbb{Z}\), 포아송분포의 \(\lambda \in \mathbb{Z}\)가 모두 같은 동점이다. 베르누이는 그 가운데 가장 단순한 경우다.

수치적으로. 먼저 쪽의 그림을 그린다.

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

ps = [0.2, 0.5, 0.8]
x = np.array([0, 1])

fig, axes = plt.subplots(1, 3, figsize=(12, 3), sharey=True)
# 막대가 둘뿐인 분포다. 높이는 1-p 와 p 이고 합은 언제나 1이다.
# p가 커질수록 질량이 0에서 1로 옮겨 갈 뿐 모양이랄 것이 없다.
for ax, p in zip(axes, ps):
    ax.bar(x, stats.bernoulli(p).pmf(x), width=0.4, alpha=0.8)
    ax.set_title(f'Bernoulli(p={p})')
    ax.set_xticks(x)
    ax.set_xlabel('x')
    ax.spines[['top', 'right']].set_visible(False)
axes[0].set_ylabel('P(X = x)')
plt.tight_layout()
plt.show()

첫째 칸은 왼쪽 막대가 높고 셋째 칸은 오른쪽 막대가 높으며 두 칸의 막대 높이가 서로 맞바뀌어 있다. 가운데 칸은 두 막대가 같은 높이다. 분수로 확인한다.

from fractions import Fraction as F
from math import sqrt
import numpy as np
from scipy import stats

for p in (F(1, 5), F(1, 2), F(4, 5)):
    q = 1 - p
    ratio = p / q                      # P(X=1)/P(X=0)
    mode = "1" if ratio > 1 else ("0" if ratio < 1 else "0 과 1 (동점)")
    pm = stats.bernoulli(float(p)).pmf(np.array([0, 1]))
    print(f"p = {str(p):>3}:  P(0) = {str(q):>3}  P(1) = {str(p):>3}  "
          f"비 P(1)/P(0) = {str(ratio):>3}  최빈값 {mode}")
    print(f"          분산 {float(p*q):.4f}   왜도 {(1-2*float(p))/sqrt(float(p*q)):+.4f}"
          f"   부동소수점 pmf(1) - pmf(0) = {pm[1]-pm[0]:+.3e}")

출력:

p = 1/5:  P(0) = 4/5  P(1) = 1/5  비 P(1)/P(0) = 1/4  최빈값 0
          분산 0.1600   왜도 +1.5000   부동소수점 pmf(1) - pmf(0) = -6.000e-01
p = 1/2:  P(0) = 1/2  P(1) = 1/2  비 P(1)/P(0) =   1  최빈값 0 과 1 (동점)
          분산 0.2500   왜도 +0.0000   부동소수점 pmf(1) - pmf(0) = +1.110e-16
p = 4/5:  P(0) = 1/5  P(1) = 4/5  비 P(1)/P(0) =   4  최빈값 1
          분산 0.1600   왜도 -1.5000   부동소수점 pmf(1) - pmf(0) = +6.000e-01

(1)의 표가 그대로 나온다. \(p = 1/5\)과 \(p = 4/5\)에서 분산이 똑같이 \(0.16\)이고 왜도가 \(+1.5\)와 \(-1.5\)로 부호만 다르다. 차이 \(\text{pmf}(1) - \text{pmf}(0)\)도 \(-0.6\)과 \(+0.6\)으로 부호만 뒤집혀 거울상임을 그대로 보인다. 비 역시 \(1/4\)과 \(4\)로 서로 역수인데, 거울상이면 \(p/(1-p)\)가 역수가 되는 것이 당연하다.

\(p = 1/2\)의 마지막 줄을 보라. 두 확률이 수학적으로 똑같이 \(1/2\)인데도 차이가 \(0\)이 아니라 \(1.11 \times 10^{-16}\)으로 나왔다. \(1/2\)은 배정도 실수로 정확히 표현되는 수인데도 그렇다. SciPy 가 PMF를 정의 1의 식 \(p^x(1-p)^{1-x}\) 꼴로 계산하기 때문이고, 실제로 pmf(0, 0.5)가 \(0.4999999999999999\)를 돌려준다. 입력이 깔끔해도 계산 경로가 반올림을 집어넣는다는 뜻이다.

이 작은 수가 중요한 까닭은 이렇다. 최빈값을 pmf.argmax()로 찾았다면 \(p = 1/2\)에서 동점을 못 보고 \(1\) 하나만 집었을 것이다. 이항분포의 보기 1에서도, 포아송분포의 보기 1에서도 같은 일이 일어났다. 동점을 가리는 믿을 만한 길은 분수뿐이고, 베르누이는 그 교훈을 가장 적은 재료로 보여 주는 경우다.

(1) 해석적으로. 본문은 미분해서 \(p = 1/2\)을 얻었지만 미분할 필요가 없다. \(p(1-p)\)를 완전제곱으로 묶으면

이다. 뺄 것이 제곱이라 언제나 \(0\) 이상이고, \(p = 1/2\)에서만 \(0\)이 된다. 따라서 최댓값은 \(1/4\)이고 그 자리는 \(p = 1/2\) 하나다. 이계도함수를 따질 일도 없고 끝점 \(p = 0, 1\)에서 \(0\)이 되는 것도 식에서 바로 읽힌다.

표준편차는 그 제곱근이다. \(u = p - \tfrac12\)로 두면

이다. \(\sqrt{1-t} \approx 1 - t/2\)를 썼다. 1차 항이 없다는 것이 "평평하다"의 뜻이다. \(p\)가 \(1/2\)에서 벗어난 만큼의 제곱으로만 줄어들므로, 조금 벗어나면 거의 손해가 없다. \(p = 0.3\)이면 \(u = -0.2\)라

이고 정확한 값 \(\sqrt{0.21} = 0.4583\)과 \(0.002\)밖에 차이 나지 않는다. 분산 쪽은 근사가 아니라 \(\frac14 - u^2\)가 정확한 등식이라는 점도 함께 기억할 만하다.

(2) 해석적으로. 95% 신뢰수준에서 오차허용폭 \(e\)를 맞추려면

이다. 여기서 문제가 생긴다. 필요한 \(n\)이 아직 모르는 \(p\)에 달려 있다. 조사를 해 봐야 \(p\)를 알고, \(p\)를 알아야 \(n\)을 정할 수 있다.

(1)의 최댓값이 이 고리를 끊는다. \(p(1-p) \le 1/4\)이므로

로 잡으면 \(p\)가 무엇이든 오차허용폭을 지킨다. \(e = 0.03\)이면

이다. 여론조사에서 "표본 1000명, 오차 \(\pm 3\%\)p"라는 문구가 늘 따라붙는 근거가 이 계산이다.

연습문제 4는 같은 계산에서 \(n \ge 1111\)을 얻는다. 어긋난 것이 아니다. 그쪽은 \(1.96/2 = 0.98\)을 \(1\)로 올려 오차한계를 \(1/\sqrt n\)으로 어림했으므로 \(n \ge 1/0.03^2 = 1111\)이 되고, 여기서는 \(1.96\)을 그대로 두어 \(1068\)이 된다. 비는 \((1/0.98)^2 = 1.041\)이라 \(4\%\) 차이다. 어느 쪽이든 보수적인 쪽으로 틀리므로 실무에서 문제가 되지 않는다.

평평함이 보장하는 것은 이 보수적 선택이 크게 낭비는 아니라는 점이다. 참값이 \(p = 0.3\)이었다면 필요한 표본은 \(1.96^2 \times 0.21/0.0009 = 896.4\), 곧 \(897\)명이다. \(1068\)은 그보다 \(0.25/0.21 = 1.19\)배이니 \(19\%\)만 더 뽑은 셈이다. 반면 \(p = 0.1\)처럼 한쪽으로 치우친 비율이라면 \(385\)명으로 족하므로 \(2.78\)배를 과하게 뽑는다. 평평함은 \(p\)가 가운데 근처일 때만 평평하다.

수치적으로. 먼저 쪽의 그림을 그린다.

import matplotlib.pyplot as plt
import numpy as np

p = np.linspace(0, 1, 200)

fig, ax = plt.subplots(figsize=(12, 3))
# 분산 p(1-p)는 위로 볼록한 포물선으로 p=0.5에서 최대 0.25 다.
# 표준편차는 그 제곱근이라 가운데가 더 평평하다.
#   -> p가 0.3이든 0.5든 표준편차가 0.458 대 0.5 로 10%밖에 차이 나지 않는다.
# 여론조사에서 p를 모르고도 표본크기를 정할 수 있는 근거가 이 평평함이다.
ax.plot(p, p * (1 - p), lw=2, label='Var = p(1-p)')
ax.plot(p, np.sqrt(p * (1 - p)), lw=2, label='SD = sqrt(p(1-p))')
ax.axvline(0.5, color='gray', ls=':', lw=1)
ax.set_xlabel('p')
ax.spines[['top', 'right']].set_visible(False)
ax.legend()
plt.show()

아래 곡선(분산)은 포물선이고 위 곡선(표준편차)은 가운데가 눌린 반원꼴이다. 점선이 지나는 \(p = 0.5\)에서 둘 다 꼭대기에 있다. 수로 재 본다.

import math

print(f"{'p':>5}{'Var':>9}{'1/4-(p-1/2)^2':>16}{'SD':>10}{'2차근사':>10}{'오차':>11}{'SD/0.5':>9}")
for p in (0.5, 0.4, 0.3, 0.2, 0.1, 0.02):
    sd = math.sqrt(p * (1 - p))
    approx = 0.5 - (p - 0.5) ** 2        # SD 의 2차 근사
    print(f"{p:>5}{p*(1-p):>9.4f}{0.25-(p-0.5)**2:>16.4f}{sd:>10.6f}"
          f"{approx:>10.6f}{sd-approx:>+11.6f}{sd/0.5:>9.4f}")

# 오차허용폭 0.03, 95% 신뢰수준에서 필요한 표본크기
z, e = 1.96, 0.03
for p in (0.5, 0.3, 0.1):
    n = z**2 * p * (1 - p) / e**2
    print(f"p = {p}:  n = z^2 p(1-p)/e^2 = {n:.1f}  ->  {math.ceil(n)}")
print(f"p=1/2 로 잡으면 p=0.3 보다 {0.25/0.21:.4f}배, p=0.1 보다 {0.25/0.09:.4f}배 많이 뽑는다")

출력:

    p      Var   1/4-(p-1/2)^2        SD      2차근사         오차   SD/0.5
  0.5   0.2500          0.2500  0.500000  0.500000  +0.000000   1.0000
  0.4   0.2400          0.2400  0.489898  0.490000  -0.000102   0.9798
  0.3   0.2100          0.2100  0.458258  0.460000  -0.001742   0.9165
  0.2   0.1600          0.1600  0.400000  0.410000  -0.010000   0.8000
  0.1   0.0900          0.0900  0.300000  0.340000  -0.040000   0.6000
 0.02   0.0196          0.0196  0.140000  0.269600  -0.129600   0.2800
p = 0.5:  n = z^2 p(1-p)/e^2 = 1067.1  ->  1068
p = 0.3:  n = z^2 p(1-p)/e^2 = 896.4  ->  897
p = 0.1:  n = z^2 p(1-p)/e^2 = 384.2  ->  385
p=1/2 로 잡으면 p=0.3 보다 1.1905배, p=0.1 보다 2.7778배 많이 뽑는다

(1)이 다 맞는다.

완전제곱. Var 열과 1/4-(p-1/2)^2 열이 여섯 줄 모두 같다. 근사가 아니라 등식이므로 당연한 결과이지만, 식을 잘못 옮기지 않았다는 확인은 된다.

평평함. SD/0.5 열이 \(p\)를 \(0.5\)에서 \(0.4\)로 옮길 때 \(0.98\), \(0.3\)으로 옮길 때 \(0.92\)다. \(p\)를 \(20\%\) 틀리게 짐작해도 표준편차는 \(2\%\)만 틀린다. 2차 근사의 오차도 \(p = 0.4\)에서 \(-0.0001\), \(p = 0.3\)에서 \(-0.0017\)로 작다. 다만 \(p = 0.1\)에서 \(-0.04\), \(p = 0.02\)에서 \(-0.13\)으로 급히 커지는데, \(u^4\) 항을 버렸기 때문이다. \(u\)가 커지면 근사가 무너진다는 것을 마지막 두 줄이 보여 준다. 본문 광고 문제의 \(p = 0.02\)가 바로 그 영역이다.

표본크기. \(1067.1 \to 1068\), \(896.4 \to 897\), \(384.2 \to 385\)가 (1)의 계산과 같고, 과잉비 \(1.1905\)와 \(2.7778\)도 \(0.25/0.21\), \(0.25/0.09\)와 맞는다. \(p = 1/2\)을 쓰는 보수적 선택은 \(p\)가 \(0.3 \sim 0.7\) 안에 있을 때만 "조금 더"이고, 그 밖에서는 표본을 두 배 넘게 낭비한다.

(1) 해석적으로. \(\hat p_n = \frac1n\sum_{i=1}^n X_i\)이고 \(X_i\)가 독립이므로

이다. \(n\)이 열 배가 되면 표준오차는 \(1/\sqrt{10} = 0.3162\)배가 된다. 본문이 말하는 "\(\sqrt{10} \approx 3.2\)배씩 줄어든다"가 이것이다.

\(n = 10^4\)가 이상한지 따져 보자. 오차의 절대크기가 줄지 않았다면 그것은 표준오차로 재서 \(z\)가 커졌다는 뜻이다. \(n = 10^3\)에서 오차 \(-0.0120\)은 \(-0.83\) SE이지만, \(n = 10^4\)에서 오차 \(-0.0113\)은 \(-2.47\) SE다. \(\lvert z \rvert = 2.47\)은 한 번 보는 데 확률이 \(2\Phi(-2.47) = 0.0137\)이니 드문 편이다. 그러나 다섯 번 보는 동안 적어도 한 번 나올 확률은 \(1 - (1 - 0.0137)^5 = 0.067\)이라 놀랄 일은 아니다. 큰수의 법칙은 오차가 매 단계 줄어든다고 약속하지 않는다. 오차의 크기 척도가 \(1/\sqrt n\)으로 줄어든다고만 약속한다.

(2) 해석적으로. 다섯 값은 서로 독립이 아니다. 뒤의 것이 앞의 것을 포함하고 있기 때문이다. \(m < n\)이면 \(S_n = \sum_{i=1}^n X_i\)에서 \(S_m\)이 \(S_n\)의 일부이므로

이다. \(\text{Cov}(S_m, S_n) = \text{Cov}(S_m, S_m + (S_n - S_m)) = \text{Var}(S_m)\)인데, 뒤 토막 \(S_n - S_m\)이 \(S_m\)과 독립이라 공분산에 기여하지 않기 때문이다. 상관계수로 바꾸면 \(p(1-p)\)가 모두 약분되어

만 남는다. \(p\)에 전혀 의존하지 않는다. 한 자리 떨어진 이웃끼리는 \(\sqrt{1/10} = 0.316\), 두 자리면 \(0.1\), 네 자리면 \(0.01\)이다. 곧 겹쳐 쌓은 평균이라 완전히 독립은 아니지만 상관이 약하다. 다섯 번을 거의 독립인 다섯 번으로 세어도 (1)의 \(0.067\) 계산이 크게 틀리지 않는 근거가 이것이다.

수치적으로. 먼저 쪽의 코드를 그대로 돌린다.

import numpy as np
from scipy import stats

np.random.seed(42)
p = 0.3
samples = stats.bernoulli(p).rvs(100_000)

# 베르누이 표본의 평균은 곧 성공 비율이다. 큰수의 법칙에 따라 p로 간다.
for n in (10, 100, 1_000, 10_000, 100_000):
    print(f"n = {n:>7}:  표본비율 = {samples[:n].mean():.4f}")
print(f"이론값 p = {p},  Var = {p*(1-p):.4f}")

출력:

n =      10:  표본비율 = 0.4000
n =     100:  표본비율 = 0.3000
n =    1000:  표본비율 = 0.2880
n =   10000:  표본비율 = 0.2887
n =  100000:  표본비율 = 0.2989
이론값 p = 0.3,  Var = 0.2100

표본이 열 배가 될 때마다 오차가 평균적으로 \(\sqrt{10} \approx 3.2\)배씩 줄어든다. 표준오차가 \(\sqrt{p(1-p)/n}\)이기 때문이다. 다만 한 번의 실행에서 수렴이 단조롭지는 않다. 위에서도 \(n = 100\)의 오차가 \(n = 1000\)보다 작은데, 우연히 잘 맞은 것일 뿐이다. 큰수의 법칙은 평균적인 경향이지 매 단계의 보장이 아니다.

이제 (1)과 (2)를 수로 확인한다.

import numpy as np
from math import sqrt
from scipy import stats

p = 0.3
np.random.seed(42)
samples = stats.bernoulli(p).rvs(100_000)
ns = (10, 100, 1_000, 10_000, 100_000)

print(f"{'n':>8}{'성공':>7}{'phat':>9}{'오차':>10}{'SE':>10}{'z':>9}")
for n in ns:
    ph = samples[:n].mean()
    se = sqrt(p * (1 - p) / n)
    print(f"{n:>8}{int(samples[:n].sum()):>7}{ph:>9.4f}{ph-p:>+10.4f}{se:>10.5f}{(ph-p)/se:>+9.3f}")

print(f"\nSE 가 한 자리 내려갈 때마다 곱해지는 비: 1/sqrt(10) = {1/sqrt(10):.4f}")

# 다섯 phat 은 겹쳐 쌓은 앞부분 평균이라 서로 독립이 아니다.
# Corr(phat_m, phat_n) = sqrt(m/n)  (m < n) 을 모의실험으로 확인한다.
rng = np.random.default_rng(0)
X = rng.binomial(1, p, size=(40_000, 1000))
cum = X.cumsum(axis=1)
for m, n in ((10, 100), (100, 1000), (10, 1000)):
    ph_m, ph_n = cum[:, m-1] / m, cum[:, n-1] / n
    print(f"Corr(phat_{m}, phat_{n}) 모의 {np.corrcoef(ph_m, ph_n)[0,1]:.4f}   이론 sqrt({m}/{n}) = {sqrt(m/n):.4f}")

# z = -2.47 정도가 다섯 번 보는 동안 한 번쯤 나오는 것이 놀라운가
one = 2 * stats.norm.sf(2.466)
print(f"\n|z| > 2.466 한 번의 확률 {one:.4f},  다섯 번 중 적어도 한 번 {1-(1-one)**5:.4f}")

출력:

       n     성공     phat        오차        SE        z
      10      4   0.4000   +0.1000   0.14491   +0.690
     100     30   0.3000   +0.0000   0.04583   +0.000
    1000    288   0.2880   -0.0120   0.01449   -0.828
   10000   2887   0.2887   -0.0113   0.00458   -2.466
  100000  29889   0.2989   -0.0011   0.00145   -0.766

SE 가 한 자리 내려갈 때마다 곱해지는 비: 1/sqrt(10) = 0.3162
Corr(phat_10, phat_100) 모의 0.3164   이론 sqrt(10/100) = 0.3162
Corr(phat_100, phat_1000) 모의 0.3210   이론 sqrt(100/1000) = 0.3162
Corr(phat_10, phat_1000) 모의 0.1009   이론 sqrt(10/1000) = 0.1000

|z| > 2.466 한 번의 확률 0.0137,  다섯 번 중 적어도 한 번 0.0665

두 유도가 다 맞는다.

표준오차. SE 열이 (1)의 표와 소수점 다섯째 자리까지 같고, 한 줄 내려갈 때마다 \(0.3162\)배가 된다.

\(z\) 열을 보라. \(+0.69\), \(0.00\), \(-0.83\), \(-2.47\), \(-0.77\)이다. 오차의 절대크기만 보면 \(n = 10^3\)과 \(n = 10^4\)가 거의 같아 "수렴이 멈춘 것처럼" 보이지만, 표준오차로 재면 \(-0.83\)에서 \(-2.47\)로 세 배 멀어진 것이다. \(n = 10^2\)에서 표본비율이 정확히 \(0.3000\)으로 떨어진 것도 \(30/100\)이라 격자가 맞았을 뿐이고, \(z = 0\)이 수렴의 증거는 아니다. 오차를 날 것으로 보지 말고 표준오차로 나누어 보라는 것이 이 표의 요점이다.

\(-2.47\)은 놀랄 값인가. 한 번 보는 확률은 \(0.0137\)이지만 다섯 번 중 적어도 한 번은 \(0.0665\)다. 열다섯 번에 한 번쯤 있는 일이고, 씨앗 하나짜리 실행에서 이것을 보았다고 큰수의 법칙을 의심할 이유는 없다. 다만 \(z\)가 한 번이라도 \(4\)를 넘었다면 이야기가 달라진다.

상관. 모의로 잰 \(0.3164\), \(0.3210\), \(0.1009\)가 유도한 \(\sqrt{m/n} = 0.3162\), \(0.3162\), \(0.1000\)과 맞는다. \(40000\)번 반복의 상관계수 표준오차가 \(1/\sqrt{40000} = 0.005\) 수준이니 \(0.3210\)의 어긋남도 그 안이다. \(p = 0.3\)이 식에 전혀 들어가지 않는다는 점도 눈여겨볼 만하다. 겹쳐 쌓은 평균들의 상관은 분포가 무엇이든 \(\sqrt{m/n}\)이고, 이것이 브라운 운동의 자기상관과 같은 모양이다.

(a) \(X\)가 0 또는 1만 가지므로 \(X^k = X\)다. 따라서 \(E[X^k] = E[X] = p\)이다. \(\square\)

(b) 중심적률은

이다. 왜도는 이를 \(\sigma^3 = \{p(1-p)\}^{3/2}\)으로 나눈

이다.

\(p = 1/2\)에서 0(대칭)이고, \(p < 1/2\)이면 양수(오른쪽 꼬리), \(p > 1/2\)이면 음수다. \(p\)가 극단으로 갈수록 왜도의 절댓값이 무한대로 커지는데, 희귀사건일수록 분포가 한쪽으로 심하게 쏠린다는 뜻이다. 이항분포의 왜도 \(\frac{1-2p}{\sqrt{np(1-p)}}\)는 여기에 \(\sqrt n\)만 나눈 꼴이며, \(n\)이 커지면 0으로 가는 것이 정규근사가 작동하는 이유다.

값이 둘뿐이므로 정의대로 더한다.

미분하면 \(M'(t) = pe^t\), \(M''(t) = pe^t\)이므로

이고 \(\text{Var}(X) = p - p^2 = p(1-p)\)이다. \(\square\)

모든 차수의 도함수가 \(pe^t\)로 같다는 점이 \(E[X^k] = p\)(연습문제 1)의 또 다른 표현이다.

이 MGF가 중요한 이유는 곱으로 쌓인다는 데 있다. 독립인 \(n\)개를 더하면 MGF가 \((1-p+pe^t)^n\)이 되는데, 이것이 바로 이항분포의 MGF다. 분포의 덧셈이 MGF의 곱셈이 된다는 성질이 다음 페이지 전체를 떠받친다.

(a) 아니다. 합이 0, 1, 2의 세 값을 가질 수 있으므로 지지집합부터 어긋난다. \(p_1 = p_2 = p\)이면 \(B(2, p)\)이고, 다르면 그것도 아니다(포아송–이항분포).

(b) 그렇다. 곱이 1이 되려면 둘 다 1이어야 하므로

이다. 지시함수로 보면 \(\mathbb{1}_A\mathbb{1}_B = \mathbb{1}_{A \cap B}\)이고 독립이면 \(P(A \cap B) = P(A)P(B)\)라는 이야기와 같다.

(c) 그렇다. 최댓값이 1이 되려면 적어도 하나가 1이면 되므로

이다. 이것도 지시함수로 보면 \(\max(\mathbb{1}_A, \mathbb{1}_B) = \mathbb{1}_{A \cup B}\)이고, 여사건의 곱으로 합집합 확률을 구하는 익숙한 계산이다.

정리하면 베르누이족은 곱과 최댓값·최솟값에 대해 닫혀 있지만 덧셈에 대해서는 닫혀 있지 않다. 0과 1만 갖는 값에 논리 연산(그리고, 또는)을 하면 다시 0과 1이지만, 산술 덧셈은 그 울타리를 벗어나기 때문이다.

\(\text{Var}(X) = p(1 - p)\)를 \(p\)에 대해 미분하면

이다. 0으로 두면 \(p = 1/2\)이다. 이계도함수가 \(-2 < 0\)이므로 최대점이며, 최대 분산은 \(1/4\)이다.

실용적 의미. 이항 비율의 신뢰구간에서 최악의 경우 분산은 \(p(1-p) \le 1/4\)이다. 보수적인 표준오차는 \(\sqrt{1/(4n)} = 1/(2\sqrt n)\)이므로, 95% 오차한계는 최대 \(1.96/(2\sqrt n) \approx 1/\sqrt n\)이다.

오차한계를 \(\le 0.03\)으로 두면 \(n \ge 1/(0.03)^2 \approx 1111\)이 되는데, 이것이 전국 여론조사에서 "n ≈ 1000" 규칙이 나온 배경이다. 실제 \(p\)는 대개 0.5에서 떨어져 있으므로 이 보수적 한계는 다소 느슨하지만, \(p\)가 무엇이든 통하는 표본크기 추정치를 준다는 것이 장점이다.

보기 2의 그림에서 보듯 표준편차 곡선이 가운데에서 평평하다는 점도 함께 기억할 만하다. \(p\)를 0.3으로 잘못 짐작해도 표준편차는 0.458로 최댓값 0.5의 92%라 표본크기 계산이 크게 어긋나지 않는다.

PMF를 한 줄로 써 둔 덕분에 가능도가 곱으로 깔끔하게 정리된다.

로그를 취하고 미분하면

이고, 정리하면 \(s(1-p) = (n-s)p\)에서

를 얻는다. 표본비율이 곧 최대가능도추정량이다.

불편성은 기댓값의 선형성에서 곧바로 나온다.

\(\square\)

분산은 \(\text{Var}(\hat p) = p(1-p)/n\)이며, 이것이 크라메르–라오 하한과 일치하므로 \(\hat p\)는 최소분산 불편추정량이다. 피셔 정보량이 \(I(p) = 1/\{p(1-p)\}\)이므로 하한이 \(1/\{nI(p)\} = p(1-p)/n\)이 되는 것을 직접 확인할 수 있다.

\(N = \sum_{i=1}^n \mathbb{1}_{A_i}\)로 쓰면 기댓값의 선형성에 의해

이다. 독립일 필요가 전혀 없다. 선형성은 결합분포에 아무 조건도 요구하지 않기 때문이다.

분산은 사정이 다르다.

둘째 항의 각 조각은 \(\text{Cov}(\mathbb{1}_{A_i}, \mathbb{1}_{A_j})\)이며, 사건들이 독립일 때만 0이 된다. 평균은 공짜로 얻지만 분산은 의존구조를 알아야 한다.

이 비대칭이 4장 전체에 되풀이해 나타난다. 초기하분포에서 평균이 이항분포와 같고 분산만 유한모집단 수정계수만큼 작았던 것이 바로 이 구조였다.

지시함수 기법의 위력을 보여 주는 고전적인 예가 일치 문제다. \(n\)명이 모자를 무작위로 도로 집어 갈 때 자기 모자를 집는 사람 수의 평균은, 각자가 자기 모자를 집을 확률이 \(1/n\)이므로 \(n \times (1/n) = 1\)이다. \(n\)과 무관하게 언제나 1이며, 사건들이 서로 독립이 아닌데도 한 줄로 끝난다.

미분하면

이고 0으로 두면 \((1-p)/p = 1\), 즉 \(p = 1/2\)다. 이계도함수가 \(-1/\{p(1-p)\} < 0\)이므로 최대이며 최댓값은 \(H(1/2) = \ln 2\)(비트로 재면 1비트)다.

왜 같은 지점인가. 둘 다 "결과가 얼마나 불확실한가"를 재는 양이고, 두 결과가 똑같이 일어날 법할 때 불확실성이 최대이기 때문이다. 대칭성만으로도 극점이 \(p = 1/2\)에 있으리라 짐작할 수 있다.

어떻게 다른가. 끝점에서의 거동이 다르다. \(p \to 0\)일 때

로 엔트로피가 분산보다 느리게 0으로 간다(\(\ln(1/p)\)만큼 더 크다). 희귀사건에서 분산은 거의 사라지지만 엔트로피는 그만큼 빨리 줄지 않는다는 뜻이며, "거의 일어나지 않는 일이 일어났다"는 사건이 담는 정보가 크기 때문이다.

또 분산은 값의 척도에 의존하지만 엔트로피는 그렇지 않다. \(X\) 대신 \(100X\)를 재면 분산은 1만 배가 되지만 엔트로피는 그대로다. 엔트로피는 결과에 붙은 이름표를 바꿔도 변하지 않는 순수한 불확실성의 척도다.

이 구별은 실무에서 나뉘는 지점이 있다. 의사결정나무의 분할 기준으로 쓰는 지니 불순도 \(2p(1-p)\)는 분산 쪽이고, 정보이득은 엔트로피 쪽이다. 둘 다 \(p = 1/2\)에서 최대라 실제 결과가 비슷하게 나오지만, 극단적으로 불균형한 자료에서는 엔트로피 쪽이 소수 범주를 더 크게 취급한다.

정보량. PMF를 \(f(x;p) = p^x(1-p)^{1-x}\)로 쓰면

이고 한 번 더 미분하면

이다. \(E[X] = p\)를 넣어 부호를 뒤집으면

를 얻는다. 분산의 역수라는 점에 주목하라.

크라메르–라오 하한. 독립인 \(n\)개 표본의 정보량은 더해지므로 \(I_n(p) = n/\{p(1-p)\}\)이고, \(p\)의 임의의 불편추정량 \(\tilde p\)에 대해

이다. 그런데 연습문제 5의 최대가능도추정량 \(\hat p = \bar X\)는

로 하한과 정확히 같다. 즉 \(\hat p\)는 유효추정량이며, 불편추정량 중에서 이보다 분산이 작은 것은 존재하지 않는다. \(\square\)

\(I(p)\)의 모양을 읽어라.

정보량이 가장 작은 곳이 \(p = 1/2\)다. 동전이 공정할수록 한 번의 던지기가 \(p\)에 대해 알려 주는 것이 적다. 연습문제 4에서 분산이 \(p = 1/2\)에서 최대였던 것과 같은 사실을 반대편에서 본 것이다.

끝점에서 발산하는 것은 착시가 아니다. \(p\)가 \(0\)에 가까우면 성공 한 번을 관측하는 것만으로 "\(p\)가 그렇게 작지는 않다"는 강한 정보가 들어온다. 다만 이 발산을 "희귀사건 추정이 쉽다"로 읽으면 안 된다. 상대적으로 보면 정반대다. \(\hat p\)의 변동계수는

로 \(p \to 0\)에서 발산한다. 절대 오차는 작아지지만 상대 오차는 커진다. \(p = 0.001\)을 20% 이내로 추정하려면 \(n\)이 수만 단위여야 한다.

이것이 8장에서 비율의 신뢰구간이 극단적인 \(\hat p\)에서 무너지는 이유의 뿌리다. \(\hat p = 0\)이면 왈드 구간의 폭이 \(0\)이 되어 버리는데, 정보량이 크다는 것과 구간이 좁아야 한다는 것은 전혀 다른 이야기다.

\(S = \sum_i X_i \sim B(n,p)\)이고 \(\hat p = S/n\)이다. 이차적률부터 구한다.

따라서

이다. 체계적으로 과소추정한다. 참값의 \(\frac{n-1}{n}\)배이므로

가 불편추정량이다. \(\square\)

\(n-1\)이 또 나왔다. 2장에서 표본분산을 \(n-1\)로 나눈 것과 정확히 같은 보정이다. 우연이 아니다. 베르누이 자료에서 표본분산을 직접 계산해 보면

로 위의 보정과 같은 식이 나온다. 멱등성 \(X_i^2 = X_i\)가 둘을 잇는다. 일반적인 \(n-1\) 보정의 베르누이 특수판인 것이다.

import numpy as np
rng = np.random.default_rng(0)
p, reps = 0.3, 200_000
print(f"참값 p(1-p) = {p*(1-p):.4f}\n")
print(f"{'n':>5}{'E[phat(1-phat)]':>18}{'이론 (n-1)/n':>15}{'보정 후':>12}")
for n in (2, 5, 10, 50):
    S = rng.binomial(n, p, reps)
    ph = S / n
    raw = np.mean(ph * (1 - ph))
    corr = np.mean(n / (n - 1) * ph * (1 - ph))
    print(f"{n:>5}{raw:>18.4f}{(n-1)/n*p*(1-p):>15.4f}{corr:>12.4f}")

출력:

참값 p(1-p) = 0.2100

    n   E[phat(1-phat)]     이론 (n-1)/n        보정 후
    2            0.1049         0.1050      0.2097
    5            0.1679         0.1680      0.2099
   10            0.1891         0.1890      0.2102
   50            0.2058         0.2058      0.2100

\(n = 2\)에서 편향이 절반이나 된다. 그리고 \(n\)이 커지면 \(\frac{n-1}{n} \to 1\)이라 편향이 사라진다. \(n = 50\)이면 이미 2% 차이라 실무에서는 대개 무시한다.

그래도 무시하면 안 되는 자리가 있다. 비율의 표준오차 \(\sqrt{\hat p(1-\hat p)/n}\)는 이 편향된 추정량을 그대로 쓰므로 표준오차를 과소평가한다. 작은 표본에서 신뢰구간이 명목 신뢰수준보다 좁아지는 원인 하나가 여기에 있으며, 8장에서 윌슨 구간이 왈드 구간보다 나은 이유와도 얽혀 있다.

자유도는 하나뿐이다. 주변분포가 고정되어 있으므로 결합분포는 \(\pi_{11} = P(X=1, Y=1)\) 하나로 완전히 결정된다.

네 칸이 모두 \(0\) 이상이어야 하므로

이다. 이것이 프레셰–회프딩 경계다. 한편 \(E[XY] = \pi_{11}\)이므로

이고, \(\pi_{11}\)의 양 끝을 넣으면 \(\rho\)의 범위가 그대로 나온다. \(q_i = 1-p_i\)로 쓰면 \(p_1 \le p_2\)이고 \(p_1+p_2 \le 1\)인 경우

이다. \(\square\)

\(p_1 = p_2 = 0.1\)이면 \(\rho = -1\)은 불가능하다. 위 식에 넣으면

이다. 아무리 강하게 음의 관계를 만들어도 \(-0.111\)보다 내려갈 수 없다. 이유는 표에서 바로 보인다. 완전한 음의 관계라면 \(X=1\)일 때 반드시 \(Y=0\), \(X=0\)일 때 반드시 \(Y=1\)이어야 하는데, 그러려면 \(P(Y=1) = P(X=0) = 0.9\)여야 한다. \(p_2 = 0.1\)과 모순이다.

규칙이 둘 보인다. 주변분포가 같으면(\(p_1 = p_2\)) \(\rho_{\max} = 1\)이 언제나 가능하다. \(X = Y\)로 두면 되기 때문이다. 그리고 \(p_1 + p_2 = 1\)이면 \(\rho_{\min} = -1\)이 가능하다. 그때는 \(Y = 1-X\)로 둘 수 있다.

3.3절과 이어 읽어라. 결합분포 절에서 "주변분포는 결합분포를 결정하지 못한다"고 했는데, 이 연습문제는 그 말의 반쪽을 보여 준다. 주변분포가 결합분포를 정하지는 못하지만 아무것이나 허용하지도 않는다. 주변분포는 결합분포가 놓일 수 있는 구간의 양 끝을 정하고, 그 안에서 의존 구조가 자유로울 뿐이다.

실무적 함의. 희귀사건 두 개의 상관계수는 구조적으로 \(0\) 근처에 갇힌다. 부도, 희귀질환, 클릭처럼 \(p\)가 작은 이진 변수들 사이에서 \(\rho = -0.05\)를 보고 "관계가 거의 없다"고 읽으면 틀릴 수 있다. 가능한 최솟값이 \(-0.111\)인 상황이라면 \(-0.05\)는 가능한 범위의 절반에 해당하는 강한 음의 관계다. 이런 이유로 이진 자료에서는 상관계수 대신 오즈비를 쓰는 일이 많다. 오즈비는 주변분포에 따라 범위가 눌리지 않기 때문이다.