コインを10回投げて,表を1裏を0とする以下のデータを得た.
data<-c(0,1,1,1,1,0,1,1,0,1)
表の出る真の確率はどれくらいだろうか?
ここでは,独立同分布のベルヌイ分布を確率モデル,ベータ分布を事前分布としたベイズ統計モデルを構築し,ベルヌイ分布のパラメータqの事後分布をMCMCによって導出することを目標とする.
x1,x2,…,xn∼iidp(x|q) p(xi|q)=qxi(1−q)1−xi q∼p(q|α,β) p(q|α,β)=qα−1(1−q)β−1B(α,β)
となる.
LL_Bern<-function(x,q) {
q^sum(x)*(1-q)^(length(x)-sum(x))
}
plot(seq(0,1,0.01),LL_Bern(data,seq(0,1,0.01)),type="l")
対数尤度は,
logp({xi}|q)=n∑i=1logp(xi|q)=n∑i=1xilogq+(n−n∑i=1xi)log(1−q)
となる.
logLL_Bern<-function(x,q) {
sum(x)*log(q)+(length(x)-sum(x))*log(1-q)
}
plot(seq(0,1,0.01),logLL_Bern(data,seq(0,1,0.01)),
type="l")
尤度・対数尤度の最大化問題を解くと,最尤推定値が得られる.
optimize(function(q) LL_Bern(data,q),c(0,1),maximum = TRUE)
## $maximum
## [1] 0.6999843
##
## $objective
## [1] 0.002223566
optimize(function(q) logLL_Bern(data,q),c(0,1),maximum = TRUE)
## $maximum
## [1] 0.7000058
##
## $objective
## [1] -6.108643
事後分布を推定するために p(q|{xi})∝p({xi}|q)p(q)=p({xi},q) という比例関係を利用して,データとパラメタの同時分布(のデータを{xi}に固定した関数)p({xi},q)からのサンプリングを行う.
ここで,具体的な事前分布としてp(q|1,1)を仮定する.これは区間[0,1]上での一様分布と等しい.つまり,p(q|1,1)=1.ここでは,今後の拡張のために一般的な形で関数を定義しておく.
prior_beta<-function(q,a,b) dbeta(q,a,b) #same as dunif(q,0,1)
joint<-function(x,q) LL_Bern(x,q)*prior_beta(q,1,1)
plot(seq(0,1,0.01),joint(data,seq(0,1,0.01)),
type="l")
MCMCサンプリングのアルゴリズムとしては,もっとも単純なメトロポリス・アルゴリズムを採用する.その手順は以下のとおり.
詳細は,DBDA, 2nd ed. ch.7 (Kruschke 2014),三中2018『統計思考の世界』13講などを参照.
Metropolis<-function(current) {
propose<-runif(1,0,1)
postOdds<-joint(data,propose)/joint(data,current)
pmove<-min(postOdds,1)
if(pmove>=runif(1,0,1)) propose else current
}
初期値を0.5として,11000回サンプリングする.
nsteps<-11000 # number of steps
mcmcsample<-rep(NA,nsteps + 1)
mcmcsample[1]<-0.5 # initial position
for (i in 1:nsteps) {
mcmcsample[i+1]<-Metropolis(mcmcsample[i])
}
最初の1000回分をバーンインとして削除して,残りの10000回をプロット.
plot(mcmcsample[-(1:1000)],type="l",col="skyblue")
hist(mcmcsample[-(1:1000)],freq =FALSE,col="skyblue")
lines(seq(0,1,0.01),
dbeta(seq(0,1,0.01),1+sum(data),1+length(data)-sum(data)),col="red")
おいしくできたね!