ほくそ笑む

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

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

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

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

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

負の二項分布(Negative Binomial Distribution)

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

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

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

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

これに対応するのは、 がガンマ分布 に従うポアソン分布です。

このガンマ分布から を生成し、各 に従うポアソン分布から乱数を生成という手順でやってみましょう。

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

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

この値をパラメータ として持つポアソン分布が中段のグラフです。

このポアソン分布からサンプリングされた値を赤い線で示し、この値からなるヒストグラムを最上段に描いています。

毎回異なるポアソン分布から値がサンプリングされており、それらの値が負の二項分布を形成していく様子がわかります。

混合ポアソン分布

同じ状況を別の角度から見てみましょう。

ガンマ分布からサンプリングされた値をパラメータ として持つポアソン分布を複数作成します。

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)

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

まとめ

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

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

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

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

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

• 確率分布いろいろ

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

補足(アニメーション作成コード)

library(animation)
library(ggplot2)
library(gridExtra)

saveGIF({
  g_nbiom <- qplot(x=NULL) + geom_histogram(binwidth=1) + xlim(0,20) + ylim(0, 1) + theme(axis.text.y=element_blank()) + xlab("負の二項分布") + ylab("")
  for(i in 1:16) {
    g_pois <- qplot(x=0:20, y=dpois(0:20, gamma_values[i])) + geom_bar(stat="identity") + theme(axis.text.y=element_blank()) + xlab("ポアソン分布") + ylab("")
    g_gamma <- qplot(x=0:20, y=dgamma(0:20, 4, 0.3/(1-0.3))) + geom_bar(stat="identity") + geom_vline(aes(xintercept = gamma_values[i], color="red"), size=2) + theme(axis.text.y=element_blank()) + xlab("ガンマ分布") + ylab("")
    grid.arrange(g_nbiom, g_pois, g_gamma, ncol=1, nrow=3)
    if(poisson_values[i] < 20) {
      g_pois <- qplot(x=0:20, y=dpois(0:20, gamma_values[i])) + geom_bar(stat="identity") + geom_vline(aes(xintercept = poisson_values[i], color="red"), size=2) + theme(axis.text.y=element_blank()) + xlab("ポアソン分布") + ylab("")
    }
    j <- i
    g_nbiom <- qplot(x=poisson_values[1:j]) + geom_histogram(binwidth=1) + xlim(0,20) + theme(axis.text.y=element_blank()) + xlab("負の二項分布") + ylab("")
    grid.arrange(g_nbiom, g_pois, g_gamma, ncol=1, nrow=3)
  }
  for(i in 17:60) {
    if(poisson_values[i] < 20) {
      g_pois <- qplot(x=0:20, y=dpois(0:20, gamma_values[i])) + geom_bar(stat="identity") + geom_vline(aes(xintercept = poisson_values[i], color="red"), size=2) + theme(axis.text.y=element_blank()) + xlab("ポアソン分布") + ylab("")
      g_nbiom <- qplot(x=poisson_values[1:i]) + geom_histogram(binwidth=1) + xlim(0,20) + theme(axis.text.y=element_blank()) + xlab("負の二項分布") + ylab("")
      g_gamma <- qplot(x=0:20, y=dgamma(0:20, 4, 0.3/(1-0.3))) + geom_bar(stat="identity") + geom_vline(aes(xintercept = gamma_values[i], color="red"), size=2) + theme(axis.text.y=element_blank()) + xlab("ガンマ分布") + ylab("")
    }
    grid.arrange(g_nbiom, g_pois, g_gamma, ncol=1, nrow=3)
  }
  for(i in seq(100, 4000, by=80)) {
    g_pois <- qplot(x=0:20, y=dpois(0:20, gamma_values[i])) + geom_bar(stat="identity") + geom_vline(aes(xintercept = poisson_values[i], color="red"), size=2) + theme(axis.text.y=element_blank()) + xlab("ポアソン分布") + ylab("")
    g_nbiom <- qplot(x=poisson_values[1:i]) + geom_histogram(binwidth=1) + xlim(0,20) + theme(axis.text.y=element_blank()) + xlab("負の二項分布") + ylab("")
    g_gamma <- qplot(x=0:20, y=dgamma(0:20, 4, 0.3/(1-0.3))) + geom_bar(stat="identity") + geom_vline(aes(xintercept = gamma_values[i], color="red"), size=2) + theme(axis.text.y=element_blank()) + xlab("ガンマ分布") + ylab("")
    grid.arrange(g_nbiom, g_pois, g_gamma, ncol=1, nrow=3)
  }
}, interval=rep(c(0.5, 0.2, 0.1, 3), c(32, 44, 48, 1)), movie.name="nbiom_animation3.gif")

library(animation)
library(ggplot2)

saveGIF({
  for(i in 1:100) {
    df <- Reduce(rbind, Map(function(lambda) data.frame(lambda, x=0:20, y=dpois(0:20, lambda)), gamma_values[1:i]))
    g <- ggplot(df, aes(x=x, y=y, fill=lambda)) + geom_bar(stat="identity") + theme(legend.position = "none", axis.text.y=element_blank()) + ggtitle(sprintf("n = %d", i))
    print(g)
  }
}, interval= c((1-(1:99/100)^2)/2, 3), movie.name="nbiom_animation2.gif")

ツイートする

Permalink | コメント(0)