이항 분포

1. 개요

이항 분포 Binomial Distribution
확률질량	$P(X=k) = \binom{n}{k}p^k(1-p)^{n-k}$, $k = 0,\dots,n$
모수	시행 수 $n$, 성공 확률 $p$
평균 · 분산	$np$ · $np(1-p)$
누적분포	정규화 불완전 베타 함수 $I_{1-p}(n-k,\,k+1)$
극한	$np$ 고정 $n\to\infty$ → 포아송, $np(1-p)\to\infty$ → 정규
켤레 사전	베타 분포
시뮬레이션	몬테카를로 고장확률 추정, 이항 도약, 합격·불합격 시험

동전을 $n$ 번 던져 앞면을 센다. 이보다 쉬운 분포는 없는데, 이 분포의 신뢰구간은 통계학 교과서가 수십 년 동안 틀리게 가르친 대표 사례다.

이항 분포(binomial distribution)는 성공 확률이 $p$ 인 독립 시행(베르누이 시행)을 $n$ 번 했을 때 성공 횟수 $X$ 의 분포로,

$$P(X = k) = \binom{n}{k}p^k(1-p)^{n-k},\qquad k = 0, 1, \dots, n$$

이다. 평균은 $np$, 분산은 $np(1-p)$ 이다. 분산이 평균보다 항상 작다는 점(과소분산)이 포아송 분포(분산 = 평균)와의 차이이고, $p \to 0$, $np = \lambda$ 고정 극한에서 포아송이 된다. 그 극한의 정확한 오차(르캉 부등식)는 포아송 분포 문서가 다뤘다.

시뮬레이션 실무에서 이항 분포가 가장 자주 등장하는 자리는 몬테카를로 결과 그 자체다. 표본 $N$ 개 중 고장이 $k$ 개 나왔으면 $k$ 는 $\mathrm{Bin}(N, p_f)$ 이고, 「고장확률이 얼마이고 얼마나 믿을 수 있나」는 이항 비율의 추정 문제다. 이 문서는 그 추정(신뢰구간), 확률질량을 계산할 때의 수치 문제, 그리고 이항 난수를 만드는 법을 다룬다. 누적분포를 불완전 베타 함수로 계산하는 법과 베이즈 A/B 예제는 베타 함수 문서가 이미 다뤘다.

2. 기본 성질과 근사

• 합. 같은 $p$ 의 $\mathrm{Bin}(n_1, p)$ 와 $\mathrm{Bin}(n_2, p)$ 의 독립 합은 $\mathrm{Bin}(n_1 + n_2, p)$ 다. $p$ 가 다르면 이항이 아니라 「포아송 이항 분포」가 된다.
• 최빈값. $\lfloor (n+1)p\rfloor$ ($(n+1)p$ 가 정수면 그것과 하나 작은 값 둘).
• 정규 근사. 중심극한정리에 의해 $X \approx \mathcal N(np, np(1-p))$. 경험칙은 $np$ 와 $n(1-p)$ 가 둘 다 5~10 이상이다. 이산 분포를 연속으로 근사하므로 연속성 보정($k$ 대신 $k + \tfrac12$)이 정확도를 크게 올린다. $n = 40$, $p = 0.3$ 에서 $P(X \le 8)$ 의 정확한 값은 0.11101 인데, 보정 없는 정규 근사는 0.08377(25% 오차), 보정하면 0.11360(2% 오차)이다.
• 포아송 근사. $p$ 가 작고 $n$ 이 크면 $\mathrm{Poisson}(np)$. 전변동 거리가 $p$ 이하로 묶인다.
• 베타 쌍대성. $P(X \ge k) = I_p(k, n-k+1)$. 「$n$ 개 균등 난수 중 $k$ 번째로 작은 값이 $p$ 이하」와 같은 사건이라 성립한다(베타 분포의 순서 통계량 절 참고). 아래 클로퍼-피어슨 구간이 이 관계로 계산된다.

3. 확률질량을 계산할 때

교과서 식을 그대로 옮기면 터진다. $n = 2000$, $k = 1000$, $p = 0.5$ 에서 $\binom{2000}{1000} \approx 2\times10^{600}$ 은 배정밀도로 표현이 안 되고(파이썬 float() 이 OverflowError), $0.5^{2000}$ 은 0 으로 언더플로한다. 정답은 $0.017839$ 로 멀쩡한 크기인데도 그렇다.

• 로그 공간. $\log P = \log\Gamma(n+1) - \log\Gamma(k+1) - \log\Gamma(n-k+1) + k\log p + (n-k)\log(1-p)$ 를 계산한 뒤 지수를 취한다. 위 예에서 $0.0178390111458544$ 로 정확히 나온다. 다만 큰 수 셋의 차이라 $n$ 이 $10^{8}$ 쯤 되면 로그 감마의 절대오차가 상대오차로 번진다.
• 편차 형태. 그래서 R 의 dbinom 등은 로더(2000)의 방법을 쓴다. 스털링 근사의 오차항과 「편차」 함수 $\mathrm{bd}_0(x, np) = x\log(x/np) + np - x$ 를 따로 정확하게 계산해 큰 수끼리의 소거를 피한다. 그래서 $n$ 이 매우 커도 상대오차가 기계 정밀도 근처에 머문다.¹
• 누적은 베타로. $\sum_{j\le k}$ 를 직접 더하지 말고 불완전 베타 함수 한 번으로(베타 함수 참고).

4. 신뢰구간 — 왈드 구간은 쓰지 마라

$k$ 번 성공을 관측했을 때 $p$ 의 95% 구간. 교과서의 첫 답은 왈드 구간

$$\hat p \pm 1.96\sqrt{\hat p(1-\hat p)/n},\qquad \hat p = k/n$$

이다. 브라운·카이·다스굽타(2001)는 이 구간의 실제 포함확률(coverage)이 「$n$ 이 충분히 크면 괜찮다」는 통념과 달리 $p$ 가 0.5 여도, $n$ 이 커도 들쭉날쭉하다는 것을 체계적으로 보였다. 정확한 포함확률은 모든 $k$ 에 대해 이항 확률을 더해 바로 계산할 수 있다.

$n$	참 $p$	왈드	윌슨	클로퍼-피어슨
50	0.05	0.920	0.962	0.988
100	0.01	0.633	0.921	0.982
1000	0.001	0.632	0.920	0.981
100	0.5	0.943	0.943	0.965

$np = 1$ 이면 왈드 구간은 $n$ 이 100 이든 1000 이든 포함확률이 63% 다. 이유는 간단하다. 확률 $(1-p)^n \approx e^{-1} = 0.37$ 로 $k = 0$ 이 나오고, 그러면 $\hat p = 0$ 이라 구간이 $[0, 0]$ 으로 쪼그라들어 참값을 놓친다. $p = 0.5$ 에서도 $n = 10, 11, \dots, 20$ 의 포함확률이 0.891, 0.935, 0.854, 0.908, 0.943, 0.882, 0.923, 0.951, 0.904, 0.936, 0.959 로 지그재그한다 — 이산성 때문에 「운 좋은 $n$」과 「운 나쁜 $n$」이 번갈아 온다.

대안은 둘이다.

• 윌슨 점수 구간(1927). 검정 통계량 $|\hat p - p|/\sqrt{p(1-p)/n} \le z$ 를 $p$ 에 대해 풀어서 얻는다(분모에 $\hat p$ 대신 $p$). 중심이 $\hat p$ 에서 $\tfrac12$ 쪽으로 당겨지고, $k = 0$ 이어도 폭이 0 이 아니다. 포함확률이 명목값 근처에서 진동하며 평균적으로 잘 맞는다. 브라운 등이 권한 기본값이다.
• 클로퍼-피어슨 「정확」 구간(1934). 베타 분위수 $[B_{0.025}(k, n-k+1),\ B_{0.975}(k+1, n-k)]$. 포함확률이 모든 $p$ 에서 95% 이상임이 보장된다. 대가는 보수성 — 표에서 보듯 실제로 98% 를 넘기며 구간이 넓다.

규제·안전 문서처럼 「95% 이상 보장」이 요구되면 클로퍼-피어슨, 일반적인 추정에는 윌슨(또는 아그레스티-쿨 「두 번 성공·두 번 실패 더하기」 근사)이 무난하다².

셋의 법칙. $n$ 번 시험에서 성공(고장)이 0 번이면 클로퍼-피어슨 단측 95% 상한은 $1 - 0.05^{1/n} \approx 3/n$ 이다. $n = 100$ 이면 0.0295, $n = 1000$ 이면 0.00299.³ 「1000 번 돌려 한 번도 안 터졌으니 고장확률 0」이 아니라 「0.3% 이하라고 말할 수 있다」가 정확한 진술이다.

5. 몬테카를로 고장확률의 상대오차

고장확률 $p_f$ 를 표본 $N$ 개로 추정하면 $\hat p_f = k/N$ 의 변동계수는

$$\mathrm{CoV}(\hat p_f) = \sqrt{\frac{1 - p_f}{N p_f}}$$

이다. $p_f = 10^{-3}$ 이면 $N = 10^4$ 에서 31.6%, $N = 10^6$ 에서 3.2% 다. 상대오차 10% 를 원하면 대략 $N \approx 100/p_f$ 개가 필요하다 — 기대 고장 횟수가 100 번은 나와야 한다. $p_f = 10^{-6}$ 급 구조 신뢰성 문제에서 순수 몬테카를로가 비현실적인 이유이고, 중요도 표본추출·부분집합 시뮬레이션이 존재하는 이유다. 이들 방법의 성능 지표는 결국 「같은 CoV 를 얻는 데 필요한 모형 호출 수가 이항 기준선보다 몇 배 적은가」로 보고된다.

6. 이항 난수 만들기

• 직접법. 균등 난수 $n$ 개를 뽑아 $p$ 미만을 센다. $O(n)$ 이라 $n$ 이 크면 느리다.
• 역변환. $P(X = 0) = (1-p)^n$ 에서 시작해 $P(X=k+1)/P(X=k) = \frac{n-k}{k+1}\frac{p}{1-p}$ 로 누적분포를 올라가며 균등 난수를 넘는 지점을 찾는다(역변환 표본추출). 기대 비용 $O(np)$ 라 $np$ 가 작을 때 좋다. 단 $n$ 이 크고 $p$ 가 크면 $(1-p)^n$ 이 언더플로해 0 에서 출발하는 버그가 생긴다 — $p > 0.5$ 면 $1-p$ 로 뒤집어 뽑고 $n - X$ 를 반환하는 것이 표준 처리다.
• BTPE. 카치트비차야누쿨·슈마이저(1988)의 「Binomial, Triangle, Parallelogram, Exponential」 방법. 분포 중심을 삼각형·평행사변형·지수 꼬리로 덮는 기각표본추출로, $np$ 와 무관하게 기대 비용이 $O(1)$ 이다. NumPy 등 다수 라이브러리가 $np$ 가 클 때(대략 30 이상) 이것을, 작을 때 역변환을 쓴다.

난수 생성기의 품질 검사에서도 이항 분포가 기본 도구다. 비트 하나하나가 공정한 동전이라면 $n$ 비트 중 1 의 개수는 $\mathrm{Bin}(n, 1/2)$ 이어야 하고, NIST 난수 검사 묶음의 첫 시험(frequency test)이 정확히 이것이다.

7. 화학 시뮬레이션의 이항 도약

길레스피 알고리즘의 타우 도약은 각 반응의 발생 횟수를 포아송으로 뽑는데, 포아송은 상한이 없어 분자 수가 음수가 될 수 있다(포아송 분포 참고). 이항 도약(binomial τ-leaping)은 반응물 분자 수 $x$ 를 시행 수로, $a\tau/x$ 를 성공 확률로 하는 $\mathrm{Bin}(x, a\tau/x)$ 를 대신 뽑는다. 평균은 같고 상한이 $x$ 로 묶여 음수가 원천적으로 불가능하다. 반응물이 둘 이상인 반응에서는 상한을 어느 종으로 잡느냐의 처리가 추가로 필요하다.

8. 관련 문서

• 포아송 분포 · 베타 분포 · 베타 함수 · 음이항 분포 · 디리클레 분포 · 감마 분포
• 중심극한정리 · 최대우도추정 · 켤레사전분포 · 신뢰구간
• 몬테카를로 방법 · 중요도 표본추출 · 부분집합 시뮬레이션 · 신뢰성 해석
• 역변환 표본추출 · 기각표본추출 · 난수 생성기 · 길레스피 알고리즘 · 스털링 수

9. Footnotes

1. 로더의 이 알고리즘은 논문이 아니라 벨 연구소 기술 메모로 돌아다니다가 R 에 들어가 사실상 표준이 됐다. 「이항 확률 하나 계산하는 데 메모 한 편」이 과해 보이지만, $n=10^9$ 에서 로그 감마 셋을 빼는 순간 유효숫자 아홉 자리가 증발하는 걸 보면 납득이 간다. ↩

2. 왈드 구간이 교과서에 오래 살아남은 이유는 한 줄로 쓸 수 있어서다. 브라운 등의 논문에 달린 토론문들을 보면 「다들 알고 있었지만 아무도 안 고쳤다」는 분위기가 역력하다. 국룰의 관성은 포함확률 63% 보다 강하다. ↩

3. 「셋의 법칙」은 의학 임상시험에서 「부작용 0 건」 보고를 해석할 때 단골로 나온다. 시뮬레이션 쪽 버전은 「몬테카를로 1000 번 돌렸는데 발산 0 번」 — 그래도 발산 확률은 0.3% 까지 남아 있다. 그 0.3% 가 학회 데모 당일에 터진다. ↩