Navigate
Courses Tutorials Practice Pricing Tools Upgrade
Lesson 4 of 8

HMC, NUTS and MCMC Diagnostics

Lesson 3 ended on a warning: a trace can look settled while exploring only half the posterior, and 20,000 steps of a crawling walk can be worth a few hundred honest draws. This lesson delivers the two things that warning demands. First, a faster engine: the physics-flavored sampler inside Stan and every serious Bayesian tool, which you will build in base R and run on Asha's order data. Second, a trust dashboard: the checks a professional runs before believing any chain.

By the end of this lesson you will be able to:

  • Measure how much information a correlated chain really holds, by computing effective sample size from the chain itself
  • Explain how Hamiltonian Monte Carlo turns sampling into physics, and why its far-away proposals still get accepted
  • Say exactly what NUTS and warmup automate, and read a divergence warning as evidence rather than noise
  • Catch a chain that looks healthy but is lying, by running several chains from dispersed starts and computing R-hat
  • Run the full pre-flight dashboard (R-hat, effective sample size, divergences, overlaid traces) before trusting any fit

Prerequisites: Lessons 1 to 3 of this course (prior, likelihood, posterior, credible intervals, and the fifteen-line Metropolis sampler you built and tuned), plus base R functions, loops and matrices.

Below is the machine from Lesson 3. Every idea in this lesson is already visible inside it: toggle the step size and watch what the trace and the histogram do. Today you learn to read those two pictures the way a professional does, and to build the engine that makes them look right by design.

Scale it up

Two unknowns, and a walk that seems fine

Lessons 2 and 3 fixed the spread of Asha's order values at $18 by decree, so that her average order value \(\mu\) was the only unknown. Nobody hands you the spread in real life. The honest model for her store treats both as unknown: \(\mu\), the true long-run average order value, and \(\sigma\) (sigma), the true dollar spread of individual orders around that average. The redesigned checkout's first full week gives her 60 orders to learn both from. Each lesson runs in a fresh interactive R session, so build the week right here (a stand-in ledger with the same shape as hers: whole-dollar orders, most between $30 and $90):

RInteractive R
set.seed(6) orders <- round(rnorm(60, mean = 54, sd = 18)) # her first full week: 60 orders head(orders, 12) #> [1] 59 43 70 85 54 61 30 67 55 35 85 33 round(c(n = length(orders), mean = mean(orders), sd = sd(orders)), 1) #> n mean sd #> 60.0 54.6 18.7 

  

  





The Metropolis rule needs nothing new. The log-numerator now takes two arguments, with a prior for each unknown: the average near $55 but loosely (a normal prior with standard deviation 10), and for the spread, which must be positive, an exponential prior that keeps it in the tens of dollars (mean $20). The sampler is the one you built in Lesson 3, except each proposal nudges BOTH unknowns at once:




  

    
RInteractive R

    
logpost <- function(mu, sigma) {
  if (sigma <= 0) return(-Inf)                   # a dollar spread cannot be negative
  sum(dnorm(orders, mu, sigma, log = TRUE)) +    # likelihood of the 60 orders
    dnorm(mu, 55, 10, log = TRUE) +              # prior: average near $55, loosely
    dexp(sigma, rate = 1/20, log = TRUE)         # prior: spread around $20
}

metropolis2 <- function(start, prop_sd, n) {
  chain <- matrix(0, n, 2, dimnames = list(NULL, c("mu", "sigma")))
  chain[1, ] <- start
  acc <- 0
  for (i in 2:n) {
    cand <- chain[i - 1, ] + rnorm(2, 0, prop_sd)     # nudge both unknowns at once
    if (log(runif(1)) < logpost(cand[1], cand[2]) - logpost(chain[i - 1, 1], chain[i - 1, 2])) {
      chain[i, ] <- cand
      acc <- acc + 1
    } else chain[i, ] <- chain[i - 1, ]
  }
  list(chain = chain, accept = acc / (n - 1))
}

set.seed(7)
fit_rw <- metropolis2(start = c(70, 30), prop_sd = c(5, 3.5), n = 20000)
keep <- fit_rw$chain[-(1:2000), ]                # drop 2,000 steps of warm-up
round(c(accept = fit_rw$accept,
        mu = mean(keep[, "mu"]), sigma = mean(keep[, "sigma"])), 2)
#> accept     mu  sigma 
#>   0.28  54.66  18.96 

round(quantile(keep[, "mu"], c(0.025, 0.975)), 1)
#>  2.5% 97.5% 
#>  49.9  59.5 

    

      
      
    

    

  

  





By every Lesson 3 standard this run is a success. Acceptance is 0.28, right in the zone Lesson 3 promised for models with several unknowns. The average order lands at $54.66 with a credible interval of ($49.9, $59.5), and the spread lands near $19, both hugging the data. Look at the trace, and at something new: with two unknowns, each draw is a PAIR, so the posterior is a cloud of plausible (average, spread) combinations:




  

    
RInteractive R

    
par(mfrow = c(1, 2))
plot(keep[1:3000, "mu"], type = "l", col = "navy",
     xlab = "iteration", ylab = "average order (mu)",
     main = "The mu trace: looks settled")
plot(keep[, "mu"], keep[, "sigma"], pch = 16, cex = 0.3,
     col = adjustcolor("navy", 0.25),
     xlab = "average order (mu)", ylab = "order spread (sigma)",
     main = "The joint posterior cloud")
par(mfrow = c(1, 1))

    

      
      
    

    

  

  





A settled trace, sensible numbers, a tidy cloud. Lesson 3 would stamp this and go home. Now look closer, because this chain is quietly wasting almost 90 percent of its work.