Hatena::ブログ(Diary)

ほくそ笑む このページをアンテナに追加 RSSフィード

2015-10-12

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

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

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

ようするに、負の二項分布は、¥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])

f:id:hoxo_m:20151012022958p:image

これに対応するのは、¥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="λがガンマ分布に従うポアソン分布")

f:id:hoxo_m:20151012022959p:image

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

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

f:id:hoxo_m:20151017134113g:image

最下段がガンマ分布であり、このガンマ分布からサンプリングされた値を赤い線で示しています。

この値をパラメータ  ¥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)

f:id:hoxo_m:20151012034835p:image

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

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

f:id:hoxo_m:20151012034836p:image

上図は 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)

f:id:hoxo_m:20151012041307p:image

この様子も無駄にアニメーションにしてみましょう。

f:id:hoxo_m:20151012042131g:image

まとめ

みどりぼんでは、個体差を持つカウントデータに対処するために GLMM(一般化線形混合モデル)を導入しています。

この本では、GLMM の例として、R の glmmML パッケージを使って「パラメータが正規分布に従う二項分布」の推定を行っています。

これに対して、負の二項分布は「パラメータがガンマ分布に従うポアソン分布」だと考えれば同じ枠組みでとらえることができると思います。

R で負の二項分布を用いた GLMM は、MASS パッケージの glm.nb() 関数によってできます。

負の二項分布がパラメータによってどのような形状をとるかの感覚を知りたい場合は、次のサイトが便利です。

みなさまの理解の一助になれば幸いです。

追記

ホクソエムの切り込み隊長 teramonagi さんが glm.nb() の使い方について詳しく書いてます。

*1:アニメーションの作成コード:https://gist.github.com/hoxo-m/3060e2fe269ff1fdf4a3

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証