R語言實現hamiltonian mcmc
MCMC步驟是：
將第一步中的提議分布改為通過哈密頓動力學來尋找建議點，就是HMC。
哈密頓動力學的微分描述是：
由以上方程可以推出哈密頓力學的三個特點：
- reversibility
- conversation of hamiltonian
- volume preservation
接下來根據哈密頓系統尋找下一個建議點：
使用leapfrog方法：
仔細觀察就可以發現這個演算法其實是5.1，5.2兩個微分方程的離散化，因此新的建議點滿足哈密頓系統的特點。
這個方法是Euler『s method一路改過來的，具體如下：
以上根據哈密頓力學找到了建議點，接下來就是計算MCMC步驟中的接受率，因為是q和p的聯合分布：
所以接受率是：
因為符合哈密頓力學，所以H不變，接受率為1.。但是因為離散化處理，所以H不會絕對相等，但是仍然接近1.
以上就是HMC的原理。
主要代碼來自Neal。
##########################################################################
#hamiltonian update
##########################################################################
HMC = function (U, grad_U, epsilon, L, current_q)
{
q = current_q
p = matrix(rnorm(length(q),0,1), ncol = 1) # independent standard normal variates
current_p = p
# Make a half step for momentum at the beginning
p = p - epsilon * grad_U(q) / 2
# Alternate full steps for position and momentum
for (i in 1:L)
{
# Make a full step for the position
q = q + epsilon * p
# Make a full step for the momentum, except at end of trajectory
if (i!=L) p = p - epsilon * grad_U(q)
}
# Make a half step for momentum at the end.
p = p - epsilon * grad_U(q) / 2
# Negate momentum at end of trajectory to make the proposal symmetric
p = -p
# Evaluate potential and kinetic energies at start and end of trajectory
current_U = U(current_q)
current_K = sum(current_p^2) / 2
proposed_U = U(q)
proposed_K = sum(p^2) / 2
# Accept or reject the state at end of trajectory, returning either
# the position at the end of the trajectory or the initial position
if (runif(1) < exp(current_U - proposed_U + current_K - proposed_K))
{
return (q) # accept
}
else
{
return (current_q) # reject
}
}
#########################################################################
#target bivariate normal
#########################################################################
mu = c(0,0)
sigma = matrix(c(1, 0.8, 0.8, 1), nrow = 2)
U_P = function(q)
{
inv_sigma = matrix(c(2.78, -2.22, -2.22, 2.78), nrow = 2)
value = t(q) %*% inv_sigma %*% q/2
return(value)
}
dU = function(q)
{
inv_sigma = matrix(c(2.78, -2.22, -2.22, 2.78), nrow = 2)
K = inv_sigma %*% q
return(K)
}
######################################################################
#simulate
######################################################################
N = 1500
q_matrix = matrix(NA, nrow = 2, ncol = N)
q_init = matrix(c(-1, 1), ncol = 1)
for (i in 1:N)
{
q_matrix[,i] = HMC(U = U_P, grad_U = dU, epsilon = 0.1, L =20,
current_q = q_init)
q_init = q_matrix[,i]
}
plot(q_matrix[1,], q_matrix[2,], type = "b", col= "red")
結果如下：
用animation包演示跳動變化：
library(animation)
saveGIF(
{
for (i in 1:N) {
plot(q_matrix[1,1:i], q_matrix[2,1:i],
xlim = c(-3,3), ylim = c(-3,3),
type = "b", col= "red")
}
}
,movie.name = "hmc.gif")
參考資料：
1.中科大。張偉平，MCMC演算法。
2.handbook of MCMC, 第5章，neal.
推薦閱讀：