負の二項分布(Negative Binomial Distribution)について学ぶ
負の二項分布とは
負の二項分布を使う機会があったのだけど、定義がいっぱいあってよくわからないので調べたことをまとめてみる。
日本語版wikipediaだと一つの定義しかないが、英語版wikipediaの項目(Negative binomial distribution)が割と充実していたので英語版wikipediaを見ながら学ぶ。
NBDと呼ばれることもあるね。
二項分布については使われる機会も多く、少し調べてみると充実した解説が山ほど出てくるが、負の二項分布についてはあまり充実した情報はないようだ。英語版wikiによると全部で5つの定義があるみたい。式書くのめんどくさいので以下スクショ。
↓リンク
Negative binomial distribution - Wikipedia
ここでは試行回数:n、成功した回数:k、失敗した回数:r、成功する確率:pと定義する。
当たり前だが k+r=n の関係が成り立つ。
5つはそれぞれ
1:r回失敗する前に何回成功したか
2:r回失敗するのに必要な試行回数
3:k回成功する前に失敗した数
4;k回成功するのに必要な試行回数
5:n回の試行で何回成功したか(二項分布)
となっている。
一番わかりやすいのは5番目の定義で、これは二項分布の定義そのもの。
n回の試行に対して何回成功したかという確率分布である。
コインを100回投げたとき何回表が出たか、みたいなやつ。
成功した回数の分布を見ているのでf(k;,n,p)とkの関数になっていることがわかる。
4つ目の定義については確か日本語版wikipediaで説明されていた気がする。
このように、二項分布の場合は定義が一つだが、負の二項分布の場合はそれぞれ定義が微妙に違うあるので、もし「負の二項分布を使って〜」みたいな説明を見たときはどの定義なのかちゃんと注意する必要がある。
これら5つの式は割と簡単に導くことができる。試しに1つ目の定義について考えてみる。
n-1 (k+r-1)回の試行でk回成功して、n回目で失敗すればいいので、
\begin{align}f(k:r,p) &=\begin{pmatrix}{r+k{-}1}\\{k}\end{pmatrix}p^{k}(1-p)^{r-1}(1-p)\\ &=\begin{pmatrix}{r+k{-}1}\\{k}\end{pmatrix}p^{k}(1-p)^{r}\end{align}
てな感じ。他も同様に導ける。そんなに難しくないね。
せっかくなのでもう少し踏み込んでみよう。
実数に拡張する
この定義式だとrもkも整数しかとれない。成功回数が何回か、みたいな議論をしてるので当然だけども。これを正の実数に拡張してみよう。そうすれば連続的な関数として取り扱うことができるようになるはずだ。全部の定義でやるのはかったるいので、1番目の定義に絞ることとする。
まずは階乗の部分をガンマ関数に
ガンマ関数についてはぐぐればいくらでも出てくるので詳しくは書かないけど、これを使えば階乗の概念を連続的な量に拡張できる。
\begin{align}\Gamma{(n+1)=n!}\end{align}
であることを考えれば
\begin{align}\begin{pmatrix}{r+k{-}1}\\{k}\end{pmatrix}&=\frac{(r+k{-}1)!}{k!(r-1)!}\\&=\frac{\Gamma{(r+k)}}{\Gamma{(k+1)}\Gamma{(r)}}\end{align}
となるね。
分布の平均値をパラメータとして使う
さて、これで一通り準備は揃ったような気もするけど、最後に分布の平均値(期待値)を計算してそいつをパラメータにしよう。
なんでそんなことするのか。
ポアソン分布にしろ、ガウス分布にしろ、平均値を表すパラメータが式の中に入ってる。なので成功確率pをそのまま使って表すよりも平均値を用いたほうが分布の形状を考えるとき扱いやすいよね。まぁ理由なんてなんでもいいんだけど、とにかく平均値を考える。
と言っても考え方は二項分布のときと同じで、試行回数にpをかければ成功回数になるし、(1-p)をかければ失敗回数になるので以下の式が成り立つはずだ。
\begin{align}\frac{\langle{k}\rangle}{p}=\frac{r}{1-p}=k+r\end{align}
\begin{align}{\langle{k}\rangle}=\frac{pr}{1-p}\equiv{m}\end{align}
平均値を表すパラメータをmと定義した。
あれ、kには期待値記号が付くのにrに付かないのはなぜ?と思う人もいるかもしれないが、現在の定義だとkについての関数になっていてrは固定値なので期待値記号は付かない。もちろん別の定義のNBDの定義だと変数も変わる。
これらを用いれば最終的には
\begin{align}f(k;r,m)&=\frac{\Gamma{(r+k)}}{\Gamma{(k+1)}\Gamma{(r)}}\left(\frac{m}{r+m}\right)^{k}\left(\frac{r}{r+m}\right)^{r} \end{align}
あるいは
\begin{align}f(k;r,m)&=\frac{\Gamma{(r+k)}}{\Gamma{(k+1)}\Gamma{(r)}}\left(\frac{m}{r}\right)^{k}\left(\frac{m}{r}+1\right)^{-(r+k)} \end{align}
とかくことができる。
上の方の式は英語版wikipediaにも載っているので参照にされたし。
そもそもなぜ”負”の二項分布なのか?
これでなんとなく負の二項分布についての定義はわかったと思うが、なぜこのような名前がついているのだろうか。実は負の二項分布には負の二項定理なるものを使えるのだ。
通常の二項定理は
\begin{align}(1+p)^{n}=\begin{pmatrix}{n}\\{k}\end{pmatrix}\sum_{k=0}^n{p^{k}}\end{align}
と書くことができる。通常nは自然数なでなければならないが、これを負に拡張して
\begin{align}(1+p)^{-n}=\begin{pmatrix}{-n}\\{k}\end{pmatrix}\sum_{k=0}^n{p^{k}}\end{align}
なるものを考える。ここでやっかいになるのが \begin{align}\begin{pmatrix}{-n}\\{k}\end{pmatrix}\end{align}であるが、これを変形すると
\begin{align}\begin{pmatrix}{-n}\\{k}\end{pmatrix}&=\frac{(-n)(-n-1)...(-n-k+1)}{k!}\\&=\frac{(-1)^k(n+k-1)(n+k-2)...(n)}{k!}\\&=(-1)^{k}\begin{pmatrix}{n+k-1}\\{n-1}\end{pmatrix}\end{align}
と変形することができる。これを負の二項分布の式に適用すると、
\begin{align}\begin{pmatrix}{r+k{-}1}\\{k}\end{pmatrix}p^{k}(1-p)^{r}=\begin{pmatrix}{-r}\\{k}\end{pmatrix}p^{k}(1-p)^{r}\end{align}
と、すっきりする形に変形することができた。この変形によってどれだけすっきりしたかイメージできるだろうか。例えば、通常の二項分布における\begin{align}\begin{pmatrix}{n}\\{k}\end{pmatrix}\end{align}の部分は、nが固定の値でkが変数なので、例えば和を取る場合など二項定理の適用がしやすい。一方で負の二項分布の最初の定理の場合だと、\begin{align}\begin{pmatrix}{k+r-1}\\{r-1}\end{pmatrix}\end{align}の部分において、rが固定、kが変数になっているので、二項定理と比べて定数、変数の場所が逆になってしまっている。変形後は\begin{align}\begin{pmatrix}{-r}\\{k}\end{pmatrix}\end{align}となっているので\begin{align}\begin{pmatrix}{定数}\\{変数}\end{pmatrix}\end{align}の形に揃えることができた。これはこのあと出てくる母関数の計算においても用いる。
母関数を計算して平均と分散を計算してみる
さて、負の二項分布の平均と分散を計算してみよう。平均値についてはさっき求めたけど、母関数を計算すればn回微分でn次のモーメントorキュムラントが求まるので便利。1次のモーメント、キュムラントは平均値にあたり、2次のモーメント、キュムラントが分散にあたる。
最後に二項分布と負の二項分布、何が違うのか比べてみる。
二項分布の場合
二項分布の式は上にあるように
\begin{align}f(k;n,p)=\begin{pmatrix}{n}\\{k}\end{pmatrix}p^{k}(1-p)^{k-1}\end{align}
である。
モーメント母関数M(θ)の定義は
\begin{align}\sum_{k}e^{\theta{k}}f\end{align}
一応確認しておくと、fは確立密度関数である。fの部分に二項分布の式をぶち込むと
\begin{align}M(\theta)=\sum_{k}\begin{pmatrix}{n}\\{k}\end{pmatrix}(p{e^\theta})^{k}(1-p)^{n-k}=(pe^{\theta}+1-p)^n\end{align}
等式の変形で二項定理を用いた。
これをn回微分すればn次のモーメントが求まる。モーメントから計算してもいいんだけど、今回はキュムラントからアプローチしてみる。キュムラント母関数の定義はK(θ)=lnM(θ)なので、モーメント母関数のlogをとるだけ。モーメント同様キュムラント母関数をn回微分すると、n次キュムラントが出てくる。n次キュムラントをCnとすると、
\begin{align}C_{n}=\frac{\partial}{\partial\theta}K(\theta)|_{\theta=0}\end{align}
となる。
平均と分散を調べたいのであれば、1階微分と2階微分を調べればいい。M(θ)のlogをとると肩のnが落っこちてくるだけなので
\begin{align}K(\theta)=n(pe^\theta+1-p)\end{align}
となる。これをθで微分していく。細かい計算は省略するけど、計算すると
\begin{align}\mu=np\end{align}
\begin{align}\sigma^2=np(1-p)\end{align}
ここでμ、σ^2はそれぞれ平均、 分散を表し、C1, C2にあたる。また、
\begin{align}\sigma^2=\mu(1-p)\end{align}
の関係が成り立つ。1-pをεと定義しなおせば
\begin{align}\sigma^2=\mu\epsilon\end{align}
というすっきりした関係が導ける。ここでpは確率を表すので、ε<1となる。つまり、
二項分布においては平均は分散より大きいということがわかる。
負の二項分布の場合
負の二項分布でも同じことをやってみよう。式は変形後の
\begin{align}f(k;r,p)=\begin{pmatrix}{-r}\\{k}\end{pmatrix}p^{k}(1-p)^{r}\end{align}
を用いる。
モーメント母関数は
\begin{align}M(\theta)=\sum_{k}\begin{pmatrix}{-r}\\{k}\end{pmatrix}(p{e^\theta})^{k}(1-p)^{r}=(\frac{1-p}{1-pe^{\theta}})^r\end{align}
等式の変形で負の二項定理を用いた。
同様にM(θ)のlogをとると肩のrが落っこちてくるだけなので
\begin{align}K(\theta)=r(\frac{1-p}{1-pe^{\theta}})\end{align}
となる。同様の計算すると
\begin{align}\mu=\frac{pr}{1-p}\end{align}
\begin{align}\sigma^2=\frac{pr}{(1-p)^2}\end{align}
ここで
\begin{align}\sigma^2=\frac{\mu}{(1-p)}\end{align}
の関係が成り立つ。二項分布の場合は平均に1-pを掛けると分散になるのに対して、負の二項分布の場合は割り算になっているという非常に綺麗な関係を導くことができた。さらに1/1-pをεと定義しなおせば
\begin{align}\sigma^2=\mu\epsilon\end{align}
と、二項分布と全く同じ形で表すこともできる。
このとき、負の二項分布においてはε>1となる。つまり、
負の二項分布においては平均は分散より小さいということがわかる。
途中まで全然違う計算をしているようで、最後に綺麗に同じような形で表すことができるのは驚きだ。ちなみに、ポアソン分布においては平均と分散は同じなので、無理やりまとめるとすれば
二項分布:ε>1 (平均>分散)
ポアソン分布:ε=1 (平均=分散)
負の二項分布:ε<1 (平均<分散)
となるね。よくできてるなぁ。
ちなみに、モーメント母関数から計算しても、別の定義の負の二項分布から計算しても同様の結果を導くことはできる。さらに言えば、平均と分散を計算しただけなら別に母関数を使わなくてもできる。ただ、俺はキュムラント母関数が便利で応用がきくので好きだ。
追記:3次以上の場合
上では平均と分散について、εを定義すれば二項分布と負の二項分布がまったく同じ関係式で表すことができるを示した。
実はこれ、3次以降のすべての次数のキュムラントにおいても一致する。
文明の利器(mathematica)を使って、キュムラント母関数を微分して6次のキュムラントまで計算してみると、どちらもまったく同じ結果になった。すっきり。
\begin{align}&C_{1}=\mu\\&C_{2}=\mu\epsilon\\&C_{3}=\mu\epsilon(2\epsilon-1)\\&C_{4}=\mu\epsilon(6\epsilon^{2}-6\epsilon+1)\\&C_{5}=\mu\epsilon(2\epsilon-1)(12\epsilon^{2}-12\epsilon+1)\\&C_{6}=\mu\epsilon(120\epsilon^{4}-240\epsilon^{3}+150\epsilon^{2}-30\epsilon+1)\end{align}
おわり。