負の二項分布

負の二項分布は分散の大きいカウントデータによく用いられる統計モデルである。パスカル分布やポアソンガンマ分布などとも呼ばれたりする。負の二項分布は、次のように捉えることができる。

  1. 試行が k 回成功するまでに必要な試行回数 x の分布
  2. 試行が k 回成功するまでに失敗した回数 y の分布

試行が k 回成功するまでの必要な試行回数 x の分布

確率 p で成功する試行があり、x 回目の試行で、 k 回目の成功となったときの x の確率分布は次のようになる。

\[ \begin{eqnarray} P(X=x) &=& \begin{pmatrix} x-1 \\ k-1 \end{pmatrix} p^{k}(1-p)^{x-k} \\ &=& \frac{\Gamma (x)}{\Gamma (k) \Gamma(x-k+1)} p^{k}(1-p)^{x-k} \end{eqnarray} \]

このとき、期待値と分散を r および p で表すと、次のようになる。

\[ E(X) = \frac{r}{p} \] \[ Var(X) = \frac{r(1-p)}{p^{2}} \]

試行が k 回成功するまでに失敗した回数 y の分布

確率 p で成功する試行があり、k 回目の成功に達するまでの失敗した回数 y の確率分布が次のようになる。

\[ \begin{eqnarray} P(Y=y) &=& \begin{pmatrix}k+y-1 \\ k-1 \end{pmatrix} p^{k}(1-p)^{y} \\ &=& \frac{\Gamma (k+y)}{\Gamma (k) \Gamma (y+1)}p^{k}(1-p)^{y} \end{eqnarray} \]

このとき、期待値と分散は次のようになる。

\[ E(Y) = \frac{k(k+1)(1-p)^{2}}{p^{2}} \] \[ V(Y) = \frac{k(1-p)}{p^{2}} \]

ポアソン分布と負の二項分布

負の二項分布の確率分布関数はポアソン分布から導ける。 確率変数 Y がポアソン分布を従うとき、Y の確率分布関数は次のように書ける。

\[ P(Y=y) = \frac{\lambda ^{y}e^{-\lambda}}{y!} \]

ポアソン分布のパラメーターは λ である。 そこで、 λ がガンマ分布 Gamma(φ-1, μφ) に従うものとすると、上記の確率分布関数は、次のようになる。

\[ \begin{eqnarray} P(Y=y|\mu, \phi) &=& \int_{0}^{\infty}P(Y|\lambda)P(\lambda |\mu,\phi)d\lambda \\ &=& \int_{0}^{\infty} \frac{\lambda^{y}e^{-\lambda}}{y!} \frac{\lambda^{\phi^{-1}-1}e^{-\lambda\mu^{-1}\phi^{-1}}}{(\mu\phi)^{\phi^{-1}}\Gamma (\phi^{-1})}d\lambda \\ &=& \frac{\Gamma(y + \phi^{-1})}{\Gamma(\phi^{-1})\Gamma(y+1)}\left(\frac{1}{1+\mu\phi}\right)^{\phi^{-1}}\left(\frac{\mu}{\phi^{-1} + \mu}\right)^{y} \end{eqnarray} \]

これが負の二項分布の確率分布関数である。生物学では次世代シーケンサーなどから得られるリードカウントデータの分布などが、この形で表される。

負の二規分布の関数

R の負の二項分布に関連する関数として、以下のようなものがある。試行回数は size、成功する確率は prob で指定する。成功する確率 prob の代わりに、直接分布の平均 mu を指定することもできる。

関数パラメーター
乱数関数rnbinom(n, size, prob, mu)
確率関数dnbinom(x, size, prob, mu, log=FALSE)
分布関数pnbinom(x, size, prob, mu, lower.tail=TRUE, log.x=FALSE)
クォンタイル関数qnbinom(x, size, prob, mu, lower.tail=TRUE, log.x=FALSE)

平均μ、分散σ2の負の二項分布に従う乱数を、1000000 個生成する例。

m <- 10             # 平均
v <- 12             # 分散
s <- m^2 / (v - m)  # 分散が12となるように負の二項分布のsizeパラメーターを設定

# 1000000個の乱数生成
x <- rnbinom(1000000, size = s, mu = m)


# 乱数の平均を求める
mean(x)
## [1] 10.00107

# 乱数の分散を求める
var(x)
## [1] 12.00969