2017-07-13

ディリクレ過程で二項分布

一昨日の記事で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回試行毎に、(コイン投げなら)コインが変わっているとして、そのコインには、表が出る確率がある。ただし、このコインはいくつかの群に分けられて、群ごとに成功率が決まっていると考えたい。ただし、何番目のコインがどの群に属するかはわからないし、何番目のコインと何番目のコインが同じ群かもわからないし、群の数もわからないし、群に属するコインの個数(割合)もわからないとする
この解析のモデルは以下の通り
- $y_i \sim Binom(n_i,p_i)$
- $p_i \sim DP(\alpha, G_0)$
- $G0 = Beta(a_1,b_1)$
- $\alpha \sim Gamma(a_0,b_0)$
ここで、観測データとして提供されるのは $y,n$
ほかに指定するべきは $p_i$ を作るための情報
その情報は $\alpha$ と $G_0$ であるが、 $G_0$ は $a_0,b_0$ を必要とする。
$\alpha$ はその値を決め打ちで与えても大丈夫だが、 $a_1.b_1$ を与えることでも生成できることが読み取れる。
結局、の３つを与えるか、を与えるかの２通りの実行方法があることとなる。
- 二通りの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)

- 何が違うか、というと、MCMCを回して記録されるのはクラスタ数、 $\alpha$ 、LPML(尤度)であるが、prior3の場合は $\alpha$ 値を固定しているので、ばらつきがないこととなる

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)

- 中央列が $\alpha$ のヒストグラム

事後分布
- 本当に知りたいのは「生起確率 $p$ の分布。もしクラスタ数を１に限ると、その値がどのあたりか、という評価になるし、クラスタ数を非常に大きくすれば、非常に自由度の高い $p$ の分布になる
- それ以外に知りたい・評価するのが適当なのは、クラスタ数がいくつなのかということと、 $\alpha$ の値がいくつなのか、ということ。また、このクラスタ数、 $\alpha$ 値は、初期値依存な状態から、MCMCで探索しながら、それなりの分布に落ち着くべきなので、クラスタ数・ $\alpha$ 値はその分布と、MCMCプロセスの推移とを評価したい。また尤度の推移も。
- それを表示するのが以下の関数である。

# 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峰性にならざるを得ない。それは $\alpha$ を小さくすればよい。その値として何が選ばれるべきかをフラットに考えるならば $a_1=b_1=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となっている。

クラスタ数が１になるよう、 $\alpha$ を小さくしたままに、 $G_0$ の分布として、特定の値に集中したベータ分布を使うと、観察データによってもそこからずれることができないことが以下の結果から読み取れる。

# 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の分布はなめらか

クラスタ数の分布。クラスタ数はとても多く、100超え

$\alpha=1,100$ の結果を比べてみる

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")

- クラスタ数がぐっと少ない方に寄り、その結果、pの分布がクラスタが作るピークを示す様態に変わっている。推定クラスタ数分布では、４〜６クラスタ辺りが多いが、pの分布が２ピークになっているのは、２つのクラスタの割合が大きいがそれ以外のクラスタに対応する峰は「飲み込まれて」しまっているから。

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