覚書：負の二項分布のパラメータ

 　確率分布のパラメータの取り方は様々で，ソフトや関数によって定義が異なってたりする．明らかな変な値が出力されるので，「あ，定義が違うんだな」ということはすぐにわかるのだが，そのたびに変換方法を調べなおすのは結構面倒である．とくに私は，メモを捨てる・失くす・読み返しても意味不明ということがよくあり，そのたびに調べなおさないといけない羽目になる．

　この問題の典型は，負の二項分布である．普段はあんまり使わないRのglm.nbという関数を使う機会があって，ちょっとイライラしたので「なくさない覚書」として，整理した内容をまとめておく．

1. ポアソンーガンマ混合分布の場合

　統計モデリングという文脈では， ポアソンーガンマ混合分布として，負の二項分布が登場する機会が多い．すなわち，離散確率変数yiが，パラメータλiをもつポアソン分布に従っており，そのパラメータλiがガンマ分布に従っているという状況である．

yi ~ Poisson (λi)

λi ~ Gamma(shape, rate)

ポアソンーガンマ混合分布は，負の二項分布と等価である．Ｒで負の二項分布の乱数を発生させたければ，以下のようにしたら良い．

n <- 10000 # 発生させる個数

rpois(n, lambda = rgamma(n, shape = 10, rate = 0.1)) # rate = 1/scale

gamma分布は，定常ポアソン過程の複数回のスパイク間隔の一般化で，rateは強度パラメータ，shapeは「何個分のスパイクの待ち時間か」に相当する．例えば上の例だと，単位時間当たり0.1回生じるイベントが10回生じるまでの時間になる．詳しくは，このサイトが丁寧で分かりやすい．

shapeとrateを用いた負の二項分布の期待値は以下の通りになる．

期待値：shape * 1 / rate　# ガンマ分布の期待値と当然同じ

分散：shape * (1 / rate)^2 * (1 + rate)  # ガンマ分布の分散に (1 + rate)掛けたもの

shapeは，dispersion parameterと呼ばれることが多い（値が小さいほど，バラつきは大きくなるので注意）．

2. RデフォルトのXnbinom関数の場合

　この関数で乱数n個を発生させる場合は以下のようになる．

rnbinom(n, size = 10, prob = 10/11) # 次の3で出てくるmuによる指定も可

　負の二項分布の説明には，①成功確率pのもとで，ベルヌーイ試行が k 回成功するまでの失敗回数 x の分布であるというものと，②成功確率pのもとで，ベルヌーイ試行が k 回成功するまでに必要な試行回数 x の分布であるという，２通りが存在する．この関数は，前者①の方に従ってパラメータをとっている．

　期待値は以下の通りになる．

期待値：size * (1 - prob) / prob

分散：size * (1 - prob) / prob^2

3.  MASSパッケージのXnegbin関数の場合

　この関数で乱数n個を発せさせる場合は以下のようになる．

rnegbin(n, mu = 100, theta = 10) 

　これは素直に期待値muとDispersion parameterであるthetaをとる．期待値は以下の通りになる．

期待値：mu

分散：mu^2/theta

4. 上の３つのパラメータ間の関係の整理

　作りたかったのはこれ．パラメータの変換はこれで完璧！

期待値muと分散varからshape = size = thetaを知りたい場合は，

theta = mu ^ 2 / var

5.  glm.nb関数の場合

　推定されるのは，平均の期待値muとDispersion parameterであるtheta．対応を考えると乱数発生をさせる場合は，3.  MASSパッケージのrnegbin関数が楽．例えば，

library(MASS)

r <- rnegbin(n, mu= 100, theta = 10) 

model <- glm.nb(r ~ 1)

exp(model$coefficients) # muの推定

model$theta # thetaの推定

6. JAGS

　混合分布ではなく，負の二項分布自体を使ってモデリングする場合は，基本glm.nbと同様（勘違いしてた。近日中に更新→更新済み）。

JAGSSの定義は，

for(i in 1:n){

    y[i] ~ dnegbin(p, size) 

}

Rのデフォルトと同じ

7. Stan

  for (i in 1:N) {

    y ~ neg_binomial(shape, rate);

  }

これは，ポアソン‐ガンマ混合分布のshapeとrateに相当．

ってか，ダルい・・・

8. 以下のサイトが役に立つ．

http://nhkuma.blogspot.com/2014/02/r-winbugs.html?m=1

https://ito-hi.blog.ss-blog.jp/2014-02-11

https://www.slideshare.net/simizu706/ss-50994149