ディリクレ過程で二項分布
- 一昨日の記事でRのDPpackageに触れた
- 今日の記事では、二項分布・ベータ分布を例にいじってみる
- 関数はDPbetabinom()である
- まずは、この関数のヘルプ記事のExamplesから
- 全部コピペして動かしてみる
library(DPpackage) # Data data(rolling) y <- cbind(rolling$y1,rolling$y2) # Prior information prior<-list(alpha=1, a1=1, b1=1) # Initial state state <- NULL # MCMC parameters mcmc <- list(nburn=5000, nsave=10000, nskip=3, ndisplay=100) # Fitting the model fit <- DPbetabinom(y=y,ngrid=100, prior=prior, mcmc=mcmc, state=state, status=TRUE) fit summary(fit) # density estimate plot(fit,output="density") # parameters plot(fit,output="param")
- 観測データは
head(y)
-
- 以下の通り、9回実行して(第２列)、成功した回数(第１列)を記録した行列
> head(y) [,1] [,2] [1,] 7 9 [2,] 4 9 [3,] 6 9 [4,] 6 9 [5,] 6 9 [6,] 6 9
hist(y[,1])
- 今、考えたいのは、9回試行毎に、(コイン投げなら)コインが変わっているとして、そのコインには、表が出る確率がある。ただし、このコインはいくつかの群に分けられて、群ごとに成功率が決まっていると考えたい。ただし、何番目のコインがどの群に属するかはわからないし、何番目のコインと何番目のコインが同じ群かもわからないし、群の数もわからないし、群に属するコインの個数(割合)もわからないとする
- この解析のモデルは以下の通り
- ここで、観測データとして提供されるのは
- ほかに指定するべきは
を作るための情報
- その情報は
と
であるが、
を必要とする。
を与えることでも生成できることが読み取れる。
- 結局、
の３つを与えるか、
を与えるかの２通りの実行方法があることとなる。
- 二通りのprior引数の記載例を以下に示す。
- MCMCの条件を改変して得られる事後分布が適当でないながら、ばらつきがよく見えるような条件で実施することにする。
prior3<-list(alpha=1, a1=1, b1=1) prior4<-list(a0=1, b0=1, a1=1, b1=1) # MCMC parameters mcmc2 <- list(nburn=500, nsave=1000, nskip=3, ndisplay=100) # Fitting the model fit3 <- DPbetabinom(y=y,ngrid=100, prior=prior3, mcmc=mcmc2, state=state, status=TRUE) fit4 <- DPbetabinom(y=y,ngrid=100, prior=prior4, mcmc=mcmc2, state=state, status=TRUE)
tasave,2,range) ncluster alpha LPML [1,] 2 1 -2562.854 [2,] 13 1 -2548.604 > apply(fit4$save.state$thetasave,2,range) ncluster alpha LPML [1,] 2 0.03358158 -2856.759 [2,] 20 6.96844599 -2409.200 par(mfrow=c(2,3)) apply(fit3$save.state$thetasave,2,hist) apply(fit4$save.state$thetasave,2,hist)
-
- 中央列が
- 中央列が
- 事後分布
# density estimate(pの値の分布) plot(fit4,output="density") # parameters(クラスタ数、alpha値の分布と推移、尤度の推移) plot(fit4,output="param")
-
- このplot()関数は使い勝手があまり良くないので、少し変える。
# density estimate(pの値の分布) plot(fit4$grid,fit4$densp.m,type="l")
# No. clusters の推移と分布 plot(fit4$save.state$thetasave[,1],type="l") hist(fit4$save.state$thetasave[,1])
# alpha の推移と分布 plot(fit4$save.state$thetasave[,2],type="l") plot(density(fit4$save.state$thetasave[,2]))
# 尤度の推移 plot(fit4$save.state$thetasave[,3],type="l")
- クラスタ数を限りなく1にすると、推定されるpの分布は1峰性にならざるを得ない。それは
が一様分布に対応する。やってみる。
# 条件を厳しくして回してみる prior5<-list(alpha=10^(-4), a1=1, b1=1) # Fitting the model fit5 <- DPbetabinom(y=y,ngrid=100, prior=prior5, mcmc=mcmc2, state=state, status=TRUE) # pの分布 plot(fit5$grid,fit5$densp.m,type="l") # No. clusters の推移と分布 plot(fit5$save.state$thetasave[,1],type="l") hist(fit5$save.state$thetasave[,1])
-
- 特定のpに鋭いピークがあり、クラスタ数も一貫して1となっている。
- クラスタ数が１になるよう、
# Fitting the model fit6 <- DPbetabinom(y=y,ngrid=100, prior=prior6, mcmc=mcmc2, state=state, status=TRUE) # pの分布 plot(fit6$grid,fit6$densp.m,type="l")
- 逆に
- 時間がかかることに注意する。新規のクラスタからのサンプルであるとみなされる場合が増え、自由度が増すことに対応する時間の増大。
# クラスタ数を巨大に prior7<-list(alpha=100, a1=1, b1=1) # Fitting the model fit7 <- DPbetabinom(y=y,ngrid=100, prior=prior7, mcmc=mcmc2, state=state, status=TRUE) # pの分布 plot(fit7$grid,fit7$densp.m,type="l") # No. clusters の推移と分布 plot(fit7$save.state$thetasave[,1],type="l") hist(fit7$save.state$thetasave[,1])
-
- pの分布はなめらか
の結果を比べてみる
prior8<-list(alpha=1, a1=1, b1=1) # Fitting the model fit8 <- DPbetabinom(y=y,ngrid=100, prior=prior8, mcmc=mcmc2, しよう state=state, status=TRUE) # pの分布 plot(fit8$grid,fit8$densp.m,type="l")
