階層離散分布のギブスサンプラー

※古い記事です。懐かしのBRugsを使っています。

問題設定

  • ２つのミカンの産地$ \tau_1$ と $\tau_2$ を考える。
  • ミカンには、甘い$y_1$ のと酸っぱい $y_2$ のがあり、年によってその割合は異なる。
  • 良作の年 $\beta_1$ は9割が甘いが、通常の年 $\beta_2$ は半々ぐらい、不作の年 $\beta_3$ は、1割しか甘いのがない。
  • 産地$\tau_1$は良作の年が多く、産地$\tau_2$は不作の年が多い。
  • 今、店で5つのミカンを買って食べたところ、｛甘い、甘い、酸っぱい、酸っぱい、甘い｝、すなわち、$y={y_1,y_1,y_2,y_2,y_1}$ だった。
  • どちらの産地のどの年のミカンが出回っているのだろうか？

尤度 $Pr(y|\tau,\beta)$ は次のとおり。

$y=1$ $y=2$
$\Pr(\tau_1)=0.5$ $\Pr(\beta_1;\tau_1)=0.7$ 0.9 0.1
$\Pr(\beta_2;\tau_1)=0.2$ 0.5 0.5
$\Pr(\beta_2;\tau_1)=0.1$ 0.1 0.9
$\Pr(\tau_2)=0.5$ $\Pr(\beta_1;\tau_2)=0.1$ 0.9 0.1
$\Pr(\beta_2;\tau_2)=0.2$ 0.5 0.5
$\Pr(\beta_3;\tau_1)=0.7$ 0.1 0.9

求めたい事後分布は、$\beta_1,\tau_1$ で例示的に示すと、

$$
\Pr(\beta_1,\tau_1|y)
=\frac{\Pr(y|\beta_1)\Pr(\beta_1|\tau_1)\Pr(\tau_1)}
{\sum_{i=1}^{3}\Pr(y|\beta_i)\Pr(\beta_i|\tau_1)\Pr(\tau_1)+\sum_{i=1}^{3}\Pr(y|\beta_i)\Pr(\beta_i|\tau_2)\Pr(\tau_2)}
$$

ただし、

$$
\Pr(y|\beta_1)=\Pr(y_1|\beta_1)\Pr(y_1|\beta_1)\Pr(y_2|\beta_1)\Pr(y_2|\beta_1)\Pr(y_1|\beta_1)
$$

$\tau$ を与えたときの $\beta$ の事後分布は、

$$
\Pr(\beta_1|\tau_1,y)
=\frac{\Pr(y|\beta_1)\Pr(\beta_1|\tau_1)\underline{\Pr(\tau_1)}}
{\sum_{i=1}^{3}\Pr(y|\beta_i)\Pr(\beta_i|\tau_1)\underline{\Pr(\tau_1)}}
$$

$$
=\frac{\Pr(y|\beta_1)\Pr(\beta_1|\tau_1)}
{\sum_{i=1}^{3}\Pr(y|\beta_i)\Pr(\beta_i|\tau_1)}
$$

$\beta$ を与えたときの $\tau$ の事後分布は、

$$
\Pr(\tau_1|\beta_1,y)
=\frac{\underline{\Pr(y|\beta_1)}\Pr(\beta_1|\tau_1)\Pr(\tau_1)}
{\underline{\Pr(y|\beta_1)}\Pr(\beta_1|\tau_1)\Pr(\tau_1)+
\underline{\Pr(y|\beta_1)}\Pr(\beta_1|\tau_2)\Pr(\tau_2)}
$$

$$
=\frac{\Pr(y|\beta_1)\Pr(\beta_1|\tau_1)\Pr(\tau_1)}
{\Pr(\beta_1|\tau_1)\Pr(\tau_1)+
\Pr(\beta_1|\tau_2)\Pr(\tau_2)}
$$

ギブスサンプラーは、上記の２種の確率で$\beta$　と　$\tau$　を発生させること。

GSのプロシジャーをRで再現してみる

各種確率を与える。

py<-rep(NA,3)
pb<-matrix(nrow=3,ncol=2)
pt <-rep(NA,2)
py[1] <-0.9
py[2] <-0.5
py[3] <-0.1
p <-cbind(py,1-py)

pb[1,1] <-0.7
pb[2,1] <-0.2
pb[3,1] <-0.1

pb[1,2] <-0.1
pb[2,2] <-0.2
pb[3,2] <-0.7

pt[1] <-0.5
pt[2] <-0.5

y <-c(1,1,2,2,1)

$\Pr(y|\beta)$

p
      py
[1,] 0.9 0.1
[2,] 0.5 0.5
[3,] 0.1 0.9

pb
     [,1] [,2]
[1,]  0.7  0.1
[2,]  0.2  0.2
[3,]  0.1  0.7

$\Pr(\tau)$

pt
[1] 0.5 0.5

$y$

y
[1] 1 1 2 2 1

τを与えたときのβの事後分布を求める関数の定義
$\Pr(β=1,2,3|\tau,y)$
$=\Pr(y|\beta=1,2,3)\Pr(β=1,2,3|\beta)/\sum_i\Pr(y|\beta=i)\Pr(\beta=i|\tau)$

pbeta <-function(y,tau){
    N <-length(y)
    p01 <-c(1,1,1)
    for(i in 1:N){
        p01 <-p01*p[,y[i]]
    }
    pb01 <-p01*pb[,tau]/
          (p01[1]*pb[1,tau]
          +p01[2]*pb[2,tau]
          +p01[3]*pb[3,tau])
    pb01
}

$\beta$を与えたときの$tau$の事後分布を求める関数の定義
$\Pr(τ=1,2|\beta,y)$
$=\Pr(\beta|\tau=1,2)\Pr(\tau=1,2)$

ptau <-function(beta){
    pt00 <-pb[beta,]*pt[]/
          (pb[beta,1]*pt[1]
          +pb[beta,2]*pt[2])
    pt00
}

βとτの初期値を(1,1)に指定。

beta00 <-1
tau00 <-1

繰り返しの回数を10000に設定

itr=10000

ギブスサンプラーの実行：o01に結果。

o01 <-matrix(nrow=itr,ncol=2)
for(i in 1:itr){
    beta01 <-sample(1:3,1,p=pbeta(y,tau00))
    tau01 <-sample(1:2,1,p=ptau(beta00))
    beta00 <-beta01
    tau00 <-tau01
    o01[i,1] <-beta01
    o01[i,2] <-tau01
}
hist(o01[,1])

BRugsでも計算してみる。

gmodel.txt

model{
for(i in 1:N){
y[i]~dcat(p[])
}

p[1] <-py[beta]
p[2] <-1-p[1]

py[1] <-0.9
py[2] <-0.5
py[3] <-0.1

pb[1,1] <-0.7
pb[2,1] <-0.2
pb[3,1] <-0.1

pb[1,2] <-0.1
pb[2,2] <-0.2
pb[3,2] <-0.7

pt[1] <-0.5
pt[2] <-0.5

beta~dcat(pb[,tau])
tau~dcat(pt[])
}

Rスクリプト

d01 <-list(y=c(1,1,2,2,1),N=5)
d02 <-list(list(beta=1,tau=1),list(beta=2,tau=2))

library(BRugs)

setwd("d:/home/doc")

modelCheck("gmodel.txt")
bugsData(d01,"./data.txt")
modelData("./data.txt")
modelCompile(numChains=2)
bugsInits(d02,numChain=2,file=c("./inits1.txt","./inits2.txt"))
modelInits(c("inits1.txt","inits2.txt"))
modelGenInits()
samplesSet(c("beta", "tau"))
summarySet(c("alpha","beta1"))
modelUpdate(10000)
samplesStats("*")

両者の結果を比較してみる。

Rスクリプト

summary(factor(o01[,1]))/length(o01[,1])
o02 <-samplesSample("beta")
summary(factor(o02))/length(o02)
summary(factor(o01[,2]))/length(o01[,2])
o03 <-samplesSample("tau")
summary(factor(o03))/length(o03)

結果

 summary(factor(o01[,1]))/length(o01[,1])
       1      2      3
  0.3073 0.6616 0.0311
o02<-samplesSample("beta")
summary(factor(o02))/length(o02)
          1          2          3
 0.31548155 0.64996500 0.03455346
summary(factor(o01[,2]))/length(o01[,2])
      1      2
 0.6043 0.3957
 o03<-samplesSample("tau")
 summary(factor(o03))/length(o03)
         1         2
 0.6079108 0.3920892

産地１の並作の年のミカンである確率が最も高い。それぞれ6割、6割5分で。