- 一昨日の記事でRのDPpackageに触れた
- 今日の記事では、二項分布・ベータ分布を例にいじってみる
- 関数はDPbetabinom()である
- まずは、この関数のヘルプ記事のExamplesから
library(DPpackage)
data(rolling)
y <- cbind(rolling$y1,rolling$y2)
prior<-list(alpha=1,
a1=1,
b1=1)
state <- NULL
mcmc <- list(nburn=5000,
nsave=10000,
nskip=3,
ndisplay=100)
fit <- DPbetabinom(y=y,ngrid=100,
prior=prior,
mcmc=mcmc,
state=state,
status=TRUE)
fit
summary(fit)
plot(fit,output="density")
plot(fit,output="param")
head(y)
-
- 以下の通り、9回実行して(第2列)、成功した回数(第1列)を記録した行列
> 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回試行毎に、(コイン投げなら)コインが変わっているとして、そのコインには、表が出る確率がある。ただし、このコインはいくつかの群に分けられて、群ごとに成功率が決まっていると考えたい。ただし、何番目のコインがどの群に属するかはわからないし、何番目のコインと何番目のコインが同じ群かもわからないし、群の数もわからないし、群に属するコインの個数(割合)もわからないとする
- この解析のモデルは以下の通り
- ここで、観測データとして提供されるのは

- ほかに指定するべきは
を作るための情報
- その情報は
と
であるが、
は
を必要とする。
はその値を決め打ちで与えても大丈夫だが、
を与えることでも生成できることが読み取れる。
- 結局、
の3つを与えるか、
を与えるかの2通りの実行方法があることとなる。
- 二通りのprior引数の記載例を以下に示す。
- MCMCの条件を改変して得られる事後分布が適当でないながら、ばらつきがよく見えるような条件で実施することにする。
prior3<-list(alpha=1,
a1=1,
b1=1)
prior4<-list(a0=1,
b0=1,
a1=1,
b1=1)
mcmc2 <- list(nburn=500,
nsave=1000,
nskip=3,
ndisplay=100)
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)
-
- 何が違うか、というと、MCMCを回して記録されるのはクラスタ数、
、LPML(尤度)であるが、prior3の場合は
値を固定しているので、ばらつきがないこととなる
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)

- 事後分布
- 本当に知りたいのは「生起確率
の分布。もしクラスタ数を1に限ると、その値がどのあたりか、という評価になるし、クラスタ数を非常に大きくすれば、非常に自由度の高い
の分布になる
- それ以外に知りたい・評価するのが適当なのは、クラスタ数がいくつなのかということと、
の値がいくつなのか、ということ。また、このクラスタ数、
値は、初期値依存な状態から、MCMCで探索しながら、それなりの分布に落ち着くべきなので、クラスタ数・
値はその分布と、MCMCプロセスの推移とを評価したい。また尤度の推移も。
- それを表示するのが以下の関数である。
plot(fit4,output="density")
plot(fit4,output="param")
-
- このplot()関数は使い勝手があまり良くないので、少し変える。
plot(fit4$grid,fit4$densp.m,type="l")

plot(fit4$save.state$thetasave[,1],type="l")
hist(fit4$save.state$thetasave[,1])


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)
fit5 <- DPbetabinom(y=y,ngrid=100,
prior=prior5,
mcmc=mcmc2,
state=state,
status=TRUE)
plot(fit5$grid,fit5$densp.m,type="l")
plot(fit5$save.state$thetasave[,1],type="l")
hist(fit5$save.state$thetasave[,1])
-
- 特定のpに鋭いピークがあり、クラスタ数も一貫して1となっている。


- クラスタ数が1になるよう、
を小さくしたままに、
の分布として、特定の値に集中したベータ分布を使うと、観察データによってもそこからずれることができないことが以下の結果から読み取れる。
fit6 <- DPbetabinom(y=y,ngrid=100,
prior=prior6,
mcmc=mcmc2,
state=state,
status=TRUE)
plot(fit6$grid,fit6$densp.m,type="l")

- 逆に
を大きくして、クラスタ数がものすごく大きくなるようにしてやってみる。まずは、
の分布は一様分布としてやってみる。
- 時間がかかることに注意する。新規のクラスタからのサンプルであるとみなされる場合が増え、自由度が増すことに対応する時間の増大。
prior7<-list(alpha=100,
a1=1,
b1=1)
fit7 <- DPbetabinom(y=y,ngrid=100,
prior=prior7,
mcmc=mcmc2,
state=state,
status=TRUE)
plot(fit7$grid,fit7$densp.m,type="l")
plot(fit7$save.state$thetasave[,1],type="l")
hist(fit7$save.state$thetasave[,1])


の結果を比べてみる
prior8<-list(alpha=1,
a1=1,
b1=1)
fit8 <- DPbetabinom(y=y,ngrid=100,
prior=prior8,
mcmc=mcmc2, しよう
state=state,
status=TRUE)
plot(fit8$grid,fit8$densp.m,type="l")
-
- クラスタ数がぐっと少ない方に寄り、その結果、pの分布がクラスタが作るピークを示す様態に変わっている。推定クラスタ数分布では、4〜6クラスタ辺りが多いが、pの分布が2ピークになっているのは、2つのクラスタの割合が大きいがそれ以外のクラスタに対応する峰は「飲み込まれて」しまっているから。



-
- クラスタ数が少ないはず、と信じれば、pは多峰性であると推定するのがよいし、クラスタ数がとても多いと信じれば、pはなめらかな分布になると推定するのがよい、ということを示している。
- じゃあ、どっちが本当なの?とい場合にはどうしたらよいのか…
の値が不明なので、それを確率変数にするというのがベイズらしいやり方。
- DPbetabinom()でもそれが用意されている