コインを繰り返し投げて、表が3回出るまでに何回投げることになるでしょうか。5回で終わることもあれば、20回以上かかることもあります。この「成功がちょうど $r$ 回に達するまでの試行回数」が従う分布が負の二項分布（negative binomial distribution）です。
負の二項分布は名前こそ馴染みにくいですが、実は非常に広い応用範囲を持っています。特に近年では、ポアソン分布では説明できない過分散（overdispersion）を持つカウントデータのモデリングにおいて不可欠な存在です。
負の二項分布を理解すると、以下のような応用が開けます。
- 生態学: 生物の空間分布（個体が集中して分布する傾向）のモデリング
- 保険数理: 事故件数のモデリング。ポアソン分布より柔軟に過分散を扱える
- RNA-seq解析: 遺伝子発現量のカウントデータは負の二項分布でモデル化される（DESeq2, edgeRなど）
- 自然言語処理: 文書中の単語の出現回数のモデリング
本記事では、負の二項分布が「幾何分布の和」として自然に導入されることを直感的に理解し、確率質量関数を導出します。さらに、ポアソン・ガンマ混合としての別表現を示し、過分散データへの応用をPythonで実装します。
本記事の内容
- 負の二項分布の直感的理解と数学的定義
- 確率質量関数の導出と二項係数の「負」の意味
- モーメント母関数、期待値、分散の計算
- ポアソン・ガンマ混合としての解釈
- 二項分布・ポアソン分布・幾何分布との関係
- Pythonによる可視化と過分散データへの応用
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 二項分布 — 負の二項分布と「逆の関係」にある分布
- 幾何分布 — 負の二項分布の $r = 1$ の場合
- ポアソン分布 — 極限として、また混合分布の基盤として
- ガンマ分布 — ポアソン・ガンマ混合の理解に必要
負の二項分布とは — 直感的な理解
幾何分布の拡張
幾何分布は「最初の成功が起きるまでの試行回数」の分布でした。成功確率 $p$ のベルヌーイ試行を繰り返すとき、初めて成功するまでに必要な試行回数 $X$ は幾何分布 $\text{Geom}(p)$ に従います。
負の二項分布は、この幾何分布を「$r$ 回目の成功まで」に一般化したものです。
たとえば、野球のバッターが打率（成功確率）$p = 0.3$ で打席に立つとします。「3本ヒットを打つまでに何回打席に立つか」という問いの答えは、$r = 3$, $p = 0.3$ の負の二項分布に従います。
最速では3打席（3打席連続ヒット）で達成できますが、ヒットとアウトが混ざるため、通常はもっと多くの打席を要します。平均的には $r/p = 3/0.3 = 10$ 打席ほどかかると期待されます。
二項分布との「逆の関係」
負の二項分布の名前に「負の」が付く理由は、二項分布と「逆の問い」に答えているからです。
- 二項分布: 試行回数 $n$ を固定し、成功回数 $X$ の分布を求める → $X \sim \text{Bin}(n, p)$
- 負の二項分布: 成功回数 $r$ を固定し、必要な試行回数 $Y$ の分布を求める → $Y \sim \text{NB}(r, p)$
つまり、「何が固定で何が確率変数か」が逆転しているのです。この逆転の関係は、確率質量関数の導出でも明確に現れます。
パラメータ化の注意
負の二項分布には教科書によって異なるパラメータ化が存在します。本記事では以下の定義を使います。
- $Y$: $r$ 回目の成功が起きるまでの総試行回数（$Y \geq r$）
- $X = Y – r$: $r$ 回目の成功が起きるまでの失敗回数（$X \geq 0$）
$Y$ で定義する流儀と $X$ で定義する流儀があります。本記事では失敗回数 $X$ を確率変数とするパラメータ化を主に使いますが、両方の表現を示します。
この直感的な理解を踏まえて、数学的定義に進みましょう。
負の二項分布の数学的定義
成功確率 $p$（$0 < p \leq 1$）のベルヌーイ試行を独立に繰り返すとき、$r$ 回目の成功が起きるまでの失敗回数 $X$ が従う分布を、パラメータ $(r, p)$ の負の二項分布といいます。
$X$ の確率質量関数（PMF）は次の式で与えられます。
$$ \begin{equation} P(X = x) = \binom{x + r – 1}{x}\,p^r\,(1-p)^x, \quad x = 0, 1, 2, \dots \end{equation} $$
ここで $\binom{x+r-1}{x} = \frac{(x+r-1)!}{x!\,(r-1)!}$ は二項係数です。
この式は「$r$ 回の成功を達成するために、ちょうど $x$ 回失敗する確率」を表しています。式の各因子の意味は次の通りです。
- $p^r$: $r$ 回の成功の確率
- $(1-p)^x$: $x$ 回の失敗の確率
- $\binom{x+r-1}{x}$: 全 $x + r$ 回の試行のうち、最後の1回が必ず成功（$r$ 回目の成功）であり、残りの $x + r – 1$ 回の中に $x$ 回の失敗と $r – 1$ 回の成功がどのように配置されるかの組み合わせ数
総試行回数 $Y = X + r$ で表す場合の確率質量関数は次のようになります。
$$ \begin{equation} P(Y = y) = \binom{y – 1}{r – 1}\,p^r\,(1-p)^{y-r}, \quad y = r, r+1, r+2, \dots \end{equation} $$
ここで $\binom{y-1}{r-1}$ は「最後（$y$ 回目）が必ず成功であり、それ以前の $y-1$ 回の試行のうち $r-1$ 回が成功である」組み合わせの数です。
次に、なぜ「負の二項分布」という名前なのかを見てみましょう。
名前の由来 — 負の二項係数
「負の二項分布」という名前は、確率質量関数が一般化二項係数の負のパラメータ版と関連していることに由来します。
一般化二項係数は任意の実数 $\alpha$ に対して次のように定義されます。
$$ \binom{\alpha}{k} = \frac{\alpha(\alpha-1)(\alpha-2)\cdots(\alpha-k+1)}{k!} $$
ここで $\alpha = -r$（負の整数）とおくと、
$$ \binom{-r}{x} = \frac{(-r)(-r-1)(-r-2)\cdots(-r-x+1)}{x!} $$
分子から $(-1)^x$ を因子として括り出すと、
$$ \binom{-r}{x} = (-1)^x \frac{r(r+1)(r+2)\cdots(r+x-1)}{x!} = (-1)^x \binom{x+r-1}{x} $$
したがって、
$$ \binom{x+r-1}{x} = (-1)^x \binom{-r}{x} $$
この関係を確率質量関数に代入すると、
$$ P(X = x) = (-1)^x \binom{-r}{x}\,p^r\,(1-p)^x = \binom{-r}{x}\,p^r\,(-1)^x(1-p)^x $$
これは $(p + (1-p))^{-r}$ の二項級数展開の各項に対応しています。
$$ (p + (1-p) \cdot (-1) \cdot (-1))^{-r} \to \sum_{x=0}^{\infty} \binom{-r}{x}\,p^r\,[-(1-p)]^x $$
つまり、負の二項分布の確率質量関数は、$(1-q)^{-r}$（$q = 1-p$）の二項級数展開（ニュートンの一般二項定理）から自然に現れるのです。これが「負の二項」という名前の由来です。
この二項級数の関係から、確率の総和が1であることも確認できます。
$$ \sum_{x=0}^{\infty} \binom{x+r-1}{x}p^r q^x = p^r \sum_{x=0}^{\infty} \binom{-r}{x}(-q)^x = p^r (1-q)^{-r} = p^r \cdot p^{-r} = 1 $$
確率質量関数の正当性が確認できました。次に、モーメント母関数とモーメントを求めましょう。
モーメント母関数と基本統計量
モーメント母関数
$X \sim \text{NB}(r, p)$（失敗回数）のモーメント母関数は次のように計算されます。
$$ M_X(t) = E[e^{tX}] = \sum_{x=0}^{\infty} e^{tx}\binom{x+r-1}{x}p^r(1-p)^x $$
$e^{tx}(1-p)^x = [(1-p)e^t]^x$ とまとめ、先ほどの二項級数を使うと、
$$ M_X(t) = p^r \sum_{x=0}^{\infty} \binom{x+r-1}{x}[(1-p)e^t]^x = p^r [1 – (1-p)e^t]^{-r} $$
ただし $|(1-p)e^t| < 1$、すなわち $t < -\ln(1-p)$ の範囲で収束します。
$$ \begin{equation} M_X(t) = \left(\frac{p}{1 – (1-p)e^t}\right)^r, \quad t < -\ln(1-p) \end{equation} $$
期待値
期待値は $M_X'(0)$ から求められますが、直感的な方法でも計算できます。
$X$ は $r$ 個の独立な幾何分布の和です。すなわち、$X = X_1 + X_2 + \cdots + X_r$ と分解でき、各 $X_i$ は「$i-1$ 回目の成功から $i$ 回目の成功までの失敗回数」であり、$X_i \sim \text{Geom}(p)$（失敗回数版）に従います。幾何分布の期待値は $(1-p)/p$ なので、
$$ \begin{equation} E[X] = r \cdot \frac{1-p}{p} = \frac{r(1-p)}{p} \end{equation} $$
総試行回数 $Y = X + r$ の期待値は $E[Y] = r/p$ です。これは直感に合います。成功確率 $p$ のとき、$r$ 回の成功を得るには平均 $r/p$ 回の試行が必要です。
分散
同様に、幾何分布の分散 $(1-p)/p^2$ を使って、
$$ \begin{equation} \text{Var}(X) = r \cdot \frac{1-p}{p^2} = \frac{r(1-p)}{p^2} \end{equation} $$
ここで重要な点があります。期待値を $\mu = r(1-p)/p$ とおくと、分散は
$$ \text{Var}(X) = \mu + \frac{\mu^2}{r} = \mu\left(1 + \frac{\mu}{r}\right) $$
と書けます。ポアソン分布では $\text{Var} = \mu$ ですが、負の二項分布では $\text{Var} > \mu$ です。この過分散（overdispersion、分散が期待値よりも大きい）の性質が、負の二項分布がポアソン分布の代替として使われる最大の理由です。
$r \to \infty$ のとき $\mu^2/r \to 0$ となり、分散が期待値に近づいてポアソン分布に収束する、ということも確認できます。
この過分散の性質をより深く理解するために、ポアソン・ガンマ混合としての解釈を見てみましょう。
ポアソン・ガンマ混合としての解釈
混合分布としての導出
負の二項分布は、ポアソン分布のパラメータ（強度）がガンマ分布に従って変動する場合の周辺分布としても導出できます。これはベイズ統計の文脈で非常に重要な解釈です。
具体的には、次の階層モデルを考えます。
$$ \begin{align} \Lambda &\sim \text{Gamma}(r,\, \beta) \\ X \mid \Lambda &\sim \text{Poisson}(\Lambda) \end{align} $$
このとき、$X$ の周辺分布（$\Lambda$ を積分消去した分布）を求めます。
$$ P(X = x) = \int_0^{\infty} P(X = x \mid \Lambda = \lambda)\,f_\Lambda(\lambda)\,d\lambda $$
ポアソン分布の確率質量関数とガンマ分布の確率密度関数を代入します。
$$ P(X = x) = \int_0^{\infty} \frac{\lambda^x e^{-\lambda}}{x!} \cdot \frac{\lambda^{r-1}e^{-\lambda/\beta}}{\beta^r\,\Gamma(r)}\,d\lambda $$
指数の部分をまとめると $e^{-\lambda(1 + 1/\beta)} = e^{-\lambda(1+\beta)/\beta}$ となり、$\lambda$ のべきは $\lambda^{x+r-1}$ です。
$$ P(X = x) = \frac{1}{x!\,\beta^r\,\Gamma(r)}\int_0^{\infty} \lambda^{x+r-1}\exp\left(-\frac{(1+\beta)\lambda}{\beta}\right) d\lambda $$
$u = \frac{(1+\beta)\lambda}{\beta}$ と変数変換すると、積分はガンマ関数 $\Gamma(x+r)$ に帰着します。
$$ P(X = x) = \frac{\Gamma(x+r)}{x!\,\Gamma(r)} \cdot \frac{\beta^x}{(1+\beta)^{x+r}} \cdot \frac{1}{\beta^r} \cdot \frac{\beta^{x+r}}{(1+\beta)^{x+r}} \cdot \frac{(1+\beta)^{x+r}}{\beta^x} $$
整理すると（$p = 1/(1+\beta)$ とおくと $1-p = \beta/(1+\beta)$）、
$$ \begin{equation} P(X = x) = \binom{x+r-1}{x}p^r(1-p)^x \end{equation} $$
が得られます。これは負の二項分布の確率質量関数そのものです。
過分散の直感的理解
ポアソン・ガンマ混合の解釈は、過分散の原因を直感的に説明してくれます。
純粋なポアソン分布では、すべての個体（または観測単位）が同じ強度 $\lambda$ でイベントを生成すると仮定します。しかし実際には、個体間で強度がばらつくことが多いです。たとえば、保険の事故件数では、慎重なドライバーと荒い運転をするドライバーでは事故率が異なります。
この個体間のばらつき（異質性）をガンマ分布で表現すると、周辺分布は負の二項分布になります。個体内のばらつき（ポアソン性）に加えて、個体間のばらつき（ガンマ性）が加わることで、分散が期待値を上回る過分散が生じるのです。
$r$ パラメータは個体間のばらつきの程度を制御しており、$r$ が小さいほどばらつきが大きく（過分散が強く）、$r \to \infty$ でばらつきがなくなってポアソン分布に収束します。
この理解は、負の二項分布をカウントデータのモデリングに使う際の最も重要な動機です。
二項分布・ポアソン分布との関係
ポアソン分布への収束
$r \to \infty$, $p \to 1$ として $r(1-p)/p \to \mu$（期待値を一定に保つ）のとき、
$$ \text{NB}(r, p) \xrightarrow{d} \text{Poisson}(\mu) $$
が成り立ちます。これは先ほどのポアソン・ガンマ混合の解釈からも理解できます。$r \to \infty$ でガンマ分布の分散がゼロに近づき、$\Lambda$ がほぼ定数 $\mu$ になるため、ポアソン分布に退化するのです。
二項分布との関係
$Y \sim \text{NB}(r, p)$（$r$ 回目の成功までの総試行回数）のとき、次の関係が成り立ちます。
$$ P(Y > n) = P(\text{Bin}(n, p) < r) $$
つまり、「$r$ 回目の成功が $n$ 回以内に起きない確率」は「$n$ 回の試行で成功が $r$ 回未満である確率」に等しいのです。これは二項分布と負の二項分布が「裏と表の関係」であることの表れです。
幾何分布との関係
$r = 1$ のとき、負の二項分布は幾何分布に一致します。
$$ \text{NB}(1, p) = \text{Geom}(p) $$
逆に言えば、負の二項分布は $r$ 個の独立な幾何分布の和です。この再生性は、モーメント母関数からも容易に確認できます。
$$ M_{\text{NB}(r,p)}(t) = \left(\frac{p}{1-(1-p)e^t}\right)^r = \left[M_{\text{Geom}(p)}(t)\right]^r $$
では、これらの理論をPythonで可視化しましょう。
Pythonによる確率質量関数の可視化
さまざまなパラメータに対する負の二項分布の形状を描画します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) rを変化（pを固定）
ax = axes[0]
p = 0.3
x = np.arange(0, 40)
for r in [1, 2, 5, 10, 20]:
pmf = stats.nbinom.pmf(x, n=r, p=p)
ax.plot(x, pmf, "o-", markersize=4, linewidth=1.5,
label=rf"$r = {r},\, p = {p}$")
ax.set_xlabel("Number of failures x", fontsize=12)
ax.set_ylabel("P(X = x)", fontsize=12)
ax.set_title(rf"Varying $r$ with fixed $p = {p}$", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (b) pを変化（rを固定）
ax = axes[1]
r = 5
for p_val in [0.2, 0.3, 0.5, 0.7, 0.9]:
x_max = max(5, int(r * (1-p_val)/p_val * 3) + 5)
x = np.arange(0, x_max)
pmf = stats.nbinom.pmf(x, n=r, p=p_val)
ax.plot(x, pmf, "o-", markersize=4, linewidth=1.5,
label=rf"$r = {r},\, p = {p_val}$")
ax.set_xlabel("Number of failures x", fontsize=12)
ax.set_ylabel("P(X = x)", fontsize=12)
ax.set_title(rf"Varying $p$ with fixed $r = {r}$", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("negbinom_pmf.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、負の二項分布の形状について以下の特徴が読み取れます。
-
左図: $r$ が大きくなるにつれて分布の中心が右に移動し、形が釣り鐘型に近づく 。$r = 1$（幾何分布）では $x = 0$ が最頻値で右に裾を引く形ですが、$r = 20$ ではピークが明確になり正規分布に類似した形状です。これは中心極限定理の効果で、$r$ 個の独立な幾何分布の和は正規分布に近づきます
-
右図: $p$ が大きい（成功確率が高い）ほど分布の中心が左に寄り、$p$ が小さいほど右に広がる 。$p = 0.9$ では少ない失敗で $r$ 回の成功に達するため、0に近い値に集中しています。$p = 0.2$ では多くの失敗を経験するため、広く分散した分布になります
次に、ポアソン分布との比較を可視化して過分散の特徴を確認しましょう。
ポアソン分布との比較 — 過分散の可視化
同じ期待値を持つポアソン分布と負の二項分布を並べて、分散の違いを視覚的に確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
mu = 10 # 期待値を固定
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) 確率質量関数の比較
ax = axes[0]
x = np.arange(0, 40)
# ポアソン分布
pmf_poisson = stats.poisson.pmf(x, mu=mu)
ax.plot(x, pmf_poisson, "bs-", markersize=5, linewidth=1.5, label=rf"Poisson($\mu={mu}$)")
# 負の二項分布（さまざまなr）
for r in [1, 3, 10]:
p_nb = r / (r + mu) # E[X] = r(1-p)/p = mu から p = r/(r+mu)
pmf_nb = stats.nbinom.pmf(x, n=r, p=p_nb)
variance = mu * (1 + mu/r)
ax.plot(x, pmf_nb, "o-", markersize=4, linewidth=1.5,
label=rf"NB($r={r}$), Var={variance:.1f}")
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("P(X = x)", fontsize=12)
ax.set_title(rf"Poisson vs Negative Binomial ($\mu = {mu}$)", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (b) シミュレーションによる分散の比較
ax = axes[1]
np.random.seed(42)
n_samples = 10000
samples_poisson = np.random.poisson(lam=mu, size=n_samples)
r_val = 3
p_nb = r_val / (r_val + mu)
samples_nb = np.random.negative_binomial(n=r_val, p=p_nb, size=n_samples)
ax.hist(samples_poisson, bins=np.arange(-0.5, 50.5, 1), density=True,
alpha=0.5, color="blue", label=f"Poisson (var={samples_poisson.var():.1f})")
ax.hist(samples_nb, bins=np.arange(-0.5, 50.5, 1), density=True,
alpha=0.5, color="red", label=f"NB (var={samples_nb.var():.1f})")
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("Density", fontsize=12)
ax.set_title(rf"Simulated samples ($\mu = {mu}$, NB: $r = {r_val}$)", fontsize=13)
ax.legend(fontsize=10)
ax.set_xlim(-1, 50)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("negbinom_vs_poisson.png", dpi=150, bbox_inches="tight")
plt.show()
この比較から、以下の重要な知見が得られます。
-
左図: 同じ期待値でも、$r$ が小さいほど分布が広がる 。$r = 1$（幾何分布の一般化）では裾が非常に重く、30以上の値も無視できない確率で出現します。一方、$r = 10$ ではポアソン分布にかなり近い形をしています。$r \to \infty$ でポアソン分布に収束するという理論が視覚的に確認できます
-
右図: ヒストグラムで見ると、負の二項分布（赤）はポアソン分布（青）よりも広く散らばっている 。分散の値を比較すると、ポアソン分布の分散はほぼ10（= $\mu$）であるのに対し、負の二項分布の分散はそれを大きく上回っています。この過分散は、実データでポアソン分布が当てはまらない場合の典型的なパターンです
次に、実際のデータへの応用として最尤推定を実装してみましょう。
過分散データへの応用 — 最尤推定
実際のカウントデータにポアソン分布と負の二項分布を当てはめ、どちらがより適切かを比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.optimize import minimize_scalar, minimize
np.random.seed(42)
# 過分散のあるカウントデータを生成
# （実際にはRNA-seq等から得られるようなデータをシミュレーション）
true_r = 3
true_mu = 8
true_p = true_r / (true_r + true_mu)
data = np.random.negative_binomial(n=true_r, p=true_p, size=200)
# --- ポアソン分布の最尤推定 ---
lambda_mle = np.mean(data)
# --- 負の二項分布の最尤推定 ---
def neg_loglik_nb(params, data):
r, p = params
if r <= 0 or p <= 0 or p >= 1:
return 1e10
return -np.sum(stats.nbinom.logpmf(data, n=r, p=p))
result = minimize(neg_loglik_nb, x0=[2.0, 0.3], args=(data,),
method="Nelder-Mead")
r_mle, p_mle = result.x
print(f"データの統計量: 平均 = {data.mean():.2f}, 分散 = {data.var():.2f}")
print(f"分散/平均比 = {data.var()/data.mean():.2f} (1なら等分散)")
print(f"ポアソンMLE: λ = {lambda_mle:.2f}")
print(f"負の二項MLE: r = {r_mle:.2f}, p = {p_mle:.3f}")
# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# ヒストグラムと当てはめ
ax = axes[0]
x_max = max(data) + 2
x = np.arange(0, x_max)
ax.hist(data, bins=np.arange(-0.5, x_max + 0.5, 1), density=True,
alpha=0.6, color="gray", edgecolor="white", label="Observed data")
pmf_poisson = stats.poisson.pmf(x, mu=lambda_mle)
ax.plot(x, pmf_poisson, "bs-", markersize=5, linewidth=1.5,
label=rf"Poisson($\lambda={lambda_mle:.1f}$)")
pmf_nb = stats.nbinom.pmf(x, n=r_mle, p=p_mle)
ax.plot(x, pmf_nb, "ro-", markersize=5, linewidth=1.5,
label=rf"NB($r={r_mle:.1f}$, $p={p_mle:.2f}$)")
ax.set_xlabel("Count", fontsize=12)
ax.set_ylabel("Density", fontsize=12)
ax.set_title("Model comparison: Poisson vs NB", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# QQプロット（ポアソン vs データ）
ax = axes[1]
sorted_data = np.sort(data)
n = len(data)
theoretical_quantiles_pois = stats.poisson.ppf(
(np.arange(1, n+1) - 0.5) / n, mu=lambda_mle)
theoretical_quantiles_nb = stats.nbinom.ppf(
(np.arange(1, n+1) - 0.5) / n, n=r_mle, p=p_mle)
ax.scatter(theoretical_quantiles_pois, sorted_data, alpha=0.5, s=20,
color="blue", label="Poisson")
ax.scatter(theoretical_quantiles_nb, sorted_data, alpha=0.5, s=20,
color="red", label="NB")
max_val = max(max(sorted_data), max(theoretical_quantiles_pois)) + 2
ax.plot([0, max_val], [0, max_val], "k--", linewidth=1, label="y = x")
ax.set_xlabel("Theoretical quantiles", fontsize=12)
ax.set_ylabel("Observed quantiles", fontsize=12)
ax.set_title("Q-Q plot", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_aspect("equal")
plt.tight_layout()
plt.savefig("negbinom_fitting.png", dpi=150, bbox_inches="tight")
plt.show()
# AICの比較
loglik_pois = np.sum(stats.poisson.logpmf(data, mu=lambda_mle))
loglik_nb = -neg_loglik_nb([r_mle, p_mle], data)
aic_pois = -2 * loglik_pois + 2 * 1 # ポアソンはパラメータ1個
aic_nb = -2 * loglik_nb + 2 * 2 # 負の二項はパラメータ2個
print(f"\nAIC比較:")
print(f"ポアソン: AIC = {aic_pois:.1f}")
print(f"負の二項: AIC = {aic_nb:.1f}")
print(f"ΔAIC = {aic_pois - aic_nb:.1f} (正なら負の二項が優位)")
このフィッティング結果から、以下のことが読み取れます。
-
左図: 負の二項分布（赤）の方がポアソン分布（青）よりもデータのヒストグラムに合っている 。ポアソン分布はデータの裾の重さを表現できず、大きな値を過小評価し、ピーク付近を過大評価しています。負の二項分布は裾まで含めてデータの形状をよく捉えています
-
右図のQQプロット: 負の二項分布（赤）の点は対角線（$y = x$）に近い のに対し、ポアソン分布（青）の点は対角線から系統的にずれています。特に上側（大きな値）でポアソン分布のずれが顕著であり、データの裾がポアソン分布の予測よりも重いことを示しています
-
AICの比較: 負の二項分布のAICがポアソン分布のAICよりも小さく、パラメータ数の増加（1個 → 2個）を考慮してもなお負の二項分布が優位であることが確認できます
連続パラメータへの拡張
ここまで $r$ は正の整数として扱ってきましたが、ポアソン・ガンマ混合の解釈を使えば、$r$ を任意の正の実数に拡張できます。この場合、$\binom{x+r-1}{x}$ は一般化二項係数 $\frac{\Gamma(x+r)}{x!\,\Gamma(r)}$ と解釈します。
$$ P(X = x) = \frac{\Gamma(x+r)}{x!\,\Gamma(r)}\,p^r\,(1-p)^x, \quad x = 0, 1, 2, \dots $$
この連続パラメータ版は、統計ソフトウェア（R、Python の scipy 等）で実装されているのが通常です。$r$ が整数でなくても「$r$ 回目の成功まで」という解釈は失われますが、ポアソン・ガンマ混合としての解釈は依然として有効です。$r$ は過分散の程度を制御するパラメータと理解できます。
まとめ
本記事では、負の二項分布の定義から導出、性質、応用までを解説しました。
- 負の二項分布は「$r$ 回目の成功までの失敗回数」の分布 であり、幾何分布の一般化です。成功確率 $p$ と成功回数 $r$ の2つのパラメータを持ちます
- 確率質量関数 $P(X = x) = \binom{x+r-1}{x}p^r(1-p)^x$ は、二項係数の負のパラメータ版（ニュートンの一般二項定理）と関連しています
- 基本統計量: $E[X] = r(1-p)/p$, $\text{Var}(X) = r(1-p)/p^2$。分散は常に期待値よりも大きく（$\text{Var} = \mu + \mu^2/r$）、過分散を示します
- ポアソン・ガンマ混合として導出でき、個体間の異質性（ばらつき）による過分散を自然に説明します
- $r \to \infty$ でポアソン分布に収束し、$r = 1$ で幾何分布に一致します
- 過分散のあるカウントデータ（RNA-seq、保険、生態学など）のモデリングに広く応用されます
次のステップとして、以下の記事も参考にしてください。