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

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

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

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

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

\[x_1,x_2,\ldots,x_{n} \sim_{iid} p(x|q) \] \[p(x_i|q) = q^{x_i}(1-q)^{1-x_i}\] \[q\sim p(q|\alpha,\beta)\] \[p(q|\alpha,\beta)=\frac{q^{\alpha-1}(1-q)^{\beta-1}}{B(\alpha,\beta)}\]

尤度・対数尤度

まず，尤度の関数を実装する．（パラメタ\(q\)についての条件付き）独立同分布の仮定より尤度は， \[\begin{align*} p(\{x_i\}|q)&=\prod_{i=1}^n p(x_i|q) \\ &=q^{\sum x_i}(1-q)^{n-\sum x_i} \end{align*}\]

となる．

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

対数尤度は， \[\begin{align*} \log p(\{x_i\}|q)&=\sum_{i=1}^n \log p(x_i|q) \\ &=\sum_{i=1}^n x_i \log q+ \left(n-\sum_{i=1}^n x_i\right)\log(1-q) \end{align*}\]

となる．

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|\{x_i\})\propto p(\{x_i\}|q)p(q)=p(\{x_i\},q) \] という比例関係を利用して，データとパラメタの同時分布（のデータを\(\{x_i\}\)に固定した関数）\(p(\{x_i\},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. 現在のパラメタの値\(\theta_{0}\)を1つサンプリングする．
  2. ランダムに別のパラメタの値\(\theta_{1}\)をとってくる．
  3. \(\theta_{0}\)と\(\theta_{1}\)の事後確率の比（事後オッズ）を計算する． \[PO(\theta_{0},\theta_{1})=\frac{p(\theta_{1}|x)}{p(\theta_{0}|x)}=\frac{p(x|\theta_{1})}{p(x|\theta_{0})}\frac{p(\theta_{1})}{p(\theta_{0})}\]
  4. 確率\(p_{move}=\mathrm{Min}(1,PO(\theta_{0},\theta_{1}))\)で\(\theta_{1}\)に移動する．確率\(1-p_{move}\)で\(\theta_{0}\)にとどまる．
  5. 新しいパラメタの値を\(\theta_{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")

おいしくできたね！