Kwangmin Kim - 음이항 분포 (Negative Binomial Distribution)

1 음이항 분포란 무엇인가

이항 분포가 “고정된 \(n\) 번 시행에서 성공 횟수”를 모형화한다면, 음이항 분포(Negative Binomial Distribution)는 관점을 뒤집어 “\(r\) 번째 성공까지 필요한 실패 횟수”를 모형화한다 (Casella & Berger, 2002, Ch.3).

이 관점 전환이 중요한 이유:

가변 표본 크기: 원하는 수의 성공을 얻을 때까지 시행을 계속하는 설계(역이항 표본 추출)를 자연스럽게 다룬다
과산포 모형: \(\text{Var}(Y) > E[Y]\) 를 자연스럽게 표현하여, 포아송의 등산포 가정이 위반될 때 대안이 된다
포아송의 일반화: \(r \to \infty\) 극한에서 포아송에 수렴하므로, 포아송을 특수한 경우로 포함한다

2 정의

2.1 형태 1: \(r\) 번째 성공이 발생하는 시행 번호

정의: 음이항 분포 (시행 수 형태)

독립 베르누이 \((p)\) 시행에서, \(r\) 번째 성공이 발생하는 시행 번호 \(X\) 의 PMF:

\[ P(X = x \mid r, p) = \binom{x-1}{r-1} p^r (1-p)^{x-r}, \quad x = r, r+1, r+2, \ldots \]

유도: \(X = x\) 라는 사건은 다음 두 조건의 교집합이다:

처음 \(x-1\) 번의 시행에서 정확히 \(r-1\) 번 성공: 확률 \(\binom{x-1}{r-1}p^{r-1}(1-p)^{x-r}\)
\(x\) 번째 시행에서 성공: 확률 \(p\)

독립성에 의해 이 두 확률을 곱하면 PMF가 도출된다.

2.2 형태 2: \(r\) 번째 성공 전 실패 횟수 (표준 형태)

\(Y = X - r\) (실패 횟수)로 정의하면:

정의: 음이항 분포 (실패 수 형태, 표준)

\[ P(Y = y \mid r, p) = \binom{r+y-1}{y} p^r (1-p)^y, \quad y = 0, 1, 2, \ldots \]

직관: 음이항 PMF \(\binom{r+y-1}{y} p^r (1-p)^y\) 의 구조

\(p^r\): \(r\) 번 성공할 확률이다.
\((1-p)^y\): \(y\) 번 실패할 확률이다.
\(\binom{r+y-1}{y}\): \(r\) 번째 성공 전에 \(y\) 번 실패가 배치될 수 있는 방법의 수이다.

핵심 아이디어: 마지막 (\(r+y\) 번째) 시행은 반드시 \(r\) 번째 성공으로 고정된다. 따라서 나머지 \(r+y-1\) 번 시행 중 \(y\) 번을 실패 위치로 선택하면 되므로, 배치 방법의 수가 \(\binom{r+y-1}{y}\) 이다.

이항 PMF와의 비교:

이항: \(\binom{n}{y}p^y(1-p)^{n-y}\) — 시행 수 \(n\) 고정, 성공 수 \(y\) 가변
음이항: \(\binom{r+y-1}{y}p^r(1-p)^y\) — 성공 수 \(r\) 고정, 실패 수 \(y\) 가변

관점을 뒤집으면 수식의 모양도 뒤집힌다.

이 형태가 Casella & Berger의 표준 형태이며, 이하에서는 이 정의를 기본으로 사용한다.

3 명칭의 유래: 음의 이항계수

“음이항(negative binomial)”이라는 이름은 이항계수의 확장에서 비롯된다.

\[ \binom{r+y-1}{y} = (-1)^y \binom{-r}{y} \]

여기서 음의 이항계수는 다음과 같이 정의된다:

\[ \binom{-r}{y} = \frac{(-r)(-r-1)(-r-2)\cdots(-r-y+1)}{y!} \]

이를 PMF에 대입하면:

\[ P(Y = y) = (-1)^y \binom{-r}{y} p^r (1-p)^y \]

이 형태가 이항 분포 \(\binom{n}{y}p^y(1-p)^{n-y}\) 와 구조적으로 유사하되, 지수 \(n\) 이 음의 정수 \(-r\) 로 대체된 것이다. 이것이 “음이항”이라는 이름의 수학적 근거이다.

3.1 정규화

\(\sum_{y=0}^{\infty} P(Y = y) = 1\) 의 확인은 음의 이항급수 확장(generalized binomial series)에 기반한다:

\[ (1-q)^{-r} = \sum_{y=0}^{\infty} \binom{r+y-1}{y} q^y, \quad |q| < 1 \]

\(q = 1-p\) 를 대입하면:

\[ \sum_{y=0}^{\infty} \binom{r+y-1}{y} p^r (1-p)^y = p^r \cdot p^{-r} = 1 \quad \checkmark \]

4 평균과 분산

4.1 기대값 유도 (축소 기법)

\[ E[Y] = \sum_{y=1}^{\infty} y \binom{r+y-1}{y} p^r (1-p)^y \]

항등식 \(y\binom{r+y-1}{y} = r\binom{r+y-1}{y-1}\) 을 적용하고 \(z = y-1\) 로 치환하면:

\[ E[Y] = r(1-p) \cdot p^{r} \cdot p^{-(r+1)} \sum_{z=0}^{\infty} \binom{(r+1)+z-1}{z} p^{r+1}(1-p)^z \]

합은 음이항 \((r+1, p)\) 의 전체 확률 합 = 1이므로:

\[ \boxed{E[Y] = \frac{r(1-p)}{p}} \]

직관: 기대값 \(E[Y] = \frac{r(1-p)}{p}\)

\(r\) 번 성공하려면 평균적으로 \(r/p\) 번 시행해야 한다 (첫 성공까지 평균 \(1/p\) 번, \(r\) 번 반복).

그 \(r/p\) 번 시행 중 성공 \(r\) 번을 제외한 나머지가 기대 실패 횟수이다:

\[\frac{r}{p} - r = r \cdot \frac{1-p}{p}\]

구체적 예시: \(p = 0.4\), \(r = 5\) 이면 \(E[Y] = 5 \times 0.6/0.4 = 7.5\). 성공 확률 40%인 시행에서 5번 성공하려면 평균 7.5번의 실패가 필요하다.

\(p\) 가 작을수록(성공이 드물수록) 기대 실패 횟수가 기하급수적으로 커진다. \(p = 0.1\), \(r = 5\) 이면 \(E[Y] = 45\) — 같은 5번 성공을 위해 평균 45번이나 실패해야 한다.

4.2 분산

유사한 축소 기법으로:

\[ \boxed{\text{Var}(Y) = \frac{r(1-p)}{p^2}} \]

직관: 분산 \(\text{Var}(Y) = \frac{r(1-p)}{p^2}\) 와 과산포

\(\text{Var}(Y) = E[Y]/p\) 이므로, \(p < 1\) 이면 분산이 항상 평균보다 크다. 이것이 음이항 분포가 포아송(평균 = 분산)보다 더 퍼진 분포인 이유이다.

왜 \(p^2\) 가 분모에 있는가: 실패 횟수 \(Y\) 의 불확실성은 두 층에서 기인한다.

각 시행의 결과가 불확실하다 (이항에서 \(p(1-p)\) 기여).
총 시행 횟수 자체도 불확실하다 (고정된 \(n\) 이 없으므로).

이 이중 불확실성이 분모에 \(p^2\) 를 만든다. 이항(고정 \(n\))의 분산 \(np(1-p)\) 보다 항상 크다.

재모수화: \(\text{Var}(Y) = \mu + \mu^2/r\) 에서 첫 항 \(\mu\) 는 포아송과 같은 기여이고, \(\mu^2/r\) 은 개체 간 이질성(heterogeneity)에서 오는 추가 변동이다. \(r\) 이 클수록 이질성이 작아져 포아송에 가까워진다.

검증: \(\text{Var}(Y) = E[Y]/p = E[Y] \cdot (1 + E[Y]/r)\) 로, 분산이 평균보다 항상 크다 ( \(p < 1\) 일 때).

5 과산포 재모수화

음이항 분포의 핵심적 실무 가치는 과산포(overdispersion)를 자연스럽게 표현하는 것이다.

\(\mu = r(1-p)/p\) (평균)로 재모수화하면:

\[ \boxed{\text{Var}(Y) = \mu + \frac{\mu^2}{r}} \]

이 표현은 분산이 평균의 이차 함수임을 보여준다.

모수	분산-평균 관계	해석
\(r\) 큼	\(\text{Var} \approx \mu\)	포아송에 가까움 (약한 과산포)
\(r\) 작음	\(\text{Var} \gg \mu\)	강한 과산포
\(r \to \infty\)	\(\text{Var} = \mu\)	포아송 극한 (등산포)

실무에서의 과산포 원인

카운트 데이터에서 과산포가 발생하는 흔한 원인:

관측되지 않은 이질성(unobserved heterogeneity): 개체마다 실제 발생률 \(\lambda\) 가 다르지만, 모형은 하나의 \(\lambda\) 를 가정
군집(clustering): 사건이 독립적이지 않고 군집적으로 발생
영과잉(zero-inflation): 0이 포아송 예측보다 훨씬 많이 관측됨

음이항 분포는 첫 번째 원인을 자연스럽게 포착한다. \(\lambda\) 가 감마 분포를 따르는 포아송 혼합(Poisson-Gamma mixture)이 바로 음이항 분포이다.

5.1 포아송-감마 혼합 해석

\(X \mid \lambda \sim \text{Poisson}(\lambda)\) 이고 \(\lambda \sim \text{Gamma}(r, \beta)\) 이면, \(\lambda\) 를 적분하여 주변 분포를 구하면:

\[ P(X = x) = \int_0^{\infty} \frac{e^{-\lambda}\lambda^x}{x!} \cdot \frac{\lambda^{r-1}e^{-\lambda/\beta}}{\Gamma(r)\beta^r} d\lambda = \binom{r+x-1}{x}\left(\frac{1}{1+\beta}\right)^r\left(\frac{\beta}{1+\beta}\right)^x \]

이는 \(p = 1/(1+\beta)\) 인 음이항 분포이다. 즉, 음이항 = 포아송의 강도 \(\lambda\) 가 감마 분포를 따르는 혼합 모형이다. 이 해석이 음이항 회귀의 이론적 기반이다.

6 포아송 극한

정리: 음이항 → 포아송 수렴

\(r \to \infty\) , \(p \to 1\) , \(r(1-p) = \lambda\) (고정)이면:

\[ \text{NegBin}(r, p) \xrightarrow{d} \text{Poisson}(\lambda) \]

적률 수렴 확인:

\[ E[Y] = \frac{r(1-p)}{p} = \frac{\lambda}{p} \to \lambda \]

\[ \text{Var}(Y) = \frac{r(1-p)}{p^2} = \frac{\lambda}{p^2} \to \lambda \]

평균과 분산이 모두 \(\lambda\) 로 수렴하여 등산포가 된다. 이는 과산포 항 \(\mu^2/r\) 이 \(r \to \infty\) 에서 0으로 소멸하기 때문이다.

직관: 음이항 → 포아송 극한 (\(r \to \infty\), \(p \to 1\))

\(r \to \infty\), \(p \to 1\), \(r(1-p) = \lambda\) 고정이라는 조건을 시나리오로 이해하면:

“성공 기준을 아주 높게 설정하되(\(r\) 매우 큼), 각 성공이 매우 쉽게(\(p \to 1\)) 달성되고, 총 기대 실패 횟수 \(\lambda\) 는 고정.”

이 극한에서 실패 사건은 드물고 독립적으로 산발적으로 발생한다 — 이것이 바로 포아송 과정의 희소성 조건이다.

수식으로: 과산포 기여 \(\mu^2/r = \lambda^2/r \to 0\) 으로 사라지므로 분산이 평균과 같아진다. 음이항의 추가 변동성이 완전히 소멸하여 포아송의 등산포 상태가 된다.

7 기하 분포와의 관계

\(r = 1\) 인 음이항이 기하 분포이다:

\[ P(Y = y \mid r=1, p) = p(1-p)^y, \quad y = 0, 1, 2, \ldots \]

기하 분포는 “첫 번째 성공까지의 실패 수”이고, 음이항은 “\(r\) 번째 성공까지의 실패 수”이다. 음이항은 기하의 자연스러운 확장이다.

합의 관계: \(Y_1, \ldots, Y_r \overset{\text{iid}}{\sim} \text{Geometric}(p)\) 이면 \(\sum_{i=1}^r Y_i \sim \text{NegBin}(r, p)\) . 즉, \(r\) 개의 독립 기하 확률변수의 합이 음이항이다.

8 응용 분야

분야	활용	구체적 예시
역학	과산포 카운트	환자별 응급실 방문 횟수 (개인차로 인한 과산포)
생태학	역이항 표본 추출	목표 종 100마리 발견까지 조사한 총 개체 수
보험	청구 건수	계약자별 연간 보험금 청구 건수
유전학	성공까지 시행 수	원하는 유전형질 \(r\) 마리 확보까지 검사 수
웹 분석	페이지 조회 수	사용자별 세션당 페이지 조회 수 (이질성)
NLP	단어 빈도	문서별 특정 단어 출현 횟수 (문서 간 변동)

9 코드 예시

9.1 Step 1: 순수 Python 구현 (원리 이해)

import math

# ----- 음이항 PMF (실패 수 형태) -----
def negbinom_pmf(y, r, p):
    """P(Y = y | r, p) = C(r+y-1, y) * p^r * (1-p)^y"""
    return math.comb(r + y - 1, y) * p**r * (1 - p)**y

# 기본 성질
r, p = 5, 0.4
mean = r * (1 - p) / p
var = r * (1 - p) / p**2
print(f"NegBin(r={r}, p={p})")
print(f"E[Y] = {mean:.4f}")
print(f"Var(Y) = {var:.4f}")
print(f"Var/Mean = {var/mean:.4f} (과산포 정도)")

# 과산포 재모수화
mu = mean
print(f"\n재모수화: Var = mu + mu^2/r = {mu} + {mu**2/r} = {mu + mu**2/r}")

# PMF
print(f"\n{'y':>3} | {'P(Y=y)':>10} | {'누적':>10}")
cum = 0.0
for y in range(20):
    prob = negbinom_pmf(y, r, p)
    cum += prob
    print(f"{y:3d} | {prob:10.6f} | {cum:10.6f}")

# ----- 포아송 극한 확인 -----
lam = 5.0
print(f"\n음이항 → 포아송 극한 (lambda={lam}):")
print(f"{'r':>6} | {'p':>8} | {'E[Y]':>8} | {'Var(Y)':>8} | {'Var/E':>8}")
for r_val in [5, 10, 50, 100, 1000]:
    p_val = r_val / (r_val + lam)
    e = r_val * (1-p_val) / p_val
    v = r_val * (1-p_val) / p_val**2
    print(f"{r_val:6d} | {p_val:8.4f} | {e:8.4f} | {v:8.4f} | {v/e:8.4f}")

9.2 Step 2: scipy/numpy 구현 (실무 활용)

import numpy as np
from scipy import stats

# ----- scipy 음이항 분포 -----
# scipy: nbinom(n, p) where n=r (성공 횟수), p=성공 확률
r, p = 5, 0.4
dist = stats.nbinom(r, p)

print(f"E[Y] = {dist.mean():.4f} (이론: {r*(1-p)/p:.4f})")
print(f"Var(Y) = {dist.var():.4f} (이론: {r*(1-p)/p**2:.4f})")

# 포아송과 PMF 비교
lam = r * (1 - p) / p
poi = stats.poisson(lam)

print(f"\n{'y':>3} | {'NegBin':>10} | {'Poisson':>10}")
for y in range(15):
    print(f"{y:3d} | {dist.pmf(y):10.6f} | {poi.pmf(y):10.6f}")

# ----- 포아송 vs 음이항 회귀 비교 -----
import numpy as np
import statsmodels.api as sm

np.random.seed(42)
n = 300

# 과산포 데이터 생성 (음이항)
x = np.random.normal(0, 1, n)
r_true = 3
mu_true = np.exp(1.0 + 0.5 * x)
p_true = r_true / (r_true + mu_true)
y = stats.nbinom(r_true, p_true).rvs()

X = sm.add_constant(x)

# 포아송 회귀
pois_model = sm.GLM(y, X, family=sm.families.Poisson()).fit()
print("=== 포아송 회귀 ===")
print(f"피어슨 분산: {pois_model.pearson_chi2/pois_model.df_resid:.3f} "
      f"(1이면 등산포)")

# 음이항 회귀
nb_model = sm.GLM(y, X, family=sm.families.NegativeBinomial()).fit()
print(f"\n=== 음이항 회귀 ===")
print(f"피어슨 분산: {nb_model.pearson_chi2/nb_model.df_resid:.3f}")
print(f"\n포아송 AIC: {pois_model.aic:.1f}")
print(f"음이항 AIC: {nb_model.aic:.1f}")
print("(AIC가 작을수록 적합도 우수)")

# ----- 포아송-감마 혼합 시뮬레이션 -----
import numpy as np
from scipy import stats

np.random.seed(42)
n_sim = 100000
r, beta = 5, 1.5  # Gamma(r, beta)

# 방법 1: 감마에서 lambda 생성 -> 포아송 생성
lambdas = stats.gamma(a=r, scale=beta).rvs(n_sim)
poisson_gamma = np.array([stats.poisson(lam).rvs() for lam in lambdas])

# 방법 2: 직접 음이항 생성
p_equiv = 1 / (1 + beta)
negbin_direct = stats.nbinom(r, p_equiv).rvs(n_sim)

print("포아송-감마 혼합 vs 음이항 직접 생성:")
print(f"  혼합 — 평균: {poisson_gamma.mean():.3f}, 분산: {poisson_gamma.var():.3f}")
print(f"  직접 — 평균: {negbin_direct.mean():.3f}, 분산: {negbin_direct.var():.3f}")
print(f"  이론 — 평균: {r*beta:.3f}, 분산: {r*beta*(1+beta):.3f}")

9.3 R 코드

# R에서의 음이항 분포
r <- 5; p <- 0.4

# R의 dnbinom(x, size, prob)
cat("E[Y]:", r*(1-p)/p, "\n")
cat("Var(Y):", r*(1-p)/p^2, "\n")

# 시뮬레이션
set.seed(42)
samples <- rnbinom(10000, size=r, prob=p)
cat("표본 평균:", mean(samples), "\n")
cat("표본 분산:", var(samples), "\n")

# 음이항 회귀 (MASS 패키지)
library(MASS)
x <- rnorm(300)
y <- rnbinom(300, size=3, mu=exp(1 + 0.5*x))
model_nb <- glm.nb(y ~ x)
summary(model_nb)

10 관련 주제

선행 지식

베르누이 분포 — 음이항의 구성 단위
이항 분포 — 관점 전환의 출발점
기대값과 분산 — 적률 계산

후속 주제

기하 분포 — \(r = 1\) 특수 경우
포아송 분포 — \(r \to \infty\) 극한
이산 분포 개요 — 분포 간 관계 정리

관련 개념

GLM — 음이항 회귀, 로그 연결
MLE — 최대우도추정
GLMM 카운트 결과 — 과산포 혼합모형