Skip to content

Instantly share code, notes, and snippets.

@yuu-ito

yuu-ito/midori_study_8.r

Last active August 29, 2015 14:19
Show Gist options
  • Save yuu-ito/367af8b857e88e307a00 to your computer and use it in GitHub Desktop.
Save yuu-ito/367af8b857e88e307a00 to your computer and use it in GitHub Desktop.
Download ZIP
久保みどり本 8章 今回は手を動かすというよりも読み物メイン。
Raw
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment