Bayesian data analysis demo 11.5

library(posterior)
library(bayesplot)
ggplot2::theme_set(bayesplot::theme_default(base_family = "sans"))

Use of posterior package for convergence diagnostics

In the Metropolis assignment it is likely that you have used a for loop for the Metropolis algorithm. Here the Metropolis algorithm has been replaced with independent draws from a normal distribution.

First illustration with a single chain

S <- 1000 # number of draws
D <- 2 # number of parameters
# 2-D matrix with dimensions ‘"iteration"’, and ‘"variable"’
thetas <- matrix(nrow=S, ncol=2)
colnames(thetas) <- c("theta1","theta2")
# Fill the matrix
for (s in 1:S) {
  thetas[s,] <- rnorm(n=2, mean=0, sd=1)
}

We can quickly get several useful summaries including Rank-normalized-Rhat and ESS-bulk and ESS-tail

summarise_draws(thetas)

## # A tibble: 2 × 10
##   variable    mean   median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>      <num>    <num> <num> <num> <num> <num> <num>    <num>    <num>
## 1 theta1   -0.0133 -0.00861 0.964 0.940 -1.62  1.59 1.00     1128.     908.
## 2 theta2    0.0122  0.0410  1.01  1.01  -1.69  1.62 0.999     905.     944.

To compute the basic (split) Rhat and basic ESS as presented in the lecture we can use additionl arguments

summarise_draws(thetas, Rhat=rhat_basic, ESS=ess_basic)

## # A tibble: 2 × 3
##   variable  Rhat   ESS
##   <chr>    <num> <num>
## 1 theta1   0.999 1130.
## 2 theta2   0.999  905.

Compute mean and 5%- and 95%-quantiles and corresponding MCSE’s

summarise_draws(thetas, mean, mcse_mean,
                ~quantile(.x, probs = c(0.05, 0.95)),
                ~mcse_quantile(.x, probs = c(0.05, 0.95)))

## # A tibble: 2 × 7
##   variable    mean mcse_mean  `5%` `95%` mcse_q5 mcse_q95
##   <chr>      <num>     <num> <num> <num>   <num>    <num>
## 1 theta1   -0.0133    0.0287 -1.62  1.59  0.0752   0.0795
## 2 theta2    0.0122    0.0335 -1.69  1.62  0.0678   0.0705

summarise_draws() function coerced the R matrix to a posterior object behind the scenes, but we can also convert it explicitlu to draws object that is better aware of number of iterations and chains. For example, CmdStanR returns posterior draws objects. posterior package provides several functions that make manipulating posterior objects easier than plain R matrix and array types.

thetas <- as_draws_df(thetas)

Illustration with several chains

N <- 1000 # number of iterations per chain
M <- 4 # number of chains
D <- 2 # number of parameters
# 3-D array with dimensions ‘"iteration"’, ‘"chain"’, and ‘"variable"’
thetas <- array(dim=c(N,M,D))
# Fill the array
for (m in 1:M) { # loop over chains
  for (n in 1:N) { # loop over iterations
    thetas[n,m,] <- rnorm(n=2, mean=0, sd=1)
  }
}

We can quickly get several useful summaries including Rank-normalized-Rhat and ESS-bulk and ESS-tail

summarise_draws(thetas)

## # A tibble: 2 × 10
##   variable   mean  median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>     <num>   <num> <num> <num> <num> <num> <num>    <num>    <num>
## 1 ...1     0.0200 0.00900 0.989 0.984 -1.62  1.66  1.00    3974.    3914.
## 2 ...2     0.0304 0.0275  0.993 0.974 -1.61  1.72  1.00    4001.    4065.

To compute the basic (split) Rhat and basic ESS as presented in the lecture we can use additionl arguments

summarise_draws(thetas, Rhat=rhat_basic, ESS=ess_basic)

## # A tibble: 2 × 3
##   variable  Rhat   ESS
##   <chr>    <num> <num>
## 1 ...1      1.00 3978.
## 2 ...2      1.00 4003.

Compute mean and 5%- and 95%-quantiles and corresponding MCSE’s

summarise_draws(thetas, mean, mcse_mean,
                ~quantile(.x, probs = c(0.05, 0.95)),
                ~mcse_quantile(.x, probs = c(0.05, 0.95)))

## # A tibble: 2 × 7
##   variable   mean mcse_mean  `5%` `95%` mcse_q5 mcse_q95
##   <chr>     <num>     <num> <num> <num>   <num>    <num>
## 1 ...1     0.0200    0.0157 -1.62  1.66  0.0259   0.0437
## 2 ...2     0.0304    0.0157 -1.61  1.72  0.0411   0.0266

summarise_draws() function coerced the R 3-D array to a posterior object behind the scenes, but we can also convert it explicitlu to draws object that is better aware of number of iterations and chains. For example, CmdStanR returns posterior draws objects. posterior package provides several functions that make manipulating posterior objects easier than plain R matrix and array types.

thetas <- as_draws_df(thetas)
variables(thetas) <- c("theta1","theta2")

bayesplot package provides Some useful plots to analyse your MCMC draws. You don’t need to include all of these in your assignment or project reports, but they are shown as seeing the plo for your own Metropolis algorithm helps to gain better understanding.

Trace plots

mcmc_trace(thetas)

Autocorrelation plots

mcmc_acf(thetas)

## Warning: The `facets` argument of `facet_grid()` is deprecated as of ggplot2 2.2.0.
## ℹ Please use the `rows` argument instead.
## ℹ The deprecated feature was likely used in the bayesplot package.
##   Please report the issue at <https://github.com/stan-dev/bayesplot/issues/>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Density plots with separate density for each chain

mcmc_dens_overlay(thetas)

Density plots using all chains

mcmc_dens(thetas)

Density plots with median and 50% central interval (see ?mcmc_areas for more options)

mcmc_areas(thetas)

Scatter plots and histograms (works also for more than two variables)

mcmc_pairs(thetas)