Metropolis algorithm + Rhat (PSRF) demonstration

ggplot2 and gridExtra are used for plotting, tidyr for manipulating data frames

``````library(ggplot2)
theme_set(theme_minimal())
library(gridExtra)
library(tidyverse)
library(gganimate)
library(MASS)
library(rstan)
library(rprojroot)
root<-has_file(".BDA_R_demos_root")\$make_fix_file()``````

Parameters of a Normal distribution used as a toy target distribution

``````y1 <- 0
y2 <- 0
r <- 0.8
S <- diag(2)
S[1, 2] <- r
S[2, 1] <- r``````

Sample from the toy distribution to visualize 90% HPD interval with ggplotâ€™s stat_ellipse()

``dft <- data.frame(mvrnorm(100000, c(0, 0), S))``

Since, implementation of the Metropolis algorithm is one of the exercises, we load here pre-computed chains and Rhat-values

The proposal distribution was intentionally selected to be slightly too small, to better illustrate convergence diagonstics

`tts' contains draws,`p1â€™ and â€˜p2â€™ contain Rhat values for t1 and t2 using 50% warm-up, `pp1â€™ and â€˜pp2â€™ contain Rhat values for t1 and t2 using 10% warm-up, Rhat-values have been computed for each time-step.

``load(root("demos_ch11","demo11_4.RData"))``

Transform the first s1 rows of the data into a â€˜tidyâ€™ format for plotting

``````s1 <- 50
dfs1 <- data.frame(iter = 1:s1, tts[1:s1, 1, ]) %>%
gather(chain, th1, -iter) %>%
within({th2 <- c(tts[1:s1, 2, ])         # xl and yl specify the previous
th2l <- c(th2[1], th2[-length(th2)])   # draw in the chain for
th1l <- c(th1[1], th1[-length(th1)])}) # plotting``````

Fix the incorrect lagged values, the lagged value of the first draw in the chain (for plotting) is the value itself (instead of the last value of the previous chain)

``````sind <- 0:9*s1+1
dfs1[sind, c('th1l','th2l')] <- dfs1[sind, c('th1','th2')]``````

Another data frame with all draws

``````inds2 <- 1:10000
dfs2 <- data.frame(iter = inds2, tts[inds2, 1, ]) %>%
gather(chain, theta1, -iter) %>%
within(theta2 <- c(tts[inds2, 2, ])) %>%
gather(var, val, -iter, -chain)``````

Third data frame with Rhat values

``````indsp <- seq(10, length(p1), 10)
dfp <- data.frame(iter = indsp,
theta1 = p1[indsp],
theta2 = p2[indsp]) %>%
gather(var,rhat,-iter) ``````

Construct a 2d-plot of the 50 first iterations of the chains

``````frame = rep(1:s1, 10)
chains1 <- ggplot(data = dfs1) +
geom_segment(aes(x = th1, xend = th1l, y = th2, yend = th2l, color = chain,
group = chain)) +
geom_point(data = dfs1[sind, ], aes(x = th1, y = th2, color = chain)) +
stat_ellipse(data = dft, aes(x = X1, y = X2), level = 0.9, color = 'black') +
coord_cartesian(xlim = c(-4, 4), ylim = c(-4, 4)) +
labs(x = 'theta1', y = 'theta2') +
scale_color_discrete(guide = "none")``````

Animate s1 first iterations of the chains. At some points some of the chains seem to halt for

``#  a moment. What really happens at that point is that``

they draw possibly a few points that are rejected (rejected points not shown) and thus the chain is not moving

``````anim <- animate(chains1 +
transition_reveal(along=iter) +

Show the animation

``anim``