Static Hamiltonian Monte Carlo

ggplot2 is used for plotting, tidyr for manipulating data frames

library(ggplot2)
theme_set(theme_minimal())
library(tidyr)
library(ggforce)
library(MASS)
library(rprojroot)
library(rstan)
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.99
S <- diag(2)
S[1, 2] <- r
S[2, 1] <- r

Draw samples from the toy distribution to visualize 90% HPD interval with ggplot’s stat_ellipse()

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

see BDA3 p. 85 for how to compute HPD for multivariate normal in 2d-case contour for 90% HPD is an ellipse, whose semimajor axes can be computed from the eigenvalues of the covariance matrix scaled by a value selected to get ellipse match the density at the edge of 90% HPD. Angle of the ellipse could be computed from the eigenvectors, but since the marginals are same we know that angle is pi/4 Starting value of the chain

t1 <- -2.5
t2 <- 2.5

Number of iterations.

M <- 5000

Insert your own HMC sampling here

# Allocate memory for the samples
tt <- matrix(rep(0, 2*M), ncol = 2)
tt[1,] <- c(t1, t2)    # Save starting point
# For demonstration load pre-computed values
# Replace this with your algorithm!
# tt is a M x 2 array, with M samples of both theta_1 and theta_2
load(root("demos_ch11","demo12_1b.RData"))
# Here we have intentionally used a very small step size for smooth
# simulations, but for more efficient simulations larger step size
# could be used

The rest is for illustration Creat data frame

df <- data.frame(id=rep(1,4000),
                 iter=rep(1:100, each=40), 
                 th1 = tt[1:4000, 1],
                 th2 = tt[1:4000, 2],
                 th1l = c(tt[1, 1], tt[1:(4000-1), 1]),
                 th2l = c(tt[1, 2], tt[1:(4000-1), 2]))

Take the first 1000 draws after warmup of 1

dfs <- data.frame(th1 = tt[seq(41,40001,by=40), 1], th2 = tt[seq(41,40001,by=40), 2])

Base for the plot

# Labels and frame indices for the plot
labs1 <- c('Samples', 'Steps of the sampler', '90% HPD')
p0 <- ggplot() +
  stat_ellipse(data = dft, aes(x = X1, y = X2, color = '3'), level = 0.9) +
  coord_cartesian(xlim = c(-3, 3), ylim = c(-3, 3)) +
  labs(x = 'theta1', y = 'theta2') +
    scale_color_manual(values = c('red', 'forestgreen','blue'), labels = labs1) +
  guides(color = guide_legend(override.aes = list(
    shape = c(16, NA, NA), linetype = c(0, 1, 1)))) +
    theme(legend.position = 'bottom', legend.title = element_blank())

Plot several iterations

for (ind in seq(40,400,by=40)) {
    pp <- p0 +   geom_point(data = df[(ind-40+1):ind,], 
                      aes(th1, th2, color ='1'), alpha=0.3, size=1) +
        geom_segment(data = df[(ind-40+1):ind,], aes(x = th1, xend = th1l, color = '2',
                                                     y = th2, yend = th2l),
                     alpha=0.5) +
        geom_point(data = df[seq(1,ind+1,by=40),], 
                   aes(th1, th2, color ='1'), size=2)
    print(pp)
}

Convergence diagnostics

samp <- tt[seq(41,40001,by=40),]
dim(samp) <- c(dim(samp),1)
samp <- aperm(samp, c(1, 3, 2))
res<-monitor(samp, probs = c(0.25, 0.5, 0.75), digits_summary = 2)
## Inference for the input samples (1 chains: each with iter = 1000; warmup = 500):
## 
##      Q5  Q50 Q95 Mean  SD  Rhat Bulk_ESS Tail_ESS
## V1 -1.5 -0.1 1.4 -0.1 0.9     1      313      322
## V2 -1.6  0.0 1.5 -0.1 0.9     1      313      322
## 
## For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
## effective sample size for bulk and tail quantities respectively (an ESS > 100 
## per chain is considered good), and Rhat is the potential scale reduction 
## factor on rank normalized split chains (at convergence, Rhat <= 1.05).
neff <- res[,'n_eff']
# both theta have owen neff, but for plotting these are so close to each
# other, so that single relative efficiency value is used
s<-dim(samp)[1]
reff <- mean(neff/(s/2))

Visual convergence diagnostics

Collapse the data frame with row numbers augmented into key-value pairs for visualizing the chains

dfb <- dfs
sb <- 1000
dfch <- within(dfb, iter <- 1:sb) %>% gather(grp, value, -iter)

Another data frame for visualizing the estimate of the autocorrelation function

nlags <- 20
dfa <- sapply(dfb, function(x) acf(x, lag.max = nlags, plot = F)$acf) %>%
  data.frame(iter = 0:(nlags)) %>% gather(grp, value, -iter)

A third data frame to visualize the cumulative averages and the 95% intervals

dfca <- (cumsum(dfb) / (1:sb)) %>%
  within({iter <- 1:sb
          uppi <-  1.96/sqrt(1:sb)
          upp <- 1.96/(sqrt(1:sb*reff))}) %>%
  gather(grp, value, -iter)

Visualize the chains

ggplot(data = dfch) +
  geom_line(aes(iter, value, color = grp)) +
  labs(title = 'Trends') +
  scale_color_discrete(labels = c('theta1','theta2')) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Visualize the estimate of the autocorrelation function

ggplot(data = dfa) +
  geom_line(aes(iter, value, color = grp)) +
  geom_hline(aes(yintercept = 0)) +
  labs(title = 'Autocorrelation function') +
  scale_color_discrete(labels = c('theta1', 'theta2')) +
  theme(legend.position = 'bottom', legend.title = element_blank())

Visualize the estimate of the Monte Carlo error estimates

# labels
labs3 <- c('theta1', 'theta2',
           '95% interval for MCMC error',
           '95% interval for independent MC')
ggplot() +
  geom_line(data = dfca, aes(iter, value, color = grp, linetype = grp)) +
  geom_line(aes(1:sb, -1.96/sqrt(1:sb*reff)), linetype = 2) +
  geom_line(aes(1:sb, -1.96/sqrt(1:sb)), linetype = 3) +
  geom_hline(aes(yintercept = 0)) +
  coord_cartesian(ylim = c(-1.5, 1.5), xlim = c(0,1000)) +
  labs(title = 'Cumulative averages') +
  scale_color_manual(values = c('red','blue',rep('black', 2)), labels = labs3) +
  scale_linetype_manual(values = c(1, 1, 2, 3), labels = labs3) +
  theme(legend.position = 'bottom', legend.title = element_blank())

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDEyLjEiCmF1dGhvcjogIkFraSBWZWh0YXJpLCBNYXJrdXMgUGFhc2luaWVtaSIKZGF0ZTogImByIGZvcm1hdChTeXMuRGF0ZSgpKWAiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdGhlbWU6IHJlYWRhYmxlCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlCi0tLQojIyBTdGF0aWMgSGFtaWx0b25pYW4gTW9udGUgQ2FybG8KCmdncGxvdDIgaXMgdXNlZCBmb3IgcGxvdHRpbmcsIHRpZHlyIGZvciBtYW5pcHVsYXRpbmcgZGF0YSBmcmFtZXMKCmBgYHtyIHNldHVwLCBtZXNzYWdlPUZBTFNFLCBlcnJvcj1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeShnZ3Bsb3QyKQp0aGVtZV9zZXQodGhlbWVfbWluaW1hbCgpKQpsaWJyYXJ5KHRpZHlyKQpsaWJyYXJ5KGdnZm9yY2UpCmxpYnJhcnkoTUFTUykKbGlicmFyeShycHJvanJvb3QpCmxpYnJhcnkocnN0YW4pCnJvb3Q8LWhhc19maWxlKCIuQkRBX1JfZGVtb3Nfcm9vdCIpJG1ha2VfZml4X2ZpbGUoKQpgYGAKClBhcmFtZXRlcnMgb2YgYSBub3JtYWwgZGlzdHJpYnV0aW9uIHVzZWQgYXMgYSB0b3kgdGFyZ2V0IGRpc3RyaWJ1dGlvbgoKYGBge3IgfQp5MSA8LSAwCnkyIDwtIDAKciA8LSAwLjk5ClMgPC0gZGlhZygyKQpTWzEsIDJdIDwtIHIKU1syLCAxXSA8LSByCmBgYAoKRHJhdyBzYW1wbGVzIGZyb20gdGhlIHRveSBkaXN0cmlidXRpb24gdG8gdmlzdWFsaXplIDkwJSBIUEQKaW50ZXJ2YWwgd2l0aCBnZ3Bsb3QncyBzdGF0X2VsbGlwc2UoKQoKYGBge3IgfQpkZnQgPC0gZGF0YS5mcmFtZShtdnJub3JtKDEwMDAwMCwgYygwLCAwKSwgUykpCmBgYAoKc2VlIEJEQTMgcC4gODUgZm9yIGhvdyB0byBjb21wdXRlIEhQRCBmb3IgbXVsdGl2YXJpYXRlIG5vcm1hbAppbiAyZC1jYXNlIGNvbnRvdXIgZm9yIDkwJSBIUEQgaXMgYW4gZWxsaXBzZSwgd2hvc2Ugc2VtaW1ham9yCmF4ZXMgY2FuIGJlIGNvbXB1dGVkIGZyb20gdGhlIGVpZ2VudmFsdWVzIG9mIHRoZSBjb3ZhcmlhbmNlCm1hdHJpeCBzY2FsZWQgYnkgYSB2YWx1ZSBzZWxlY3RlZCB0byBnZXQgZWxsaXBzZSBtYXRjaCB0aGUKZGVuc2l0eSBhdCB0aGUgZWRnZSBvZiA5MCUgSFBELiBBbmdsZSBvZiB0aGUgZWxsaXBzZSBjb3VsZCBiZQpjb21wdXRlZCBmcm9tIHRoZSBlaWdlbnZlY3RvcnMsIGJ1dCBzaW5jZSB0aGUgbWFyZ2luYWxzIGFyZSBzYW1lCndlIGtub3cgdGhhdCBhbmdsZSBpcyBwaS80ClN0YXJ0aW5nIHZhbHVlIG9mIHRoZSBjaGFpbgoKYGBge3IgfQp0MSA8LSAtMi41CnQyIDwtIDIuNQpgYGAKCk51bWJlciBvZiBpdGVyYXRpb25zLgoKYGBge3IgfQpNIDwtIDUwMDAKYGBgCgpJbnNlcnQgeW91ciBvd24gSE1DIHNhbXBsaW5nIGhlcmUKCmBgYHtyIH0KIyBBbGxvY2F0ZSBtZW1vcnkgZm9yIHRoZSBzYW1wbGVzCnR0IDwtIG1hdHJpeChyZXAoMCwgMipNKSwgbmNvbCA9IDIpCnR0WzEsXSA8LSBjKHQxLCB0MikgICAgIyBTYXZlIHN0YXJ0aW5nIHBvaW50CiMgRm9yIGRlbW9uc3RyYXRpb24gbG9hZCBwcmUtY29tcHV0ZWQgdmFsdWVzCiMgUmVwbGFjZSB0aGlzIHdpdGggeW91ciBhbGdvcml0aG0hCiMgdHQgaXMgYSBNIHggMiBhcnJheSwgd2l0aCBNIHNhbXBsZXMgb2YgYm90aCB0aGV0YV8xIGFuZCB0aGV0YV8yCmxvYWQocm9vdCgiZGVtb3NfY2gxMSIsImRlbW8xMl8xYi5SRGF0YSIpKQojIEhlcmUgd2UgaGF2ZSBpbnRlbnRpb25hbGx5IHVzZWQgYSB2ZXJ5IHNtYWxsIHN0ZXAgc2l6ZSBmb3Igc21vb3RoCiMgc2ltdWxhdGlvbnMsIGJ1dCBmb3IgbW9yZSBlZmZpY2llbnQgc2ltdWxhdGlvbnMgbGFyZ2VyIHN0ZXAgc2l6ZQojIGNvdWxkIGJlIHVzZWQKYGBgCgpUaGUgcmVzdCBpcyBmb3IgaWxsdXN0cmF0aW9uCkNyZWF0IGRhdGEgZnJhbWUKCmBgYHtyIH0KZGYgPC0gZGF0YS5mcmFtZShpZD1yZXAoMSw0MDAwKSwKICAgICAgICAgICAgICAgICBpdGVyPXJlcCgxOjEwMCwgZWFjaD00MCksIAogICAgICAgICAgICAgICAgIHRoMSA9IHR0WzE6NDAwMCwgMV0sCiAgICAgICAgICAgICAgICAgdGgyID0gdHRbMTo0MDAwLCAyXSwKICAgICAgICAgICAgICAgICB0aDFsID0gYyh0dFsxLCAxXSwgdHRbMTooNDAwMC0xKSwgMV0pLAogICAgICAgICAgICAgICAgIHRoMmwgPSBjKHR0WzEsIDJdLCB0dFsxOig0MDAwLTEpLCAyXSkpCmBgYAoKVGFrZSB0aGUgZmlyc3QgMTAwMCBkcmF3cyBhZnRlciB3YXJtdXAgb2YgMQoKYGBge3IgfQpkZnMgPC0gZGF0YS5mcmFtZSh0aDEgPSB0dFtzZXEoNDEsNDAwMDEsYnk9NDApLCAxXSwgdGgyID0gdHRbc2VxKDQxLDQwMDAxLGJ5PTQwKSwgMl0pCmBgYAoKQmFzZSBmb3IgdGhlIHBsb3QKCmBgYHtyIH0KIyBMYWJlbHMgYW5kIGZyYW1lIGluZGljZXMgZm9yIHRoZSBwbG90CmxhYnMxIDwtIGMoJ1NhbXBsZXMnLCAnU3RlcHMgb2YgdGhlIHNhbXBsZXInLCAnOTAlIEhQRCcpCnAwIDwtIGdncGxvdCgpICsKICBzdGF0X2VsbGlwc2UoZGF0YSA9IGRmdCwgYWVzKHggPSBYMSwgeSA9IFgyLCBjb2xvciA9ICczJyksIGxldmVsID0gMC45KSArCiAgY29vcmRfY2FydGVzaWFuKHhsaW0gPSBjKC0zLCAzKSwgeWxpbSA9IGMoLTMsIDMpKSArCiAgbGFicyh4ID0gJ3RoZXRhMScsIHkgPSAndGhldGEyJykgKwogICAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGMoJ3JlZCcsICdmb3Jlc3RncmVlbicsJ2JsdWUnKSwgbGFiZWxzID0gbGFiczEpICsKICBndWlkZXMoY29sb3IgPSBndWlkZV9sZWdlbmQob3ZlcnJpZGUuYWVzID0gbGlzdCgKICAgIHNoYXBlID0gYygxNiwgTkEsIE5BKSwgbGluZXR5cGUgPSBjKDAsIDEsIDEpKSkpICsKICAgIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICdib3R0b20nLCBsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCkpCmBgYAoKUGxvdCBzZXZlcmFsIGl0ZXJhdGlvbnMKCmBgYHtyIH0KZm9yIChpbmQgaW4gc2VxKDQwLDQwMCxieT00MCkpIHsKICAgIHBwIDwtIHAwICsgICBnZW9tX3BvaW50KGRhdGEgPSBkZlsoaW5kLTQwKzEpOmluZCxdLCAKICAgICAgICAgICAgICAgICAgICAgIGFlcyh0aDEsIHRoMiwgY29sb3IgPScxJyksIGFscGhhPTAuMywgc2l6ZT0xKSArCiAgICAgICAgZ2VvbV9zZWdtZW50KGRhdGEgPSBkZlsoaW5kLTQwKzEpOmluZCxdLCBhZXMoeCA9IHRoMSwgeGVuZCA9IHRoMWwsIGNvbG9yID0gJzInLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHkgPSB0aDIsIHllbmQgPSB0aDJsKSwKICAgICAgICAgICAgICAgICAgICAgYWxwaGE9MC41KSArCiAgICAgICAgZ2VvbV9wb2ludChkYXRhID0gZGZbc2VxKDEsaW5kKzEsYnk9NDApLF0sIAogICAgICAgICAgICAgICAgICAgYWVzKHRoMSwgdGgyLCBjb2xvciA9JzEnKSwgc2l6ZT0yKQogICAgcHJpbnQocHApCn0KYGBgCgojIyMgQ29udmVyZ2VuY2UgZGlhZ25vc3RpY3MKCmBgYHtyIH0Kc2FtcCA8LSB0dFtzZXEoNDEsNDAwMDEsYnk9NDApLF0KZGltKHNhbXApIDwtIGMoZGltKHNhbXApLDEpCnNhbXAgPC0gYXBlcm0oc2FtcCwgYygxLCAzLCAyKSkKcmVzPC1tb25pdG9yKHNhbXAsIHByb2JzID0gYygwLjI1LCAwLjUsIDAuNzUpLCBkaWdpdHNfc3VtbWFyeSA9IDIpCm5lZmYgPC0gcmVzWywnbl9lZmYnXQojIGJvdGggdGhldGEgaGF2ZSBvd2VuIG5lZmYsIGJ1dCBmb3IgcGxvdHRpbmcgdGhlc2UgYXJlIHNvIGNsb3NlIHRvIGVhY2gKIyBvdGhlciwgc28gdGhhdCBzaW5nbGUgcmVsYXRpdmUgZWZmaWNpZW5jeSB2YWx1ZSBpcyB1c2VkCnM8LWRpbShzYW1wKVsxXQpyZWZmIDwtIG1lYW4obmVmZi8ocy8yKSkKYGBgCgojIyMgVmlzdWFsIGNvbnZlcmdlbmNlIGRpYWdub3N0aWNzCkNvbGxhcHNlIHRoZSBkYXRhIGZyYW1lIHdpdGggcm93IG51bWJlcnMgYXVnbWVudGVkCmludG8ga2V5LXZhbHVlIHBhaXJzIGZvciB2aXN1YWxpemluZyB0aGUgY2hhaW5zCgpgYGB7ciB9CmRmYiA8LSBkZnMKc2IgPC0gMTAwMApkZmNoIDwtIHdpdGhpbihkZmIsIGl0ZXIgPC0gMTpzYikgJT4lIGdhdGhlcihncnAsIHZhbHVlLCAtaXRlcikKYGBgCgpBbm90aGVyIGRhdGEgZnJhbWUgZm9yIHZpc3VhbGl6aW5nIHRoZSBlc3RpbWF0ZSBvZgp0aGUgYXV0b2NvcnJlbGF0aW9uIGZ1bmN0aW9uCgpgYGB7ciB9Cm5sYWdzIDwtIDIwCmRmYSA8LSBzYXBwbHkoZGZiLCBmdW5jdGlvbih4KSBhY2YoeCwgbGFnLm1heCA9IG5sYWdzLCBwbG90ID0gRikkYWNmKSAlPiUKICBkYXRhLmZyYW1lKGl0ZXIgPSAwOihubGFncykpICU+JSBnYXRoZXIoZ3JwLCB2YWx1ZSwgLWl0ZXIpCmBgYAoKQSB0aGlyZCBkYXRhIGZyYW1lIHRvIHZpc3VhbGl6ZSB0aGUgY3VtdWxhdGl2ZSBhdmVyYWdlcwphbmQgdGhlIDk1JSBpbnRlcnZhbHMKCmBgYHtyIH0KZGZjYSA8LSAoY3Vtc3VtKGRmYikgLyAoMTpzYikpICU+JQogIHdpdGhpbih7aXRlciA8LSAxOnNiCiAgICAgICAgICB1cHBpIDwtICAxLjk2L3NxcnQoMTpzYikKICAgICAgICAgIHVwcCA8LSAxLjk2LyhzcXJ0KDE6c2IqcmVmZikpfSkgJT4lCiAgZ2F0aGVyKGdycCwgdmFsdWUsIC1pdGVyKQpgYGAKClZpc3VhbGl6ZSB0aGUgY2hhaW5zCgpgYGB7ciB9CmdncGxvdChkYXRhID0gZGZjaCkgKwogIGdlb21fbGluZShhZXMoaXRlciwgdmFsdWUsIGNvbG9yID0gZ3JwKSkgKwogIGxhYnModGl0bGUgPSAnVHJlbmRzJykgKwogIHNjYWxlX2NvbG9yX2Rpc2NyZXRlKGxhYmVscyA9IGMoJ3RoZXRhMScsJ3RoZXRhMicpKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gJ2JvdHRvbScsIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkKYGBgCgpWaXN1YWxpemUgdGhlIGVzdGltYXRlIG9mIHRoZSBhdXRvY29ycmVsYXRpb24gZnVuY3Rpb24KCmBgYHtyIH0KZ2dwbG90KGRhdGEgPSBkZmEpICsKICBnZW9tX2xpbmUoYWVzKGl0ZXIsIHZhbHVlLCBjb2xvciA9IGdycCkpICsKICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gMCkpICsKICBsYWJzKHRpdGxlID0gJ0F1dG9jb3JyZWxhdGlvbiBmdW5jdGlvbicpICsKICBzY2FsZV9jb2xvcl9kaXNjcmV0ZShsYWJlbHMgPSBjKCd0aGV0YTEnLCAndGhldGEyJykpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAnYm90dG9tJywgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKQpgYGAKClZpc3VhbGl6ZSB0aGUgZXN0aW1hdGUgb2YgdGhlIE1vbnRlIENhcmxvIGVycm9yIGVzdGltYXRlcwoKYGBge3IgfQojIGxhYmVscwpsYWJzMyA8LSBjKCd0aGV0YTEnLCAndGhldGEyJywKICAgICAgICAgICAnOTUlIGludGVydmFsIGZvciBNQ01DIGVycm9yJywKICAgICAgICAgICAnOTUlIGludGVydmFsIGZvciBpbmRlcGVuZGVudCBNQycpCmdncGxvdCgpICsKICBnZW9tX2xpbmUoZGF0YSA9IGRmY2EsIGFlcyhpdGVyLCB2YWx1ZSwgY29sb3IgPSBncnAsIGxpbmV0eXBlID0gZ3JwKSkgKwogIGdlb21fbGluZShhZXMoMTpzYiwgLTEuOTYvc3FydCgxOnNiKnJlZmYpKSwgbGluZXR5cGUgPSAyKSArCiAgZ2VvbV9saW5lKGFlcygxOnNiLCAtMS45Ni9zcXJ0KDE6c2IpKSwgbGluZXR5cGUgPSAzKSArCiAgZ2VvbV9obGluZShhZXMoeWludGVyY2VwdCA9IDApKSArCiAgY29vcmRfY2FydGVzaWFuKHlsaW0gPSBjKC0xLjUsIDEuNSksIHhsaW0gPSBjKDAsMTAwMCkpICsKICBsYWJzKHRpdGxlID0gJ0N1bXVsYXRpdmUgYXZlcmFnZXMnKSArCiAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGMoJ3JlZCcsJ2JsdWUnLHJlcCgnYmxhY2snLCAyKSksIGxhYmVscyA9IGxhYnMzKSArCiAgc2NhbGVfbGluZXR5cGVfbWFudWFsKHZhbHVlcyA9IGMoMSwgMSwgMiwgMyksIGxhYmVscyA9IGxhYnMzKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gJ2JvdHRvbScsIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfYmxhbmsoKSkKYGBgCgo=