# MINUS LOG POSTERIOR DENSITY FUNCTION. Arguments are the parameter vector # and the data matrix. log_pi <- function (par, data) { # Start with log of prior density. d <- dnorm(par[1],0,2,log=TRUE) + dnorm(par[2],0,2,log=TRUE) + dnorm(par[3],0,10,log=TRUE) + sum (dnorm (par[-(1:3)],0,exp(par[2]),log=TRUE)) # Find predicted values for data points. pred <- par[3] + data[,-1,drop=FALSE] %*% par[-(1:3)] e <- data[,1] - pred # Add log likelihood to log prior density. d <- d + sum(dnorm(e,0,exp(par[1]),log=TRUE)) -d } # MINUS GRADIENT OF LOG POSTERIOR. grad_log_pi <- function (par, data) { # Start with gradient of log prior density. g <- c( -par[1]/2^2, -par[2]/2^2 - 2 + sum(par[-(1:3)]^2)*exp(-2*par[2]), -par[3]/10^2, -par[-(1:3)]*exp(-2*par[2]) ) # Find predicted values for data points. pred <- par[3] + data[,-1,drop=FALSE] %*% par[-(1:3)] e <- data[,1] - pred # Add gradient of log likelihood to gradient of log prior density. g[1] <- g[1] + sum (e^2*exp(-2*par[1]) - 1) g[3] <- g[3] + sum (e*exp(-2*par[1])) g[4] <- g[4] + sum (data[,2]*e*exp(-2*par[1])) g[5] <- g[5] + sum (data[,3]*e*exp(-2*par[1])) -g } # CHECK THE GRADIENT COMPUTATION USING FINITE DIFFERENCES. check_grad <- function (U, grad_U, w, epsilon=0.001, ...) { grd <- rep(0,length(w)) for (i in 1:length(w)) { w1 <- w w1[i] <- w1[i]-epsilon w2 <- w w2[i] <- w2[i]+epsilon E1 <- U(w1,...) E2 <- U(w2,...) grd[i] <- (E2-E1) / (2*epsilon) } gr <- grad_U(w,...) list (grd, gr, grd-gr) } # HAMILTONIAN MONTE CARLO UPDATE. HMC = function (U, grad_U, masses, epsilon, L, current_q, ...) { q = current_q p = rnorm(length(q),0,sqrt(masses)) 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 / masses # 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*masses)) proposed_U = U(q,...) proposed_K = sum (p^2 / (2*masses)) # 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)) { accept_count <<- accept_count + 1 return (q) # accept } else { return (current_q) # reject } }