knitr::opts_chunk$set( collapse = TRUE, comment = "#>", cache = TRUE )
library(gentleHMC) library(tidyverse) library(rstan) library(tidyverse) options(mc.cores = parallel::detectCores()) rstan_options(auto_write = TRUE)
Animate an HMC sampler based on a two-variable banana-shaped distribution
b <- banana(a = 1.25, b = .5, r = .95) int_length <- 1.9 n_steps <- 14 n_samples <- 100 dpi <- 300 set.seed(100) out <- animate_HMC(b, int_length/n_steps, n_steps, n_samples, c(0, 4), c(1, 1)) # gganimate::gganimate(out$gg, filename = here::here('vignettes/part_1_files/ani.mp4'), # saver = 'mp4', title_frame = FALSE, interval = 1/n_steps, # ani.height = 6/16*9, ani.width = 6, # ani.dev = function(...) grDevices::png(..., res = dpi, units = "in", # type = 'cairo', bg = 'white'), # other.opts = '-c:v libx264 -profile:v high -pix_fmt yuv420p -crf 1')
Estimate the same model using Stan
library(rstan) sm <- stan_model(model_code = ' data {} transformed data { real a; real b; real r; a = 1.25; b = .5; r = .95; } parameters { real x; real y; } model { target += (((a^8*b^2 + x^2*(1 + b*x*(2*r + b*x)) - 2*a^6*b*y - 2*a^2*x*(r + b*x)*y + a^4*(2*b*x*(r + b*x) + y^2))/(a^2*(-1 + r^2)) - 2*log(pi()) - log(4 - 4*r^2))/2); } ') fit <- sampling(sm, control = list(adapt_delta = .9)) ssp <- rstan::get_sampler_params(fit, inc_warmup = FALSE) purrr::map(seq_along(ssp), ~dplyr::mutate(as.data.frame(ssp[[.]]), chain = .)) %>% dplyr::bind_rows() %>% tibble::as_tibble() %>% dplyr::summarise_at(dplyr::vars(stepsize__, n_leapfrog__), dplyr::funs(mean)) pd <- fit@sim$samples pddf <- purrr::map(seq_along(pd), ~dplyr::mutate(as.data.frame(pd[[.]]), Chain = .)) %>% dplyr::bind_rows() pddf %>% dplyr::select(-lp__) %>% bayesplot::mcmc_scatter(alpha = .1) + ggplot2::theme_minimal() + ggplot2::labs(x = expression(theta[1]), y = expression(theta[2])) pddf %>% bayesplot::mcmc_acf(pars = c('x', 'y')) pddf %>% bayesplot::mcmc_combo(pars = c('x', 'y')) pddf %>% bayesplot::mcmc_hex(pars = c('x', 'y'))
Animate a single HMC particle circling around a Gaussian distribution
n_steps <- 150 int_length <- 2*pi * (n_steps - 1)/n_steps dpi <- 300 set.seed(100) U <- function(q) {(1/2)*(q[1]^2 + q[2]^2) + log(2) + log(pi)} grad_U <- function(q)q track <- function(p = c(1, .5), q = c(.5, .5)) { out <- gentleHMC:::single_HMC_sample(U, grad_U, int_length/n_steps, n_steps, q, p) td <- bind_cols(out %>% pluck('intermediate_q') %>% map(~as.data.frame(t(.))) %>% bind_rows() %>% rename(x.q = V1, y.q = V2), out %>% pluck('intermediate_p') %>% map(~as.data.frame(t(.))) %>% bind_rows() %>% rename(x.p = V1, y.p = V2)) %>% mutate(z = sqrt(x.p^2 + y.p^2), step = row_number()) %>% as_tibble %>% select(-y.p, -x.p) td } td <- track() bd <- crossing(x.q = seq(-1, 1, by = .01)*2, y.q = seq(-1, 1, by = .01)*2) %>% mutate(z = exp(-.5 * (x.q^2 + y.q^2))) %>% crossing(step = seq_len(n_steps + 1)) g <- ggplot(NULL, aes(x = x.q, y = y.q, frame = step)) + geom_tile(aes(fill = z), bd) + geom_point(aes(size = -z, alpha = -z), td, colour = 'black', fill = 'black') + scale_fill_gradient(low = 'white', high = gray(.25)) + scale_alpha_continuous(range = c(.5, 1)) + theme_minimal() + coord_cartesian(xlim = c(-1, 1)*2, ylim = c(-1, 1) * 9/16 * 2) + guides(alpha = 'none', fill = 'none', colour = 'none', size = 'none') + labs(x = expression(theta[1]), y = expression(theta[2])) # gganimate::gganimate(g, filename = here::here('vignettes/part_1_files/one_particle.mp4'), cumulative = TRUE, # saver = 'mp4', title_frame = FALSE, interval = 1/16, # ani.height = 6/16*9, ani.width = 6, # ani.dev = function(...) grDevices::png(..., res = dpi, units = "in", # type = 'cairo', bg = 'white'), # other.opts = '-c:v libx264 -profile:v high -pix_fmt yuv420p -crf 1')
Create a sequence of animations showing HMC sampling from a bivariate Gaussian distribution
f <- function(t, q, p) { int_length <- pi/3 n_steps <- 20 n_samples <- 1 dpi <- 300 b <- list(U = function(q) {(1/2)*(q[1]^2 + q[2]^2) + log(2) + log(pi)}, grad_U = function(q)q) out <- gentleHMC:::animate_HMC(b, int_length/n_steps, n_steps, n_samples, q, p) gg <- out$gg + coord_cartesian(xlim = c(-1, 1), ylim = c(-1, 1)) # gganimate::gganimate(gg, filename = here::here(str_glue('hmc_t_{t}.mp4')), # saver = 'mp4', title_frame = FALSE, interval = 1/n_steps, # ani.height = 6/16*9, ani.width = 6, # ani.dev = function(...) grDevices::png(..., res = dpi, units = "in", # type = 'cairo', bg = 'white'), # other.opts = '-c:v libx264 -profile:v high -pix_fmt yuv420p -crf 1') } f(1, c(0, 0), c(.5, -.5)) f(2, c(0.433, -0.433), c(-0.00997, 0.616)) f(3, c(0.208, 0.317), c(.5, .15)) f(4, c(0.537, 0.288), c(-.5, .15))
Create an animation depicting inefficient HMC sampling from a bivariate Gaussian distribution
int_length <- .99 * 2 * pi n_steps <- 20 n_samples <- 10 dpi <- 300 b <- list(U = function(q) {(1/2)*(q[1]^2 + q[2]^2) + log(2) + log(pi)}, grad_U = function(q)q) out <- animate_HMC(b, int_length/n_steps, n_steps, n_samples, c(.5, .5), c(-.5, .25)) gg <- out$gg + coord_cartesian(xlim = c(-2, 2), ylim = c(-2, 2)) # gganimate::gganimate(gg, filename = here::here('vignettes/part_1_files/badL.mp4'), # saver = 'mp4', title_frame = FALSE, interval = 1/n_steps, # ani.height = 6/16*9, ani.width = 6, # ani.dev = function(...) grDevices::png(..., res = dpi, units = "in", # type = 'cairo', bg = 'white'), # other.opts = '-c:v libx264 -profile:v high -pix_fmt yuv420p -crf 1')
Generate an animation showing a side-by-side comparision of HMC and RW M-H
b <- banana(a = 1.25, b = .5, r = .99) int_length <- 1.9 n_steps <- 14 n_samples <- 50 dpi <- 300 set.seed(100) out <- gentleHMC:::animate_HMC(b, .5 * int_length/n_steps, n_steps, n_samples, start_q = c(1, 1), start_p = NULL) hmc <- out$data %>% filter(subframe > 1) set.seed(100) out <- gentleHMC:::animate_HMC(b, 1.0225 * int_length/n_steps, 1, n_steps * n_samples, start_q = c(1, 1), start_p = NULL) rw <- out$data %>% filter(subframe > 1) pd <- bind_rows(hmc %>% mutate(type = 'Hamiltonian'), rw %>% mutate(type = 'Metropolis')) %>% group_by(type) %>% mutate(frame = dense_rank(frame)) %>% ungroup gg <- ggplot2::ggplot(pd, ggplot2::aes(x = x, y = y, frame = frame)) + ggplot2::geom_point(ggplot2::aes(shape = accept_path), data = pd %>% dplyr::filter(accept_point | reject_point), show.legend = FALSE, fill = '#eb3454', colour = gray(.5)) + ggplot2::geom_path(ggplot2::aes(group = factor(draw):factor(frame), colour = cl), show.legend = FALSE) + ggplot2::labs(x = expression(theta[1]), y = expression(theta[2])) + viridis::scale_colour_viridis(alpha = 1, begin = 0, end = 1, direction = 1, option = 'D', discrete = FALSE) + ggplot2::facet_wrap(~type) + ggplot2::theme_minimal() + ggplot2::geom_blank(ggplot2::aes(shape = accept_path), data = dplyr::mutate(pd[1:2,], accept_path = c(FALSE, TRUE))) + ggplot2::scale_shape_manual(values = c(4, 21)) + ggplot2::theme(strip.text = element_text(colour = '#eb3454', face = 'bold', size = rel(1.2))) # gganimate::gganimate(gg, filename = here::here('vignettes/part_1_files/hmc.rw.mp4'), # saver = 'mp4', title_frame = FALSE, interval = 2/n_steps, # ani.height = 6/16*9, ani.width = 6, # ani.dev = function(...) grDevices::png(..., res = dpi, units = "in", # type = 'cairo', bg = 'white'), # other.opts = '-c:v libx264 -profile:v high -pix_fmt yuv420p -crf 1')
Generate HCM and RW trace plots (using Stan).
sm <- stan_model(model_code = ' functions { real pc_exponential_rate(real U, real a) {return - log(a)/U;} } data { int<lower = 0> nr; int<lower = 0> nc; matrix[nr, nc] y; } transformed data { real sigma_rate; sigma_rate = pc_exponential_rate(.5, .1); } parameters { vector[nr] theta_; vector<lower = 0>[nr] sigma_; } transformed parameters { vector[nr] theta; vector[nr] sigma; sigma = sigma_ / sigma_rate; theta = theta_ + 10; } model { theta_ ~ normal(0, 1); sigma_ ~ exponential(1); for (i in 1:nr) { y[i,] ~ normal(theta[i], sigma[i] / fabs(theta[i])); } } ') nc <- 30 nr <- 1 sigma <- matrix(rexp(nr) / (- log(.1)/.5), nrow = nr, ncol = nc) theta <- matrix(10 + rnorm(nr), nrow = nr, ncol = nc) y <- matrix(rnorm(nr * nc), nrow = nr) * sigma / abs(theta) + theta fit <- sampling(sm, data = list(y = y, nc = nc, nr = nr), control = list(adapt_delta = .9), chains = 1) fit2 <- sampling(sm, data = list(y = y, nc = nc, nr = nr), algorithm = 'HMC', control = with(list(e = 1/1000), list(adapt_engaged = TRUE, stepsize = e, int_time = 1e-64, metric = 'unit_e')), iter = 20000, thin = 10, chains = 1) ggplot_the_par_df <- function(fit) { get_par_df(fit) %>% mutate(Iteration = row_number()) %>% select(Iteration, `1` = sigma.1., `2` = theta.1.) %>% gather(par, val, -Iteration) %>% ggplot(aes(x = Iteration, y = val)) + geom_line(colour = '#eb3454') + facet_wrap(~par, scale = 'free_y', ncol = 1, labeller = label_bquote(cols = theta[.(par)])) + labs(y = NULL) } ggplot_the_par_df(fit) # ggsave(here::here('vignettes/part_1_files/trace_good.png'), width = 6, height = 3.5) ggplot_the_par_df(fit2) # ggsave(here::here('vignettes/part_1_files/trace_bad.png'), width = 6, height = 3.5)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.