読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

久保(2012)マルコフ連鎖モンテカルロ(MCMC)法入門のための例題

概要

生態学データ解析 - 本/データ解析のための統計モデリング入門

『データ解析のための統計モデリング入門』第 8 章の例題を R でやる.

ふつうの最尤推定

二項分布

\displaystyle P(X=k)={n\choose k}p^k(1-p)^{n-k}\quad\mbox{for}\ k=0,1,2,\dots,n

のパラメータ, 成功確率 p最尤推定量は解析的にかんたんに求まる.

サンプルサイズを m とすると

\displaystyle \hat p = \frac{\sum x_i}{nm}

N <- 8
dat <- c(4,3,4,5,5,2,3,1,4,0,1,5,5,6,5,4,4,5,3,4)

phat <-sum(dat)/(N*length(dat))
phat
#[1] 0.45625

plot(table(dat),xlab="",ylab="")
points(0:6,dbinom(0:6,N,phat)*length(dat),pch=4,type="b",cex=2,lty=2)

f:id:abrahamcow:20150725090752p:plain
図は例題のデータを集計した棒グラフ. 破線は二項分布の当てはめ.

ふらふら試行錯誤による最尤推定

効率が悪く, 精度もよくない方法.

  1. 成功確率を離散化する(ここでは p=0.01 から 0.99 まで 0.01 刻みの値をとることにする)
  2. パラメータ q の初期値を選ぶ
  3. q を増やすか減らすかをランダムに決める(新しく選んだ q の値を q' とする)
  4. q’ において尤度が大きくなるなら q の値をq'に変更する
  5. 3, 4 を繰り返す
L <- function(p){
  sum(dbinom(dat, N, p, log=TRUE))
}

q <- seq(0.01,0.99,by=0.01) #離散化
lv <- sapply(q,L)
#nihongo()
plot(q,lv, xlab="成功確率 p", ylab="対数尤度")

#おまけ:最尤推定
#res <- optimize(L,c(0,1), maximum=TRUE)

f:id:abrahamcow:20150725090949p:plain
図は離散化した生存確率と対応する対数尤度.

furafura <- function(qp0){
  pv <- numeric(100)
  for(i in 1:100){
    pv[i] <- q[qp0]
    lv0 <- lv[qp0]
    qp1 <- qp0 +sample(c(1,-1),1)
    lv1 <- lv[qp1]
    if(lv1 > lv0){
      qp0 <- qp1
    }
  }
  pv
}

pv1 <- furafura(30)
pv2 <- furafura(60)
matplot(cbind(pv1,pv2),type="l",lwd=2,xlab="試行錯誤のステップ数",ylab="成功確率 p")

f:id:abrahamcow:20150725091235p:plain
図は試行錯誤による p の変化. 黒い実線は p=0.30 からスタートした試行錯誤. 赤い破線は p = 0.60 からスタートした試行錯誤. 0.46 付近に収束していく.

MCMCアルゴリズムのひとつ:メトロポリス

ふらふら試行錯誤の手順に,

  • q’ で尤度が小さくなる場合であっても,確率 rq の値を変更する

を追加する.
ここで r は尤度比 L(q’)/ L(q)

MH <- function(iter){
  pv <- numeric(iter)
  for(i in 1:iter){
    pv[i] <- q[qp0]
    lv0 <- lv[qp0]
    qp1 <- qp0 +sample(c(1,-1),1)
    lv1 <- lv[qp1]
    if(lv1 > lv0){
      qp0 <- qp1
    }else{
      r <- exp(lv1-lv0) #対数尤度比を尤度比に
      if(runif(1) < r) qp0 <- qp1
    }
  }
  pv
}
qp0 <- 30
pv1 <- MH(100)
plot(pv1,type="l")

plotfunc <- function(pv1){
  #デフォルトを保存
  def.par <- par(no.readonly = TRUE)
  nf <- layout(matrix(c(1,2),2,2,byrow = TRUE), c(3,1), c(1,3), TRUE)
  
  # plot1
  plot(pv1, type="l",xlab="MCMC step 数", ylab="成功確率 p")
  
  #plot2
  par(mai = c(0.85, 0, 0.68, 0.35))
  hist1 <- hist(pv1, plot=FALSE)
  barplot(hist1$count, horiz=TRUE, space=0)
  
  #もとにもどす
  par(def.par)
}
pv1 <- MH(100)
plotfunc(pv1)

pv1 <- MH(1000)
plotfunc(pv1)

pv1 <- MH(10000)
plotfunc(pv1)

f:id:abrahamcow:20150725091149p:plain
f:id:abrahamcow:20150725091210p:plain
f:id:abrahamcow:20150725091216p:plain
図はメトロポリス法による MCMC サンプリングの例. 定常分布に収束していく様子がわかる.