久保みどり本 8章 今回は手を動かすというよりも読み物メイン。

midori_study_8.r

library("ggplot2")
library("magrittr")
library("dplyr")

# midori 8
# http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html
data <- c(4,3,4,5,5,2,3,1,4,0,1,5,5,6,5,4,4,5,3,4)
hist(data)

loglik <- function(q){
  # 二項分布対数尤度 logL(q)
  # 20個体ぶんの「データが得られる確率」の積
  sum(data * log(q) + (8 - data) * log(1-q) + log(choose(8,data)))
}

loglik(0.29)
loglik(0.30)
loglik(0.31)

getsteps <-function(init.q,step=100,alpha=0.01){
  iter.length <- step
  q <- init.q
  x <- loglik(q)
  for(i in 2:iter.length){
    sign <- sample(c(alpha,-alpha),1)
    if(x[i-1] < loglik(q[i-1]+sign)){
      q[i] <- q[i-1]+sign
      x[i] <- loglik(q[i-1]+sign)
    }else{
      q[i] <- q[i-1]
      x[i] <- x[i-1]
    }
  }
  data.frame(init.q,step=1:step,q,loglik=x)
}
q1 <- getsteps(init.q = 0.5,alpha=.01)
q2 <- getsteps(init.q = 0.3,alpha=.01)
qplot(data=rbind(q1,q2),x=step,y=q,col=factor(init.q))+geom_line() # p.174 8.5

# p.177 metropolis
metropolis<-function(init.q,step=100,alpha=0.01){
  iter.length <- step
  q <- init.q
  x <- loglik(q)
  for(i in 2:iter.length){
    # 両隣を更新先の候補とする
    sign <- sample(c(alpha,-alpha),1) 

    # rの確率（尤度比）で値qを選択。r = 尤度_新/尤度_前 = exp(尤度_新 - 尤度_前)
    # 前回よりも尤度が大きければ必ずそちらを選択するが
    # 尤度が小さくても(r<1)の場合でもrの確率で選択する場合がある。
    r <- exp(loglik(q[i-1]+sign) - x[i-1])
    if(runif(1) < r){
      q[i] <- q[i-1]+sign
      x[i] <- loglik(q[i-1]+sign)
    }else{
      q[i] <- q[i-1]
      x[i] <- x[i-1]
    }
  }
  data.frame(init.q,step=1:step,q,loglik=x)
}
step.n <-100
q1 <- metropolis(init.q = 0.3,alpha=.01,step = step.n)
# p.179
qplot(data=q1,x=step,y=q,col=factor(init.q))+
  geom_line() +
  ggtitle(paste0("mcmc step:",step.n))
qplot(data=q1,x=q,binwidth=0.01)+ggtitle(paste0("mcmc step:",step.n))

step.n <-1000
q2 <-metropolis(init.q = 0.3,alpha=.01,step = step.n)
qplot(data=q2,x=step,y=q,col=factor(init.q))+
  geom_line() +
  ggtitle(paste0("mcmc step:",step.n))
qplot(data=q2,x=q,binwidth=0.01)+ggtitle(paste0("mcmc step:",step.n))

step.n <-100000
q3 <-metropolis(init.q = 0.3,alpha=.01,step = step.n)
qplot(data=q3,x=step,y=q,col=factor(init.q))+
  geom_line() +
  ggtitle(paste0("mcmc step:",step.n))
qplot(data=q3,x=q,binwidth=0.01)+ggtitle(paste0("mcmc step:",step.n))
summary(q3$q)

# p.181
step.n <-500
q1 <- metropolis(init.q = 0.3,step = step.n)
q2 <- metropolis(init.q = 0.4,step = step.n)
q3 <- metropolis(init.q = 0.5,step = step.n)
q4 <- metropolis(init.q = 0.6,step = step.n)
qplot(data=rbind(q1,q2,q3,q4),x=step,y=q,col=factor(init.q),alpha=0.5)+geom_line()