StochMCMC.jl¶
Author: | Al-Ahmadgaid B. Asaad (alasaadstat@gmail.com | http://alstatr.blogspot.com/) |
---|---|
Requires: | julia releases 0.4.1 or later |
Date: | May 12, 2017 |
License: | MIT |
Website: | https://github.com/alstat/StochMCMC.jl |
(Under heavy construction, Target Finish Date on Monday 12pm - Philippine Time)
A julia package for Stochastic Gradient Markov Chain Monte Carlo. The package is part of my master’s thesis entitled Bayesian Autoregressive Distributed Lag via Stochastic Gradient Hamiltonian Monte Carlo or BADL-SGHMC, under the supervision of Dr. Joselito C. Magadia of School of Statistics, University of the Philippines Diliman. This work aims to accommodate other Stochastic Gradient MCMCs in the near future.
Installation¶
To install the package, run the following
Pkg.clone("https://github.com/alstat/StochMCMC.jl")
And to load the package, run
using StochMCMC
Hamiltonian Monte Carlo¶
Setup the necessary paramters including the gradients. The potential energy is the negative logposterior given by U
, the gradient is dU
; the kinetic energy is the standard Gaussian function given by K
, with gradient dK
. Thus,
U <- function(theta) - logpost(theta)
K <- function(p, Sigma = diag(length(p))) (t(p) %*% solve(Sigma) %*% p) / 2
dU <- function(theta, alpha = a, b = eye_mat[1, 1]) {
c(
- alpha * sum(y - (theta[1] + theta[2] * x)),
- alpha * sum((y - (theta[1] + theta[2] * x)) * x)
) + b * theta
}
dK <- function (p, Sigma = diag(length(p))) solve(Sigma) %*% p
Run the MCMC:
set.seed(123)
HMC_object <- HMC(U, K, dU, dK, c(0, 0), 2)
chain2 <- mcmc(HMC_object, leapfrog_params = c(eps = .09, tau = 20), r = 10000)
Extract the estimate
est2 <- colMeans(chain2[seq((burn_in + 1), nrow(chain2), by = thinning), ])
est2
# [1] -0.2977521 -0.5158439
Stochastic Gradient Hamiltonian Monte Carlo¶
Define the gradient noise and other parameters of the SGHMC:
dU_noise <- function(theta, alpha = a, b = eye_mat[1, 1]) {
c(
- alpha * sum(y - (theta[1] + theta[2] * x)),
- alpha * sum((y - (theta[1] + theta[2] * x)) * x)
) + b * theta + matrix(rnorm(2), 2, 1)
}
Run the MCMC:
set.seed(123)
SGHMC_object <- SGHMC(dU_noise, dK, diag(2), diag(2), diag(2), init_est = c(0, 0), 2)
chain3 <- mcmc(SGHMC_object, leapfrog_params = c(eps = .09, tau = 20), r = 10000)
Extract the estimate:
est3 <- colMeans(chain3[seq((burn_in + 1), nrow(chain3), by = thinning), ])
est3
# [1] -0.2920243 -0.4729136
Plot it
p0 <- xyplot(y ~ x, type = c("p", "g"), col = "black") %>%
update(xlab = "x", ylab = "y")
p1 <- histogram(chain3[, 1], col = "gray50", border = "white") %>%
update(xlab = expression(paste("Chain Values of ", w[0]))) %>%
update(panel = function (x, ...) {
panel.grid(-1, -1)
panel.histogram(x, ...)
panel.abline(v = w0, lty = 2, col = "black", lwd = 2)
})
p2 <- histogram(chain3[, 2], col = "gray50", border = "white") %>%
update(xlab = expression(paste("Chain Values of ", w[1]))) %>%
update(panel = function (x, ...) {
panel.grid(-1, -1)
panel.histogram(x, ...)
panel.abline(v = w1, lty = 2, col = "black", lwd = 2)
})
p3 <- xyplot(chain3[, 1] ~ 1:nrow(chain3[, ]), type = c("g", "l"), col = "gray50", lwd = 1) %>%
update(xlab = "Iterations", ylab = expression(paste("Chain Values of ", w[0]))) %>%
update(panel = function (x, y, ...) {
panel.xyplot(x, y, ...)
panel.abline(h = w0, col = "black", lty = 2, lwd = 2)
})
p4 <- xyplot(chain3[, 2] ~ 1:nrow(chain3[,]), type = c("g", "l"), col = "gray50", lwd = 1) %>%
update(xlab = "Iterations", ylab = expression(paste("Chain Values of ", w[1]))) %>%
update(panel = function (x, y, ...) {
panel.xyplot(x, y, ...)
panel.abline(h = w1, col = "black", lty = 2, lwd = 2)
})
p5 <- xyplot(chain3[, 2] ~ chain3[, 1]) %>%
update(type = c("p", "g"), pch = 21, fill = 'white', col = "black") %>%
update(xlab = expression(paste("Chain Values of ", w[0]))) %>%
update(ylab = expression(paste("Chain Values of ", w[1]))) %>%
update(panel = function (x, y, ...) {
panel.xyplot(x, y, ...)
})
p6 <- xyplot(y ~ x, col = "black", fill = "gray80", cex = 1.3, type = "p", pch = 21) %>%
update(xlim = c(-1.1, 1.1), ylim = c(-1.1, 1.1), panel = function(x, y, ...) {
panel.grid(h = -1, v = -1)
xseq <- seq(-1, 1, length.out = 100)
for (i in seq((burn_in + 1), nrow(chain3), by = thinning)) {
yhat <- chain3[i, 1] + chain3[i, 2] * xseq
panel.xyplot(xseq, yhat, type = "l", col = "gray50")
}
panel.xyplot(x, y, ...)
panel.xyplot(xseq, est3[1] + est3[2] * xseq, type = "l", col = "black", lwd = 2)
})
acf1 <- acf(chain3[seq((burn_in + 1), nrow(chain3), by = thinning), 1], plot = FALSE)
acf2 <- acf(chain3[seq((burn_in + 1), nrow(chain3), by = thinning), 2], plot = FALSE)
p7 <- xyplot(acf1$acf ~ acf1$lag, type = c("h", "g"), lwd = 2, col = "black") %>%
update(xlab = "Lags", ylab = expression(paste("Autocorrelations of ", w[1])))
p8 <- xyplot(acf2$acf ~ acf2$lag, type = c("h", "g"), lwd = 2, col = "black") %>%
update(xlab = "Lags", ylab = expression(paste("Autocorrelations of ", w[1])))
grid.arrange(p0, p1, p2, p3, p4, p5, p6, p7, p8, ncol = 3)