1次元ハイブリッドモンテカルロを実装してみた

1次元HMCのコードを書いたのでメモ

動機はNUTSに挑戦する前にパラメータと尤度、ハミルトニアンの振る舞いについて理解しておきたかったので

偉大なる先輩のblogはこちら
ハイブリッドモンテカルロを実装してみた
http://d.hatena.ne.jp/koh_ta/20121116/1353074113

########################################################
#HMC
#initial value
epsilon<-0.25
N<-1000
L<-4

p<-rep(NA, N)
q<-rep(NA, N)
s<-rep(NA, N)
h<-rep(NA, N)

p[1]<-rnorm(0, 1, n=1)
q[1]<-3
h[1]<-hamiltonian(q[1], p[1])

SD<-2
MU<-1
#q obey to N(3, 2.2) f:function itself, df:differentiation function
f<-function(x){dnorm(x, mean=MU, sd=SD)}
df<-function(x){dnorm(x, mean=MU, sd=SD) * (-1/(2*SD^2)*(2*x-2*MU) )}
hamiltonian<-function(q,p){log(exp(f(q))-(p^2)/2)}
#ploting function and differentiation function
plot(df, -1, 3)
par(new=T)
plot(f, -1, 3, col="red")


leapfrog<-function(q_var, p_var, epsilon){
  p_var<-p_var-(epsilon*df(q_var)/2)
  q_var<-q_var+epsilon*p_temp
  p_var<-p_var-(epsilon*df(q_var)/2)
  cat(q_var, p_var, "\r\n")
  c(q_var, p_var)
}

for(i in 2:N){
 p_temp<-p[i-1]
 q_temp<-q[i-1]
 
 for(j in 1:L){
   qp_temp<-leapfrog(q_temp, p_temp, epsilon)
   q_temp<-qp_temp[1]
   p_temp<-qp_temp[2]
 }
 #accept by hamiltonian ratio
 alpha<-min(1, hamiltonian(q_temp, p_temp)/hamiltonian(q[1], p[1]), na.rm=TRUE)
 if(alpha>runif(1)){
   q[i]<-q_temp
   p[i]<- (-1)*p_temp
   s[i]<-1
 }else{
   q[i]<-q[i-1]
   p[i]<-p[i-1]
   s[i]<-0
 }
 h[i]<-hamiltonian(q[i], p[i])
}

ylim<-c(-2, 4)

plot(p,type="l",ylim=ylim)
par(new=T)
plot(q, col="red",type="l",ylim=ylim)
par(new=T)
plot(s, col="blue", ylim=ylim)
par(new=T)
plot(h, col="green",type="l",ylim=ylim)