음이항분포¶

\(Y \sim \text{NegBin}(r = 3, p = 0.3)\)이므로

이다. 표준편차가 \(\sqrt{23.3} = 4.83\)으로 평균의 절반에 가깝다. 대기 문제는 원래 흔들림이 크다.

정확히 7번 만에 끝나려면 처음 6번 중에 성공이 2번, 7번째가 성공이어야 한다.

평균이 10인데 최빈값은 그보다 작다. 오른쪽으로 치우친 분포라 평균보다 빨리 끝나는 경우가 절반을 넘는다.

(1) 해석적으로. 실패 횟수 판본의 PMF를 먼저 적어 둔다. 정의 1은 시행 횟수 \(Y\)를 썼으나 SciPy 는 실패 횟수 \(K = Y - r\)을 쓰므로 \(k = m - r\)을 대입하면

다. 이웃한 두 확률의 비를 본다. \(q = 1-p\)라 두면

이다. \(\binom{k+r-1}{k}\big/\binom{k+r-2}{k-1} = \frac{(k+r-1)!}{(k+r-2)!}\cdot\frac{(k-1)!}{k!} = \frac{k+r-1}{k}\)를 썼고 \(p^r\)은 약분되었다. 이 비가 \(1\) 이상인 조건은

다. \(\frac{k+r-1}{k} = 1 + \frac{r-1}{k}\)는 \(k\)에 대해 감소하므로 비가 \(1\)을 지나는 자리는 한 곳뿐이고, 봉우리도 하나다. 따라서 문턱을 \(m = \dfrac{q(r-1)}{p}\)라 하면

이다. \(m\)이 정수이면 \(k = m\)에서 비가 정확히 \(1\)이 되어 \(p(m) = p(m-1)\)이므로 최빈값이 둘이다.

\(r = 1\)이면 \(m = 0\)이라 최빈값이 \(0\)이고 PMF가 처음부터 단조 감소한다. 기하분포로 환원되는 것이 식에서 그대로 읽힌다.

\(r = 5\), \(p = 0.4\)를 넣어 보면

으로 정수다. 곧 이 쪽의 기본 보기가 바로 동점인 경우이고, 최빈값은 \(5\)와 \(6\) 둘이다. 평균은 \(rq/p = 7.5\)라 최빈값보다 크다. 본문 문제에서 "평균이 10인데 최빈값은 그보다 작다"고 한 것과 같은 치우침이다.

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

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

# scipy의 음이항분포는 **실패 횟수**를 세는 판본이다.
# nbinom(r, p).pmf(k) = "r번째 성공이 나오기까지 실패가 k번 일어날 확률"
# 따라서 총 시행 횟수는 k + r 이고, x가 0부터 시작한다.
# 기하분포는 r = 1 인 특수한 경우다.
r, p = 5, 0.4
x = np.arange(0, 30)

fig, ax = plt.subplots(figsize=(12, 3))
# 기하분포와 달리 봉우리가 생긴다. 실패가 너무 적어도(운이 좋아도)
# 너무 많아도 확률이 낮기 때문이다.
ax.bar(x, stats.nbinom(r, p).pmf(x), alpha=0.7, label=f'NegBin(r={r}, p={p})')
ax.set_xlabel('k (number of failures before r-th success)')
ax.spines[['top', 'right']].set_visible(False)
ax.legend()
plt.show()

\(k = 5\)와 \(k = 6\)의 막대가 같은 높이로 나란히 서 있고 양옆으로 주저앉는다. 오른쪽 꼬리가 왼쪽보다 길다. 분수로 확인한다.

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

r, p = 5, F(2, 5)        # p = 0.4 를 분수로
q = 1 - p

# 비 p(k)/p(k-1) = q(k+r-1)/k
for k in range(1, 10):
    ratio = F(k + r - 1, k) * q
    mark = ">1 (오름)" if ratio > 1 else ("=1 (동점)" if ratio == 1 else "<1 (내림)")
    print(f"p({k})/p({k-1}) = {str(ratio):>6} = {float(ratio):.4f}   {mark}")

thr = q * (r - 1) / p
print(f"문턱 q(r-1)/p = {thr}  (정수인가? {thr.denominator == 1})")

def exact(k):            # C(k+r-1, k) p^r q^k
    return F(comb(k + r - 1, k)) * p**r * q**k

print(f"P(5) = {exact(5)},  P(6) = {exact(6)},  같은가? {exact(5) == exact(6)}")
rv = stats.nbinom(5, 0.4)
pm = rv.pmf(np.arange(0, 30))
print(f"부동소수점: pmf(5) - pmf(6) = {pm[5] - pm[6]:+.3e}")
print(f"동점으로 잡히는 자리 = {np.flatnonzero(pm == pm.max())}")
print(f"평균 rq/p = {float(r*q/p)},  P(K >= 30) = {rv.sf(29):.3e}")

출력:

p(1)/p(0) =      3 = 3.0000   >1 (오름)
p(2)/p(1) =    9/5 = 1.8000   >1 (오름)
p(3)/p(2) =    7/5 = 1.4000   >1 (오름)
p(4)/p(3) =    6/5 = 1.2000   >1 (오름)
p(5)/p(4) =  27/25 = 1.0800   >1 (오름)
p(6)/p(5) =      1 = 1.0000   =1 (동점)
p(7)/p(6) =  33/35 = 0.9429   <1 (내림)
p(8)/p(7) =   9/10 = 0.9000   <1 (내림)
p(9)/p(8) =  13/15 = 0.8667   <1 (내림)
문턱 q(r-1)/p = 6  (정수인가? True)
P(5) = 979776/9765625,  P(6) = 979776/9765625,  같은가? True
부동소수점: pmf(5) - pmf(6) = -5.551e-17
동점으로 잡히는 자리 = [6]
평균 rq/p = 7.5,  P(K >= 30) = 3.211e-04

비가 \(k = 5\)에서 \(27/25\)로 아직 \(1\)보다 크고 \(k = 6\)에서 정확히 \(1\)을 찍은 뒤 \(k = 7\)에서 \(33/35\)로 내려간다. \(27/25\)가 \(1\)에 아주 가깝다는 점에도 주의할 것. 봉우리가 \(k = 4\)부터 \(k = 7\)까지 거의 평평하다는 뜻이다.

분수로 재면 \(P(5)\)와 \(P(6)\)이 한 치도 다르지 않다(\(979776/9765625\)). 그런데 부동소수점으로는 \(5.6 \times 10^{-17}\)만큼 어긋나 flatnonzero가 \(6\) 하나만 집는다. 이항분포의 보기 1에서 본 것과 똑같은 일이다. \(p = 0.4\)가 이진수로 딱 떨어지지 않는다는 것 하나로 동점이 가려진다.

꼬리는 \(P(K \ge 30) = 3.21 \times 10^{-4}\)이다. 포아송의 \(3.45 \times 10^{-7}\)보다 세 자리 크지만 그림에서 버려도 좋을 만큼은 작다. 음이항분포의 꼬리는 기하분포와 마찬가지로 지수적으로만 줄어들어 두꺼운데, 평균 \(7.5\)에서 \(30\)은 표준편차 \(\sqrt{18.75} = 4.33\)으로 재어 \(5.2\) 표준편차 밖이다.

(1) 해석적으로. 최빈값은 보기 1의 식에 넣으면 된다. 문턱 \(m = q(r-1)/p\)를 \(q = 0.6\), \(p = 0.4\)에서 계산하면

이다. \(r = 1\)은 \(m = 0\)이라 오르는 구간이 없어 최빈값 \(0\) 하나이고 기하분포가 된다. \(r = 3\)은 \(m = 3\)이 정수라 최빈값이 \(\{2, 3\}\) 둘이다. \(r = 10\)은 \(m = 13.5\)가 정수가 아니라 \(\lfloor 13.5 \rfloor = 13\) 하나다.

왜도는 기하분포의 합이라는 구조에서 바로 나온다. \(K = \sum_{i=1}^r K_i\)이고 \(K_i\)가 독립이므로 3차 중심적률(= 3차 누적률)이 더해진다.

따라서

왜도는 \(1/\sqrt r\)로 줄어든다. 이것이 바로 중심극한정리가 작동하는 모습이다. 독립인 것을 \(r\)개 더하면 3차 누적률은 \(r\)배가 되지만 표준편차는 \(\sqrt r\)배가 되어 \(r^{3/2}\)로 나뉘므로, 비대칭이 \(r^{-1/2}\) 속도로 씻겨 나간다.

\(r = 1\)인 기하분포(실패 횟수)의 왜도는 \(\gamma_1 = (2-p)/\sqrt{q}\)이므로 \(p = 0.4\)에서

이다. 세 \(r\)에 넣으면 \(2.0656\), \(1.1926\), \(0.6532\)다. \(r\)을 열 배로 해야 왜도가 3분의 1이 된다. "\(r\)이 커지면 대칭에 가까워진다"는 말은 맞지만 그 속도가 매우 느리다는 것이 수가 말해 주는 바다. 왜도를 \(0.2\) 아래로 내리려면 \(r \ge (2.0656/0.2)^2 = 107\)이 필요하다.

(2) 해석적으로. 왜도가 정규근사 오차의 주항을 쥐고 있다. 격자분포에 연속성 수정을 넣으면 에지워스 전개의 첫 보정항이 남아

이다(\(\varphi\)는 표준정규밀도). 오차 \(\Phi(z) - P(K \le k)\)의 크기는 \(\frac{\gamma_1}{6}\lvert z^2-1 \rvert\varphi(z)\)이고, 이 함수는 \(z = 0\)에서 최대 \(\varphi(0) = 1/\sqrt{2\pi}\)를 갖는다(\(z = \pm\sqrt3\)에서 두 번째 봉우리가 있으나 \(2\varphi(\sqrt3) = 0.178\)로 더 작다). 따라서

이고, \(z = 0\)에서 \((z^2-1) < 0\)이므로 오차의 부호는 \(-\)다. 곧 최대오차가 평균 부근에서 나타나며 정규근사가 누적확률을 덜 잡는다. 세 \(r\)에 넣으면 \(-0.137\), \(-0.079\), \(-0.043\)으로 예측된다.

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

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

p = 0.4
x = np.arange(0, 40)

fig, ax = plt.subplots(figsize=(12, 3))
# r = 1 이면 기하분포다. 단조 감소하며 최빈값이 0이다.
# r이 커질수록 봉우리가 생기고 오른쪽으로 옮겨 가며 대칭에 가까워진다.
# 독립인 기하확률변수를 r개 더한 것이므로 중심극한정리가 작동한다.
for r in [1, 3, 10]:
    ax.plot(x, stats.nbinom(r, p).pmf(x), 'o-', markersize=4,
            label=f'r={r} (mean={r*(1-p)/p:.1f})')
ax.set_xlabel('k (failures before r-th success)')
ax.spines[['top', 'right']].set_visible(False)
ax.legend()
plt.show()

\(r = 1\) 곡선은 \(k = 0\)에서 출발해 단조 감소하고, \(r = 3\)은 낮은 봉우리가 왼쪽에 있으며 \(r = 10\)은 \(k = 13\) 근처에 종 모양 봉우리가 있다. 다만 세 곡선 모두 오른쪽 꼬리가 왼쪽보다 길다. \(r = 10\)조차 눈으로는 대칭처럼 보이지만 왜도가 \(0.65\)나 된다. 수로 확인한다.

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

p = F(2, 5)
q = 1 - p
k = np.arange(0, 2000)

print(f"{'r':>3}{'문턱 q(r-1)/p':>15}{'정수?':>7}{'최빈값(식)':>12}{'argmax':>8}"
      f"{'평균':>8}{'분산':>9}{'왜도':>9}{'왜도x√r':>10}")
for r in (1, 3, 10):
    thr = q * (r - 1) / p                  # 분수로 계산해야 정수 판정이 맞는다
    integral, m = thr.denominator == 1, int(thr)
    modes = f"{{{m-1}, {m}}}" if integral and m >= 1 else f"{m}"
    rv = stats.nbinom(r, float(p))
    skew = float(rv.stats("s"))
    print(f"{r:>3}{str(thr):>15}{str(integral):>7}{modes:>12}{rv.pmf(k).argmax():>8}"
          f"{float(r*q/p):>8.2f}{float(r*q/p**2):>9.2f}{skew:>9.4f}{skew*sqrt(r):>10.4f}")

# 동점은 분수로만 보인다. r = 3 에서 P(2) 와 P(3) 을 정확히 비교한다.
def exact(r, j):
    return F(comb(j + r - 1, j)) * p**r * q**j

for r, a, b in ((3, 2, 3), (5, 5, 6)):
    pm = stats.nbinom(r, float(p)).pmf(np.arange(0, 40))
    print(f"r={r}: P({a}) = {exact(r,a)},  P({b}) = {exact(r,b)},  같은가? {exact(r,a)==exact(r,b)}"
          f"   부동소수점 차 {pm[a]-pm[b]:+.3e}")

# 정규근사(연속성 수정)의 CDF 최대오차가 에지워스 예측과 맞는가.
print(f"{'r':>5}{'왜도 g':>9}{'최대오차':>11}{'-g/(6√(2π))':>14}{'오차/g':>10}")
for r in (1, 3, 10, 30, 100):
    rv = stats.nbinom(r, float(p))
    mu, sd = float(r*q/p), sqrt(float(r*q/p**2))
    e = stats.norm(mu, sd).cdf(k + 0.5) - rv.cdf(k)
    i = np.abs(e).argmax()
    g = float(rv.stats("s"))
    print(f"{r:>5}{g:>9.4f}{e[i]:>+11.4f}{-g/(6*sqrt(2*pi)):>+14.4f}{e[i]/g:>10.4f}")

출력:

  r    문턱 q(r-1)/p    정수?      최빈값(식)  argmax      평균       분산       왜도     왜도x√r
  1              0   True           0       0    1.50     3.75   2.0656    2.0656
  3              3   True      {2, 3}       2    4.50    11.25   1.1926    2.0656
 10           27/2  False          13      13   15.00    37.50   0.6532    2.0656
r=3: P(2) = 432/3125,  P(3) = 432/3125,  같은가? True   부동소수점 차 +1.110e-16
r=5: P(5) = 979776/9765625,  P(6) = 979776/9765625,  같은가? True   부동소수점 차 -5.551e-17
    r     왜도 g       최대오차   -g/(6√(2π))      오차/g
    1   2.0656    -0.1400       -0.1373   -0.0678
    3   1.1926    -0.0801       -0.0793   -0.0672
   10   0.6532    -0.0435       -0.0434   -0.0665
   30   0.3771    -0.0251       -0.0251   -0.0665
  100   0.2066    -0.0137       -0.0137   -0.0665

(1)과 (2)가 모두 맞는다.

왜도 \(\times \sqrt r\). 마지막 열이 세 줄 모두 정확히 \(2.0656\)이다. 왜도가 \(r^{-1/2}\)에 비례한다는 것을 소수점 넷째 자리까지 확인한 셈이고, \(r = 1\)의 왜도 \((2-p)/\sqrt q = 1.6/\sqrt{0.6}\)이 그 비례상수다.

최빈값과 동점. 식이 준 \(0\), \(\{2,3\}\), \(13\) 가운데 argmax는 \(0\), \(2\), \(13\)을 준다. \(r = 3\)에서 동점을 놓친 것인데, 분수로 재면 \(P(2) = P(3) = 432/3125\)로 정확히 같고 부동소수점 차가 \(1.1 \times 10^{-16}\)이다. 문턱을 float로 계산하면 더 나쁘다. \(0.6 \times 2 / 0.4\)가 배정도로는 \(2.9999999999999996\)이 되어 int()가 \(2\)를 돌려주고, 정수 판정 자체가 실패한다. 코드에서 Fraction을 쓴 이유가 그것이다.

정규근사 오차. 오차/g 열이 \(-0.0678\), \(-0.0672\), \(-0.0665\), \(-0.0665\), \(-0.0665\)로 예측한 \(-1/(6\sqrt{2\pi}) = -0.06649\)에 눕는다. \(r \ge 10\)에서는 소수점 넷째 자리까지 맞고 \(r = 1, 3\)에서 \(2\%\)쯤 어긋나는데, 버린 2차 보정항(\(\gamma_2\) 항)이 그만큼 일하기 때문이다. 부호가 모두 \(-\)인 것도 유도한 대로다.

쓸 만한가. \(r = 100\)에서도 CDF 최대오차가 \(0.0137\)이다. 이항분포가 \(np = 20\)에서 \(0.0039\)였던 것과 견주면 같은 "큰 표본"에서도 음이항 쪽이 네 배 가까이 나쁘다. 왜도가 \(1/\sqrt r\)로만 줄기 때문이고, 과대산포 자료에 정규근사를 선뜻 들이대면 안 되는 까닭이다. 대신 에지워스 보정항을 더하거나 아예 정확한 음이항 CDF를 쓰는 편이 낫다.

(1) 해석적으로. 평균 쪽은 쉽다. \(\text{Var}(K) = rq/p^2 = 5 \times 0.6/0.16 = 18.75\)이므로 \(N = 2 \times 10^5\)에서

이다. \(P(K = 10)\) 쪽도 쉽다. 표본비율이므로 \(\theta = P(K=10) = 0.0620\)에서

이다. 분산 쪽이 일이 많다. \(\text{Var}(S^2) \approx (\mu_4 - \sigma^4)/N\)이고 음이항분포의 \(\mu_4\)가 필요한데, 닫힌 꼴을 외우는 대신 누적률로 간다.

\(K = \sum_{i=1}^r K_i\)이고 \(K_i\)가 독립이므로 누적률이 더해진다.

그리고 4차 중심적률과 누적률의 관계는 \(\mu_4 = \kappa_4 + 3\kappa_2^2\)이므로

이다. 이항분포의 보기 3에서 쓴 것과 똑같은 공식인데, 거기서는 베르누이를 \(n\)개 더했고 여기서는 기하를 \(r\)개 더했을 뿐이다. 독립인 것을 더하는 구조는 분포가 무엇이든 같은 길을 낸다.

남은 것은 기하분포(실패 횟수) 한 토막의 중심적률이다. \(p = 0.4\), \(q = 0.6\)에서

이므로

이다. 이것을 넣으면

이다. 분산이 평균보다 일곱 배 넘게 더 흔들린다. 꼬리가 두꺼운 분포라서 제곱의 평균이 특히 불안하다.

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

import numpy as np
from scipy import stats

np.random.seed(42)
r, p = 5, 0.4

# 정의를 그대로 실행한다. 기하분포(시행 번호 판본) r개를 더한다.
# scipy의 geom은 1부터 시작하는 시행 번호를 주므로 합의 최솟값이 r 이다.
trials = stats.geom(p).rvs(size=(r, 200_000)).sum(axis=0)
failures = trials - r          # scipy의 nbinom과 맞추려면 실패 횟수로 바꾼다

print(f"평균  : 표본 {failures.mean():.4f},  이론 {r*(1-p)/p:.4f}")
print(f"분산  : 표본 {failures.var():.4f},  이론 {r*(1-p)/p**2:.4f}")
print(f"P(K=10): 표본 {np.mean(failures == 10):.4f},  이론 {stats.nbinom(r, p).pmf(10):.4f}")

출력:

평균  : 표본 7.5096,  이론 7.5000
분산  : 표본 18.8118,  이론 18.7500
P(K=10): 표본 0.0624,  이론 0.0620

기하분포 \(r\)개의 합이 음이항분포라는 사실이 세 수치에서 모두 확인된다.

이 "확인"이 얼마나 단단한지 (1)의 표준오차로 재 본다. 유도한 \(\mu_4\)도 열거로 검산한다.

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

r, p, Nsim = 5, 0.4, 200_000
q = 1 - p
mean, var = r * q / p, r * q / p**2

# 기하(실패 횟수) 한 토막의 중심적률. 열거로 정확히 구하고 식과 맞춘다.
k = np.arange(0, 600)
g = stats.nbinom(1, p).pmf(k)
mg = (k * g).sum()
mu2_g, mu4_g = (((k - mg)**2) * g).sum(), (((k - mg)**4) * g).sum()
print(f"기하: mu2 = {mu2_g:.4f} (q/p^2 = {q/p**2:.4f}),  "
      f"mu4 = {mu4_g:.4f} (q(1+7q+q^2)/p^4 = {q*(1+7*q+q*q)/p**4:.4f})")
kappa4_g = mu4_g - 3 * mu2_g**2
print(f"기하: kappa4 = mu4 - 3 mu2^2 = {kappa4_g:.4f}")

# 누적률은 더해진다. K = 기하 r 개의 합.
mu4_K = r * kappa4_g + 3 * (r * mu2_g)**2
mu4_enum = (((k - mean)**4) * stats.nbinom(r, p).pmf(k)).sum()
print(f"NB:   mu4 = r*kappa4 + 3(r*mu2)^2 = {mu4_K:.4f}   열거 {mu4_enum:.4f}")
print(f"NB:   mu4 - sigma^4 = {mu4_K - var**2:.4f}")

np.random.seed(42)
trials = stats.geom(p).rvs(size=(r, Nsim)).sum(axis=0)
failures = trials - r

theta = stats.nbinom(r, p).pmf(10)
se_mean = sqrt(var / Nsim)
se_var = sqrt((mu4_K - var**2) / Nsim)
se_pmf = sqrt(theta * (1 - theta) / Nsim)
for name, value, theo, se in (("평균   ", failures.mean(), mean, se_mean),
                              ("분산   ", failures.var(), var, se_var),
                              ("P(K=10)", np.mean(failures == 10), theta, se_pmf)):
    print(f"{name} 표본 {value:.4f}   이론 {theo:.4f}   SE {se:.4f}   z = {(value - theo) / se:+.3f}")

출력:

기하: mu2 = 3.7500 (q/p^2 = 3.7500),  mu4 = 130.3125 (q(1+7q+q^2)/p^4 = 130.3125)
기하: kappa4 = mu4 - 3 mu2^2 = 88.1250
NB:   mu4 = r*kappa4 + 3(r*mu2)^2 = 1495.3125   열거 1495.3125
NB:   mu4 - sigma^4 = 1143.7500
평균    표본 7.5096   이론 7.5000   SE 0.0097   z = +0.995
분산    표본 18.8118   이론 18.7500   SE 0.0756   z = +0.818
P(K=10) 표본 0.0624   이론 0.0620   SE 0.0005   z = +0.753

(1)이 다 맞는다.

\(\mu_4\)가 두 길에서 같다. 누적률 가법성으로 계산한 \(1495.3125\)와 PMF를 \(k = 0\)부터 \(599\)까지 열거해 더한 \(1495.3125\)가 소수점 넷째 자리까지 같다. 기하 한 토막의 \(\mu_4\)도 열거값과 식 \(q(1+7q+q^2)/p^4\)가 일치한다. 유도가 맞았다는 확인이 이것이다. 열거 쪽은 꼬리를 자른 근사이지만 \(P(K \ge 600)\)이 배정도 밑으로 떨어질 만큼 작아 차이가 드러나지 않는다.

세 \(z\). \(+1.00\), \(+0.82\), \(+0.75\)다. 모두 \(1\) SE 근처이니 몬테카를로 오차 범위 안이고, 세 값이 다 양수인 것도 \(2 \times 10^5\)개 표본 하나에서 나온 값들이라 서로 상관되어 있으니 이상하지 않다. 셋 다 같은 방향으로 조금 큰 쪽에 앉은 표본을 받았을 뿐이다.

날 것의 차이는 오해를 부른다. 분산의 차이 \(18.8118 - 18.75 = 0.0618\)은 평균의 차이 \(0.0096\)보다 여섯 배 크다. 그러나 표준오차로 나누면 \(+0.82\) 대 \(+1.00\)으로 분산 쪽이 오히려 더 잘 맞은 것이다. 꼬리가 두꺼운 분포에서 분산의 표준오차가 크다는 것을 모르면 거꾸로 읽게 된다.

\(Z = k\)이려면 처음 \(k - 1\)번의 시행에서 정확히 \(r - 1\)번 성공하고 \(k\)번째 시행에서 성공해야 한다:

\(Z\)는 독립인 기하 확률변수 \(Y_1, \ldots, Y_r\)(각 성공까지의 대기 시간)의 합으로 쓸 수 있다. 따라서:

\(r = 1\)이면 기하분포로 환원된다. 음이항분포는 기하분포를 여러 번의 성공으로 일반화하며, 과대산포된 계수 자료의 모형에 바탕이 된다.

조건부 확률에 사전밀도를 곱해 적분한다.

남은 적분은 감마적분이므로 \(\Gamma(y+r)\,(1+1/\theta)^{-(y+r)}\)이다. \(p = 1/(1+\theta)\)로 두면 \(1-p = \theta/(1+\theta)\)이고, 정리하면

를 얻는다. 이것이 (실패 횟수를 세는 판본의) 음이항분포이다. \(\square\)

\(r\)이 정수일 필요도 없다. 이항계수 대신 감마함수로 쓴 덕분에 \(r > 0\)인 실수면 된다.

왜 중요한가. 적률을 계산하면

이다. 분산이 평균보다 \((1+\theta)\)배 크다. 포아송분포는 평균과 분산이 같아야 하는데, 실제 계수 자료는 거의 언제나 분산이 더 크다. 이를 과대산포라 한다.

이 유도가 그 원인을 말해 준다. 개체마다 사건 발생률 \(\lambda\)가 다르면(관측되지 않은 이질성), 전체를 뭉뚱그린 분포는 포아송이 아니라 그 혼합이 되고 분산이 부풀어 오른다. 보험 가입자마다 사고 성향이 다르고, 지역마다 감염 위험이 다르며, 유전자마다 발현량이 다르다.

실무적 귀결은 분명하다. 과대산포된 자료에 포아송 회귀를 쓰면 계수 추정은 그런대로 나와도 표준오차가 심하게 과소평가되어 있지도 않은 유의성이 쏟아진다. 음이항 회귀를 쓰면 \(\theta\)가 이 여분의 산포를 흡수한다. RNA 시퀀싱 자료 분석 도구들이 하나같이 음이항 모형을 쓰는 이유이며, \(\theta \to 0\)이면 포아송으로 돌아가므로 포아송을 특수한 경우로 포함한다.

시행 횟수 판본을 \(Y\), 실패 횟수 판본을 \(K\)라 하면 \(Y = K + r\)이다. 상수만큼 평행이동한 것이므로 분산은 같고 평균만 \(r\) 차이 난다.

SciPy는 실패 횟수 판본이다. stats.nbinom(r, p).pmf(k)가 "\(r\)번째 성공까지 실패가 \(k\)번"일 확률을 준다. mean()을 찍어 보면 \(r/p\)가 아니라 \(r(1-p)/p\)가 나오는 것으로 곧바로 확인된다. \(r = 5\), \(p = 0.4\)이면 12.5가 아니라 7.5다.

기하분포에서 stats.geom이 시행 번호 판본이었던 것과 반대라는 점이 특히 헷갈린다. 같은 라이브러리 안에서 규약이 다르므로, 두 분포를 함께 쓰는 코드에서는 \(r=1\)일 때 geom과 nbinom의 평균이 1만큼 어긋난다는 사실을 기억해야 한다.

from scipy import stats
print(stats.geom(0.4).mean())      # 2.5   = 1/p        (시행 번호)
print(stats.nbinom(1, 0.4).mean()) # 1.5   = (1-p)/p    (실패 횟수)

\(q = 1-p\)로 두고 더한다.

여기서 음이항급수

를 쓰면 합이 \(p^r(1-q)^{-r} = p^r p^{-r} = 1\)이 된다. \(\square\)

이름의 유래가 바로 이 급수다. 뉴턴의 이항정리를 음의 지수 \(-r\)로 확장하면

이므로 PMF를 \(\binom{-r}{k}p^r(-q)^k\)로도 쓸 수 있다. 이항계수의 위 첨자가 음수라는 데서 "음이항"이 왔다.

이 표현은 형식적인 이름짓기에 그치지 않는다. \(r\)이 정수일 필요가 없다는 사실이 여기서 따라 나온다. 이항계수를 감마함수로 쓰면 \(r > 0\)인 어떤 실수에서도 PMF가 잘 정의되며, 이것이 계수 자료 회귀에서 \(r\)을 연속 모수로 추정하는 근거다(연습문제 2의 혼합 표현과 같은 이야기다).

\(q = 1-p\)이고 평균이 \(rq/p = \lambda\)로 고정되어 있으므로 \(q = \lambda/(r+\lambda)\), \(p = r/(r+\lambda)\)로 두면 된다. 고정된 \(k\)에 대해

이다. 세 조각으로 나누어 극한을 취한다.

• (가)의 분자는 \(r\)에 대한 \(k\)차 다항식이므로 \(r^k\)로 나누면 1로 간다. 즉 (가) \(\approx r^k/k!\).
• \(\left(\frac{\lambda}{r+\lambda}\right)^k \cdot r^k = \lambda^k\left(\frac{r}{r+\lambda}\right)^k \to \lambda^k\).
• \(\left(\frac{r}{r+\lambda}\right)^r = \left(1 + \frac\lambda r\right)^{-r} \to e^{-\lambda}\).

셋을 곱하면

이다. \(\square\)

뜻. 성공이 거의 확실해지면(\(p \to 1\)) 실패는 드문 사건이 되고, 드문 사건의 개수는 포아송을 따른다. 이항분포에서 \(p \to 0\)으로 포아송에 닿았던 것과 대칭을 이루는 길이다.

모수로 보면 이 극한은 \(r \to \infty\)에서 과대산포가 사라지는 과정이기도 하다. 음이항의 분산은 평균의 \(1/p = 1 + \lambda/r\)배인데, \(r \to \infty\)면 이 배율이 1로 가서 포아송의 "평균 = 분산"이 회복된다.

판단. 포아송이면 평균과 분산이 같아야 하는데 \(s^2/\bar x = 11.8/4.2 = 2.81\)로 세 배 가깝다. 뚜렷한 과대산포이므로 포아송은 부적절하다.

적률추정. 실패 횟수 판본에서 평균과 분산은

이다. 둘째 식에서 곧바로

이고, 첫째 식에 넣으면

이다. \(\hat r\)이 정수가 아니어도 된다는 것은 연습문제 4에서 보았다.

해석. 혼합 관점(연습문제 2)으로 읽으면 \(\hat r = 2.32\)가 이질성의 크기를 말해 준다. \(r\)이 작을수록 날짜별 발생률의 편차가 크다는 뜻이고, \(r \to \infty\)면 포아송으로 돌아간다. 여기서는 \(r\)이 작으니 "어떤 날은 조용하고 어떤 날은 폭주하는" 자료다.

다만 과대산포의 원인을 모형으로 덮기 전에 먼저 확인할 것이 있다. 요일 효과나 사건(제품 출시, 장애)처럼 설명 가능한 변동이 있으면 그것을 설명변수로 넣는 편이 낫다. 설명변수를 넣고도 남는 산포가 있을 때 음이항이 제 역할을 한다.

생성함수. \(Z\)는 독립인 기하확률변수 \(Y_1, \ldots, Y_r\)의 합이고(연습문제 1), 기하분포의 확률생성함수는 \(ps/(1-qs)\)이다(기하분포 문서 연습문제 8). 독립인 합에서는 생성함수가 곱해지므로

이다. \(q = 1-p\)로 썼다.

(a) 적률. \(\log G_Z(s) = r\log p + r\log s - r\log(1-qs)\)를 미분하면

이고 \(s = 1\)에서 \(G_Z(1) = 1\)이므로

이다. 한 번 더 미분해 계승적률 \(E[Z(Z-1)] = G_Z''(1)\)을 구하고 정리하면 \(\operatorname{Var}(Z) = rq/p^2\)을 얻는다. 기하분포(\(r=1\))의 값에 정확히 \(r\)을 곱한 것으로, 독립인 합이므로 당연한 결과다.

(b) 가법성. \(Z_1 \sim \text{NB}(r_1, p)\), \(Z_2 \sim \text{NB}(r_2, p)\)가 독립이면

이므로 \(Z_1 + Z_2 \sim \text{NB}(r_1+r_2,\, p)\)이다. \(\square\)

\(p\)가 같아야 한다는 조건이 핵심이다. \(p_1 \ne p_2\)이면 지수를 합칠 수 없고, 합은 음이항분포가 아니다. 이항분포가 \(p\)가 같을 때만 더해지는 것과 같은 사정이다.

정수가 아니어도 성립한다. 위 계산은 \(r_1, r_2\)가 정수라는 사실을 어디에서도 쓰지 않았다. 지수법칙 \(a^{r_1}a^{r_2} = a^{r_1+r_2}\)만 썼을 뿐이다. 연습문제 4에서 본 대로 \(r > 0\)인 실수에서 PMF가 잘 정의되므로, 가법성은 실수 \(r\)에서 그대로 성립한다.

이것이 실무에서 쓰이는 방식이 있다. 계수 자료를 합칠 때 모형이 유지된다는 뜻이기 때문이다. 일별 사고 건수가 \(\text{NB}(r, p)\)이면 주간 합계는 \(\text{NB}(7r, p)\)다. 단, 날짜들이 독립이고 \(p\)가 같아야 한다. 포아송분포가 \(\text{Poisson}(\lambda_1) + \text{Poisson}(\lambda_2) = \text{Poisson}(\lambda_1+\lambda_2)\)로 닫혀 있는 것의 과대산포판이다.

두 사건이 같은 사건이다. 같은 베르누이 시행열을 놓고 두 확률변수를 함께 정의하자. \(Z\)는 \(r\)번째 성공이 일어나는 시행 번호이고, \(X\)는 처음 \(n\)번의 시행에서 나온 성공 횟수다. 그러면

이다. 같은 표본점의 집합이므로 확률도 같다. \(\square\)

계산이 전혀 없다는 점이 요점이다. 두 분포의 PMF를 더해 맞추려 들면 이항계수 항등식과 씨름하게 되지만, "\(r\)번째 성공이 \(n\)번 안에 일어난다"와 "\(n\)번 안에 성공이 \(r\)번 이상 나온다"가 말만 다른 같은 문장임을 알아채면 한 줄로 끝난다. 확률론에서 자주 통하는 수법이다.

import numpy as np
from scipy import stats

r, p = 4, 0.35
print("쌍대성 확인:  P(Z<=n) = P(X>=r),  X~B(n,p)")
print(f"{'n':>4}{'P(Z<=n) [nbinom]':>20}{'P(X>=r) [binom]':>19}")
for n in (4, 6, 10, 15, 25):
    lhs = stats.nbinom(r, p).cdf(n - r)          # 실패 k=n-r 이하
    rhs = 1 - stats.binom(n, p).cdf(r - 1)
    print(f"{n:>4}{lhs:>20.8f}{rhs:>19.8f}")

출력:

쌍대성 확인:  P(Z<=n) = P(X>=r),  X~B(n,p)
   n    P(Z<=n) [nbinom]    P(X>=r) [binom]
   4          0.01500625         0.01500625
   6          0.11742391         0.11742391
  10          0.48617298         0.48617298
  15          0.82730351         0.82730351
  25          0.99031529         0.99031529

소수점 여덟 자리까지 일치한다. nbinom이 실패 횟수 판본이므로 \(Z \le n\)을 실패 \(\le n-r\)로 옮겨 넣었다(연습문제 3).

어디에 쓰이는가.

• 표본크기 설계. "성공 사례 \(r\)건을 모으려는데 \(n\)번 안에 끝날 확률"을 묻는 임상시험·품질검사 문제가 이항분포의 꼬리확률로 바뀐다. 익숙한 쪽으로 옮겨서 풀 수 있다.
• 축차검정(중간 분석). 목표 사건 수에 도달할 때까지 관측하는 설계에서, 정해진 \(n\)에서 멈추는 설계와 확률을 견줄 때 이 항등식이 다리가 된다.
• 수치 계산. 한쪽 분포의 꼬리가 수치적으로 불안정할 때 다른 쪽으로 바꿔 계산한다.

더 중요한 개념적 함의가 있다. 같은 자료를 "\(n\)을 고정하고 성공 수를 본다"(이항)와 "\(r\)을 고정하고 시행 수를 본다"(음이항)로 볼 수 있으며, 이 둘은 멈추는 규칙만 다르다. 확률은 위처럼 맞아떨어지지만 추론은 달라질 수 있다. 9장에서 다룰 정지규칙 문제이며, 같은 자료에 대해 두 설계가 다른 \(p\)값을 주는 유명한 예가 여기서 나온다.

불편성 증명. 정의대로 더한다.

계수를 정리하는 것이 전부다.

이므로

이다. \(j = k-1\)로 바꾸면 \(j\)는 \(r-1\)부터 시작하고 \(k-r = j-(r-1)\)이므로

이다. 마지막 합이 \(\text{NB}(r-1, p)\)의 확률을 모두 더한 것이라 \(1\)이 된다. \(\square\)

\(r/Z\)의 편향은 위쪽이다. \(g(z) = r/z\)는 볼록함수이므로 옌센 부등식에 의해

이다(등호는 \(Z\)가 상수일 때만). 즉 \(r/Z\)는 \(p\)를 체계적으로 과대추정한다.

import numpy as np
rng = np.random.default_rng(0)
p, reps = 0.3, 500_000
print(f"참값 p = {p}\n")
print(f"{'r':>4}{'E[(r-1)/(Z-1)]':>18}{'E[r/Z]':>12}{'편향':>10}")
for r in (2, 3, 5, 10):
    Z = rng.negative_binomial(r, p, reps) + r      # 시행 횟수 판본
    unb = np.mean((r - 1) / (Z - 1))
    naive = np.mean(r / Z)
    print(f"{r:>4}{unb:>18.5f}{naive:>12.5f}{naive - p:>+10.5f}")

출력:

참값 p = 0.3

   r    E[(r-1)/(Z-1)]      E[r/Z]        편향
   2           0.29989     0.41479  +0.11479
   3           0.30023     0.37638  +0.07638
   5           0.30010     0.34487  +0.04487
  10           0.30006     0.32183  +0.02183

\(r = 2\)에서 \(r/Z\)의 편향이 \(+0.115\)나 된다. 참값 \(0.3\)을 \(0.41\)로 읽는 것이니 38% 과대추정이다. \(r\)이 커지면 편향이 줄지만 사라지지는 않는다. 반면 \(\tilde p\)는 모든 \(r\)에서 소수점 셋째 자리까지 \(0.300\)이다.

\(r \ge 2\)가 필요한 이유. \(r = 1\)이면 \(\tilde p = 0/(Z-1) = 0\)으로 쓸모가 없다. \(r = 1\)(기하분포)에서 \(p\)의 불편추정량이 없는 것은 아니지만, 쓸 만한 것이 없다. \(E[\delta(Z)] = \sum_{k\ge1}\delta(k)p\,q^{k-1} = p\)가 모든 \(p\)에서 성립하려면 멱급수의 계수를 맞추어야 하므로 \(\delta(1) = 1\), \(\delta(k) = 0\ (k \ge 2)\)뿐이다. 즉 유일한 불편추정량은 지시함수 \(\mathbb{1}\{Z = 1\}\)이고, 이것은 \(0\) 아니면 \(1\)만 내놓는다. 불편성만으로는 좋은 추정량이 보장되지 않는다는 것을 보여 주는 표준적인 예이며, 6장의 추정량 품질 논의와 이어진다.

역표본추출이 쓰이는 자리. 희귀사건의 비율을 추정할 때 \(n\)을 먼저 정하면 성공이 하나도 안 나와 \(\hat p = 0\)이 되는 사고가 난다. "성공 \(r\)건이 모일 때까지" 관측하면 그 일이 원천적으로 일어나지 않고, 덤으로 위의 불편추정량까지 얻는다. 감염병 감시나 희귀 결함 검사에서 쓰는 설계이며, 연습문제 8의 정지규칙 논의와 바로 이어진다.

(a) \(r\)을 알 때. 로그가능도는 상수항을 빼면

이고, 미분해 \(0\)으로 두면

를 얻는다.

적률추정량과 같다. 실패 횟수 판본의 평균이 \(E[K] = r(1-p)/p\)이므로, 적률법은 \(\bar k = r(1-\hat p)/\hat p\)를 풀어 \(\hat p = r/(r+\bar k)\)를 준다. 정확히 같은 식이다. \(r\)이 알려져 있으면 모수가 하나뿐이고 \(\bar k\)가 완전충분통계량이라, 두 방법이 일치한다.

(b) \(r\)도 모를 때. 이제 \(\hat p = r/(r+\bar k)\)를 로그가능도에 되넣어 \(r\)만의 함수로 만든다(프로파일 가능도). PMF를 감마함수로 쓰면

이고 \(r\)로 미분하면 디감마함수 \(\psi = (\log\Gamma)'\)가 나온다.

즉 우도방정식은

이다.

왜 닫힌 형태가 없는가. \(\psi\)는 초등함수로 표현되지 않고, 더구나 \(\psi(k_i + r)\)처럼 자료마다 다른 점에서 평가된 \(\psi\)의 합이 들어 있다. \(r\)을 한쪽으로 분리할 방법이 없다. 그래서 실무에서는 적률추정량을 출발점으로 수치적으로 푼다.

import numpy as np
from scipy import stats, optimize, special

rng = np.random.default_rng(1)
r_true, p_true, n = 3.0, 0.4, 500
K = rng.negative_binomial(r_true, p_true, n)     # 실패 횟수 판본
kbar, s2 = K.mean(), K.var(ddof=1)

# r 을 안다고 할 때: MLE = 적률법
p_known = r_true / (r_true + kbar)
print(f"r 을 알 때   MLE p_hat = {p_known:.4f}   적률법 p_hat = {r_true/(r_true+kbar):.4f}  (동일)")

# r 도 모를 때
p_mom = kbar / s2
r_mom = kbar * p_mom / (1 - p_mom)
def negll(theta):
    r, p = theta
    if r <= 0 or not (0 < p < 1): return np.inf
    return -np.sum(stats.nbinom.logpmf(K, r, p))
res = optimize.minimize(negll, [r_mom, p_mom], method='Nelder-Mead')
r_mle, p_mle = res.x
print(f"\n참값        r = {r_true:.4f}   p = {p_true:.4f}")
print(f"적률법      r = {r_mom:.4f}   p = {p_mom:.4f}")
print(f"최대가능도  r = {r_mle:.4f}   p = {p_mle:.4f}")
print(f"\nr 의 우도방정식 잔차: "
      f"{np.mean(special.digamma(K+r_mle))-special.digamma(r_mle)+np.log(p_mle):.2e}")

출력:

r 을 알 때   MLE p_hat = 0.3935   적률법 p_hat = 0.3935  (동일)

참값        r = 3.0000   p = 0.4000
적률법      r = 3.1492   p = 0.4051
최대가능도  r = 3.3411   p = 0.4195

r 의 우도방정식 잔차: 2.10e-06

\(r\)을 알 때는 두 방법이 소수점 넷째 자리까지 같고, 모를 때는 갈린다. 잔차가 \(10^{-6}\) 수준인 것은 수치해가 위 우도방정식을 실제로 만족한다는 확인이다.

\(r\)은 추정하기 어려운 모수다. 자료에 담긴 \(r\)에 대한 정보가 적어서, \(n = 500\)인데도 참값 \(3.0\)에 대해 적률법 \(3.15\), 최대가능도 \(3.34\)로 제법 흩어진다. 꼬리의 모양을 결정하는 모수라 극단값 몇 개에 크게 휘둘리기 때문이다. 실무에서 음이항 회귀의 산포모수 신뢰구간이 늘 넓은 것이 이 때문이며, 표본이 작으면 \(r\)을 아예 고정하고 쓰는 편이 안정적일 때도 있다.

어느 쪽을 쓸 것인가. 최대가능도가 점근적으로 효율적이므로 기본은 최대가능도다. 다만 적률법은 닫힌 형태라 계산이 즉시 끝나고 수렴 실패가 없으므로, 최대가능도의 출발점으로 쓰는 것이 위 코드가 보여 주는 표준적인 조합이다.