[data analysis] 음이항분포(Negative Binomial Distribution)

음이항분포(Negative Binomial Distribution)

음이항 분포는 베르누이 시행을 반복하면서 r번째 성공까지의 확률의 변화를 나타내는 분포이며 파스칼 분포 (Pascal distribution)라고도 합니다. 예를 들어 동전의 앞면이 r회가 될 때까지 동전던지는 시행을 계속하는 경우로서 동전던지는 총 횟수가 확률변수 x가 됩니다. 이 확률변수는 동전 던지기 1회 당 앞면이 나오는 확률(p)가 앞면의 수(r)에 의존합니다. 이를 일반화하면 식 1과 같이 음이항분포의 확률변수는 성공횟수(r)와 베르누이 확률(p)을 매개변수로하는 확률분포로 나타낼 수 있습니다.

X ~ NB(r, p)

(식 1)

음이항 분포에서 최종 시행은 성공으로 시행이 종료되므로 식 2와 같이 마지막 시행을 제외한 경우에서 r - 1 횟수를 선택하는 이항확률로 간주할 수 있습니다(r은 성공횟수).

\begin{align} f(x) &= P(X = x)\\&=\binom{x-1}{r-1}p^r(1-p)^{x-r}\\ x&=r, r+1, \cdots \end{align}

(식 2)

식 2에서 성공횟수 r = 1인 경우는 기하분포가 됩니다.

예 1)

동전 던지기에서 4번 앞면이 발생하고 그 총 횟수가 10회가 되는 사건 A의 확률을 계산합니다.

사건 A의 확률은 식 3.1.32와 같이 나타낼 수 있습니다.

A = {X = 10}, 10번 시행에 4번 앞면	(식 3.1.32)
B = {X = 9}, 9번 시행에 3번 앞면
C = {X = 1}, 마지막 1번 시행이 앞면
P(A) = P(B ∩ C)

위 과정에서 P(C) = p(단일 시행에서 성공확률)가 되며 P(B)는 이항분포가 되므로 식 3과 같이 나타낼 수 있습니다.

$$P(B) = \binom{9}{3}p^3(1-p)^6$$

(식 3)

결국 A의 확률질량함수는 식 3에서 나타낸 함수에 최후의 성공확률을 고려한 것으로 식 4와 같이 계산됩니다.

\begin{align}P(B) &= \binom{9}{3}p^3(1-p)^6p\\&= \binom{9}{3}p^4(1-p)^6\end{align}

(식 4)

식 4를 일반화하면 총 횟수(n), 성공횟수(r), 시행당 확률(p)를 모수로하여 확률질량 함수는 식 5과 같이 나타낼 수 있습니다.

$$\text{pmf} = \binom{n-1}{r-1} p^r(1 − p)^{n-r}$$

(식 5)

다음 코드는 scipy.stats 모듈의 nbinom.pmf(x-r, r, p)를 적용하여 계산한 결과입니다.

p=1/2
p_A=round(stats.nbinom.pmf(6, 4, p), 3);p_A

0.082

round(special.comb(9, 3)*p**4*(1-p)**6, 3)

0.082

그림 1 모수 p에 따른 음이항분포의 변화.

x=range(3, 30)
p1=stats.nbinom.pmf(x, 3, 0.1)
p3=stats.nbinom.pmf(x, 3, 0.3)
p5=stats.nbinom.pmf(x, 3, 0.5)
p7=stats.nbinom.pmf(x, 3, 0.7)
fig, ax=plt.subplots(figsize=(4,3))
p=[p1,p3,p5,p7]
col=['g','b','r','k']
nme=[0.1, 0.3, 0.5, 0.7]
for i in range(len(p)):
    ax.plot(x, p[i], color=col[i], label=f"NB(3, {nme[i]})")
ax.set_xlabel("x")
ax.set_ylabel("Probability")
ax.legend(loc="best")
plt.show()

음이항분포는 기하분포의 합으로 나타낼 수 있으며 식 6과 같이 기대값과 분산을 계산할 수 있습니다.

\begin{align}X &=X_1+X_2+\cdots+X_r\\E(X)&= E(X_1)+E(X_2)+\cdots+E(X_r)\\&=r\frac{1}{p}\\ Var(X) &=Var(X_1)+Var(X_2)+\cdots+Var(X_r)\\&=r\frac{1-p}{p^2}\end{align}

(식 6)

예 2)

석유회사의 탐사 유정이 석유를 발견할 확률은 20%가 되어야 한다는 지질학적 연구를 근거로 다음을 결정합시다.

1. 세번째로 탐한 유정에서 석유가 처음 발견될 확률
2. 3개의 생산 유정을 설치하려는 경우 시추해야 하는 유정 수의 평균과 분산

a. 조건은 X = 3, r = 1, p = 0.2이므로 식 7과 같이 계산할 수 있습니다.

$$P(X = 3) = \binom{3-1}{1-1}0.2^1(1-0.2)^{3-1}$$

(식 7)

round(stats.nbinom.pmf(2, 1, 0.2), 3)

0.128

b. 식 3.1.36에 조건 r = 3, p = 0.2를 적용하여 계산합니다(식 8).

\begin{align}E(X) &= \frac{3}{0.2}=15\\ Var(x)&=3\frac{1-0.2}{0.2^2}=60 \end{align}

(식 8)

scipy.stats.nbinom 클래스의 다양한 메서드의 매개변수 loc은 데이터의 위치이동을 위해 지정합니다. 기본값은 0으로 표준화된 데이터로 인식합니다. 즉, x축 상에서 확률변수의 이동을 나타내지만 그에 대응하는 y축(각 x에 대한 확률)의 변화는 없습니다. 그러므로 pmf, cmf 등은 loc과는 관계없지만 기대값은 이 인자를 고려해 주어야 합니다. 다음 코드의 결과에서 기대값(평균)은 위 결과와 차이를 보입니다.

mu,var=stats.nbinom.stats(3, 0.2, moments="mv")
mu, var

(array(12.), array(60.))

이 차이는 loc=0 값에 의한 것입니다. 예제는 3개의 성공 유정을 위한 것이므로 위치를 3으로 이동하여야 합니다.

mu,var=stats.nbinom.stats(3, 0.2, loc=3, moments="mv")
mu, var

(array(15.), array(60.))