Loading [MathJax]/jax/output/HTML-CSS/jax.js

コインを10回投げて表の出る真の確率を推定する

コインを10回投げて,表を1裏を0とする以下のデータを得た.

data<-c(0,1,1,1,1,0,1,1,0,1)

表の出る真の確率はどれくらいだろうか? 

ここでは,独立同分布のベルヌイ分布を確率モデル,ベータ分布を事前分布としたベイズ統計モデルを構築し,ベルヌイ分布のパラメータqの事後分布をMCMCによって導出することを目標とする.

x1,x2,,xniidp(x|q) p(xi|q)=qxi(1q)1xi qp(q|α,β) p(q|α,β)=qα1(1q)β1B(α,β)

尤度・対数尤度

まず,尤度の関数を実装する.(パラメタqについての条件付き)独立同分布の仮定より尤度は, p({xi}|q)=ni=1p(xi|q)=qxi(1q)nxi

となる.

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)=ni=1logp(xi|q)=ni=1xilogq+(nni=1xi)log(1q)

となる.

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サンプリングのアルゴリズムとしては,もっとも単純なメトロポリス・アルゴリズムを採用する.その手順は以下のとおり.

  1. 現在のパラメタの値θ0を1つサンプリングする.
  2. ランダムに別のパラメタの値θ1をとってくる.
  3. θ0θ1の事後確率の比(事後オッズ)を計算する. PO(θ0,θ1)=p(θ1|x)p(θ0|x)=p(x|θ1)p(x|θ0)p(θ1)p(θ0)
  4. 確率pmove=Min(1,PO(θ0,θ1))θ1に移動する.確率1pmoveθ0にとどまる.
  5. 新しいパラメタの値をθ0として,(1) に戻る.

詳細は,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")

おいしくできたね!