음이항 분포

1. 개요

음이항 분포 Negative binomial distribution
정의 1	성공확률 $p$ 시행에서 $r$ 번째 성공 전까지의 실패 횟수
정의 2	비율이 감마 분포인 포아송의 혼합
모수 (생태·회귀)	평균 $\mu$, 산포 $k>0$ (작을수록 뭉침)
분산	$\mu+\mu^2/k$
극한	$k\to\infty$ 이면 포아송, $r=1$ 이면 기하 분포
쓰이는 곳	과대산포 계수 자료, RNA-seq, 생태 조사, 사고 횟수, 버스트 유전자 발현

포아송으로 맞췄더니 0 이 4800 배 많이 나왔다 — 그럴 때 꺼내는 분포.

음이항 분포(negative binomial distribution)는 성공확률이 $p$ 인 독립 시행을 반복할 때 $r$ 번째 성공이 나오기 전까지의 실패 횟수 $K$ 의 분포다.

$$P(K = k) = \binom{k + r - 1}{k}\,p^r(1-p)^k, \qquad k = 0, 1, 2, \dots$$

평균은 $r(1-p)/p$, 분산은 $r(1-p)/p^2$ 로 분산이 항상 평균보다 크다. 이름의 「음」은 이 확률들이 음의 지수 이항급수 $(1 - q)^{-r} = \sum_k \binom{k+r-1}{k}q^k$ 의 항이라서 붙었다. 이항 분포가 「$n$ 번 던져 성공 몇 번」이라면 음이항은 「성공 $r$ 번 볼 때까지 실패 몇 번」으로, 세는 대상과 고정하는 대상이 뒤바뀐 쌍이다. $r = 1$ 이면 기하 분포, $r$ 이 정수일 때는 파스칼 분포라고도 부른다.

그런데 시뮬레이션과 통계 실무에서 음이항을 쓰는 이유는 이 동전 던지기 해석이 아니다. $r$ 을 양의 실수 $k$ 로 풀어 주면 음이항은 「비율이 들쭉날쭉한 포아송 분포」, 즉 감마 분포로 섞은 포아송이 되고, 이것이 과대산포(overdispersion) 계수 자료의 표준 모형이다. 이 문서는 그쪽 — 평균-산포 모수화, 0 과 꼬리, 산포 모수 추정의 함정, 계산 — 을 다룬다. 감마-포아송 혼합의 유도 자체는 감마 분포 문서가 이미 했으니 여기선 결과만 쓴다.

2. 두 모수화 — 실수가 가장 많이 나는 곳

같은 분포를 분야마다 다르게 쓴다. 생태학·회귀·생물정보학은 평균 $\mu$ 와 산포 모수 $k$(또는 그 역수 $\alpha = 1/k$)를 쓴다:

$$P(K = k') = \frac{\Gamma(k' + k)}{\Gamma(k)\,k'!}\Big(\frac{k}{k+\mu}\Big)^{k}\Big(\frac{\mu}{k+\mu}\Big)^{k'}, \qquad \operatorname{Var} = \mu + \frac{\mu^2}{k}.$$

대응은 $r = k$, $p = k/(k+\mu)$ 다. 분산식의 첫 항 $\mu$ 는 포아송 잡음, 둘째 항 $\mu^2/k$ 는 비율 자체의 흔들림(감마의 분산)이다. $k \to \infty$ 면 포아송으로 돌아가고, $k$ 가 작을수록 「뭉침(clumping)」이 심하다. 생태학의 블리스·피셔(1953)가 정착시킨 표기이고, 그린우드·율(1920)이 공장 노동자의 사고 횟수를 「사람마다 다른 사고 경향(accident proneness)」으로 설명하며 처음 이 혼합을 썼다.¹

라이브러리마다 기본 모수가 다르다. scipy.stats.nbinom(n, p) 는 「$n$ 번 성공 전 실패 수, 성공확률 $p$」 정의라 평균 $\mu$, 산포 $k$ 로 쓰려면 nbinom(k, k/(k+mu)) 로 넣어야 한다. R 의 dnbinom 은 size·prob 와 size·mu 를 둘 다 받고, 회귀 패키지는 $\alpha = 1/k$ 를 「dispersion」이라 부르기도 한다. 같은 이름 「산포」가 어떤 데선 $k$, 어떤 데선 $1/k$ 다 — 숫자를 옮길 때 분산식으로 한 번 검산하는 게 안전하다.²

3. 0 과 꼬리 — 같은 평균, 전혀 다른 분포

평균 $\mu = 10$ 으로 고정하고 $k$ 만 바꿔 보면(scipy 로 계산):

모형	분산	$P(K=0)$	$P(K \ge 30)$
포아송	10	$4.5\times10^{-5}$	$2.5\times10^{-7}$
음이항 $k=100$	11	$7.3\times10^{-5}$	$1.3\times10^{-6}$
음이항 $k=10$	20	$9.8\times10^{-4}$	$5.3\times10^{-4}$
음이항 $k=2$	60	0.028	0.025
음이항 $k=0.5$	210	0.218	0.088

$k = 0.5$ 에서 0 의 확률은 포아송의 약 4800 배, 30 이상의 확률은 약 35 만 배다. 0 은 $(k/(k+\mu))^k$ 라 $k$ 가 작으면 천천히만 줄어들고, 꼬리는 기하급수로 꺾이는 대신 훨씬 느리게 내려온다. 그래서 음이항은 0 과잉과 두꺼운 꼬리를 모수 하나로 동시에 설명한다. 「0 이 너무 많으니 영과잉 모형(zero-inflated)부터」로 직행하기 전에 음이항이 0 을 이미 설명하는지 먼저 봐야 하는 이유다 — 영과잉은 「구조적으로 셀 수 없는 0」이라는 별도의 서사가 있을 때 쓰는 것이다.³

포아송 대비 판정은 포아송 분포 문서의 산포 지수(분산/평균)로 하고, 음이항에서는 그 값이 $1 + \mu/k$ 다. 평균이 커질수록 산포 지수도 선형으로 커진다는 것이 음이항의 지문이다. 산포 지수가 평균과 무관하게 일정하면(분산 $= \phi\mu$) 오히려 준포아송이나 NB1 이 맞다.

4. 어디서 나오는가

• 이질성. 개체마다 비율 $\lambda$ 가 다르고 그 분포가 감마로 근사되면, 개체를 섞어 센 횟수는 음이항이다. 사고 경향, 숙주마다 다른 기생충 수, 지점마다 다른 포획 수가 이 꼴이다.
• 버스트 유전자 발현. 단백질(또는 mRNA)이 속도 $a$ 로 한꺼번에 여러 개씩 — 평균 $b$ 인 기하 분포 크기의 버스트로 — 만들어지고 개당 속도 $\gamma$ 로 분해되면, 화학 마스터 방정식의 정상해가 정확히 $r = a/\gamma$, $p = 1/(1+b)$ 인 음이항이다(샤레자이·스웨인 2008). $a = 2$, $b = 5$, $\gamma = 1$ 에서 상태 400개로 잘라 정상분포를 고유벡터로 구하면 평균 10.000, 음이항 PMF 와의 최대 차이 $1.7\times10^{-16}$ 이다. 즉 단일세포 계수의 과대산포는 「측정 잡음」이 아니라 전사 동역학의 지문일 수 있다.
• RNA-seq. 생물학적 반복 사이의 발현 차이 때문에 리드 수가 포아송보다 넓게 퍼진다. DESeq2·edgeR 은 유전자마다 음이항 일반화 선형 모형을 맞추고, 반복이 3개뿐인 실험에서 유전자별 산포를 그대로 추정하면 너무 흔들리므로 유전자들 사이에서 산포 추정을 공유·수축한다(경험적 베이즈).
• 기다림. 원래 정의대로 「$r$ 번째 성공까지 실패 수」. 테스트를 통과할 때까지의 재시도 횟수, 패킷 재전송 수 같은 것.

5. 산포 모수 추정 — 적률법은 깨질 수 있다

적률법은 분산식을 뒤집어 $\hat k = \bar x^2/(s^2 - \bar x)$ 로 쓴다. 한 줄이라 편하지만 문제가 둘 있다.

1. 표본분산이 평균보다 작으면 $\hat k$ 가 음수(또는 무한대)가 된다. 과대산포가 약하거나 표본이 작으면 실제로 일어난다.
2. $k$ 가 작을 때 효율이 떨어진다. 앤스콤(1950)이 따진 것이 이것이다.

진짜 $\mu = 2$, $k = 0.5$ 인 감마-포아송 표본($n = 50$)을 2000번 뽑아 비교하면:

추정량	10% 분위	중앙값	90% 분위	$\log(\hat k/k)$ 의 RMS
적률법	0.354	0.575	0.916	0.402
최대우도	0.353	0.528	0.792	0.329

최대우도 쪽이 중앙값이 참값에 가깝고 위쪽 꼬리가 짧다. $n = 30$ 에서는 2000번 중 한 번 적률법이 아예 실패했고($s^2 \le \bar x$), 같은 표본에서 최대우도는 $\hat k > 10^3$ 으로 달아났다 — 「과대산포가 안 보인다」를 각자의 방식으로 말한 것이다. 이때는 경계 $k = \infty$(포아송)를 정직하게 보고하는 것이 맞다.

최대우도추정의 점수방정식은 $\mu$ 에 대해서는 $\hat\mu = \bar x$ 로 닫히지만 $k$ 에 대해서는

$$\sum_i\big[\psi(x_i + k) - \psi(k)\big] + n\log\frac{k}{k + \bar x} = 0$$

이라 디감마 함수 $\psi$ 가 든 1차원 근 찾기가 된다. $\log k$ 를 변수로 두고 브렌트 방법이나 유계 1차원 최소화를 쓰면 안정적이다. $\psi(x_i + k) - \psi(k) = \sum_{j=0}^{x_i - 1} 1/(k + j)$ 라 정수 자료에서는 디감마 없이 합으로도 계산할 수 있지만, $x_i$ 가 수천인 RNA-seq 에서는 디감마가 훨씬 싸다.

6. 계산 — 감마 함수는 로그로

• PMF 는 로그 공간에서. $\Gamma(k' + k)/(\Gamma(k)\,k'!)$ 를 그대로 계산하면 $k'$ 이 170 만 넘어도 배정밀도 계승이 넘친다. gammaln 으로 $\log P$ 를 계산하고 마지막에 지수화한다. 우도 계산은 애초에 로그 PMF 만 필요하다.
• $(k/(k+\mu))^k$ 의 소실. $k$ 가 크고 $\mu/k$ 가 작으면 $k\log(1 - \mu/(k+\mu))$ 를 log1p 로 써야 포아송 극한 쪽에서 자릿수를 안 잃는다.
• 표본 생성. 감마-포아송 혼합을 그대로 쓴다: $\lambda \sim \mathrm{Gamma}(k, \mu/k)$ 를 뽑고 $\mathrm{Poisson}(\lambda)$ 를 뽑는다. 실수 $k$ 에도 그대로 되고 numpy 의 negative_binomial 도 이 방식이다. 정수 $r$ 에서 「실패를 세며 동전을 던지는」 방식은 $p$ 가 작으면 기대 반복 수가 $r/p$ 라 느리다.

7. 관련 문서

• 포아송 분포 · 감마 분포 · 이항 분포 · 베타-이항 분포 · 기하 분포
• 일반화 선형 모형 · 최대우도추정 · 디감마 함수 · 계층 베이즈 모형 · 영과잉 모형
• 화학 마스터 방정식 · 단일세포 RNA 시퀀싱 · 지수족

8. Footnotes

1. 같은 자료를 「한 번 사고 난 사람은 다음 사고가 쉬워진다(전염)」로 설명해도 음이항이 나온다(폴리아 과정). 분포만 보고는 이질성과 전염을 구별할 수 없다는 것이 이 바닥의 오래된 결론이다 — 사고 난 직원이 억울해할 만한 대목. ↩

2. 논문 A 의 「dispersion 0.1」과 논문 B 의 「dispersion 10」이 같은 데이터일 수 있다. 서로 상대가 틀렸다고 리뷰 코멘트를 다는 건 국룰이고, 분산식 한 줄이면 화해가 된다. ↩

3. 리뷰어가 「0 이 많으니 ZIP 쓰세요」라고 하면 일단 음이항부터 돌려 보자. 0 이 다 설명되면 영과잉 성분의 추정치가 0 근처로 붙고, 모수 하나 더 쓴 대가로 AIC 만 나빠진다. ↩