可視化で理解する「負の二項分布」

みどりぼんでカウントデータの過分散対策のために使われると書かれている負の二項分布ですが、Wikipediaの説明を読んでもよく分かりません。

そこでおススメなのが、このスライドです。

負の二項分布について from Hiroshi Shimizu

ようするに、負の二項分布は、 $\lambda$ がガンマ分布に従うようなポアソン分布だと思えばだいたい OK みたいです。

今日はこれを可視化してみます。

負の二項分布(Negative Binomial Distribution)

負の二項分布はパラメータを 2つ持ちます。成功回数を表す $r$ と成功確率を表す $p$ です。

統計言語 R には負の二項分布に従う乱数を生成する関数 rnbinom() があり、これらのパラメータはそれぞれ引数 size と prob に対応しています。

したがって、 $r=4$ , $p=0.3$ の負の二項分布は次のようにして描画できます。

negative_binom_values <- rnbinom(10000, size = 4, prob = 0.3)
barplot(table(negative_binom_values)[1:21])

これに対応するのは、 $\lambda$ がガンマ分布 ${\rm Gamma}(r, p/(1-p))$ に従うポアソン分布です。
このガンマ分布から $\lambda$ を生成し、各 $\lambda$ に従うポアソン分布から乱数を生成という手順でやってみましょう。

gamma_values <- rgamma(10000, 4, 0.3/(1-0.3))
poisson_values <- sapply(gamma_values, function(lambda) rpois(1, lambda))
barplot(table(poisson_values)[1:21], main="λがガンマ分布に従うポアソン分布")

見事に同じような分布の形状になりました。

少しわかりにくいので、この様子をアニメーションにしてみましょう。*1

最下段がガンマ分布であり、このガンマ分布からサンプリングされた値を赤い線で示しています。
この値をパラメータ $\lambda$ として持つポアソン分布が中段のグラフです。
このポアソン分布からサンプリングされた値を赤い線で示し、この値からなるヒストグラムを最上段に描いています。
毎回異なるポアソン分布から値がサンプリングされており、それらの値が負の二項分布を形成していく様子がわかります。

混合ポアソン分布

同じ状況を別の角度から見てみましょう。
ガンマ分布からサンプリングされた値をパラメータ $\lambda$ として持つポアソン分布を複数作成します。

library(ggplot2)

df <- Reduce(rbind, Map(function(lambda) data.frame(lambda, x=0:20, y=dpois(0:20, lambda)), gamma_values[1:16]))
ggplot(df, aes(x=x, y=y, fill=lambda)) + geom_bar(stat="identity") + facet_wrap(~ lambda, nrow=4)

これらの分布を合成したものが負の二項分布です。

ggplot(df, aes(x=x, y=y, fill=lambda)) + geom_bar(stat="identity") + theme(legend.position = "none")

上図は 16個のポアソン分布を合成したものですが、1000個ほど合成すればほとんど負の二項分布と一致することがわかります。

library(gridExtra)

df2 <- Reduce(rbind, Map(function(lambda) data.frame(lambda, x=0:20, y=dpois(0:20, lambda)), gamma_values[1:1000]))
g1 <- ggplot(df2, aes(x=x, y=y, fill=lambda)) + geom_bar(stat="identity") + theme(legend.position = "none") + xlab("混合ポアソン分布(n=1000)") + ylab("")
g2 <- qplot(x=0:20, y=dnbinom(0:20, 4, 0.3)) + geom_bar(stat="identity") + xlab("負の二項分布") + ylab("")
grid.arrange(g1, g2, ncol=1)