R-hat 对抗迭代 RStan

Question

我正在尝试生成如下所示的类似图，以显示 R-hat 在迭代过程中的变化：

我尝试了以下选项：

summary(fit1)$summary : 给出 R-hat 所有链都被合并
summary(fit1)$c_summary : 分别为每个链提供 R-hat

你能帮我得到给定参数每次迭代的R-hat吗？

Answer 1

rstan 提供 Rhat() function, which takes a matrix of iterations x chains and returns R-hat. We can extract this matrix from the fitted model and apply Rhat() cumulatively over it. The code below uses the 8 schools model as an example (copied from the getting started guide).

library(tidyverse)
library(purrr)
library(rstan)
theme_set(theme_bw())

# Fit the 8 schools model.
schools_dat <- list(J = 8, 
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))
fit <- stan(file = 'schools.stan', data = schools_dat)

# Extract draws for mu as a matrix; columns are chains and rows are iterations.
mu_draws = as.array(fit)[,,"mu"]

# Get the cumulative R-hat as of each iteration.
mu_rhat = map_dfr(
  1:nrow(mu_draws),
  function(i) {
    return(data.frame(iteration = i,
                      rhat = Rhat(mu_draws[1:i,])))
  }
)

# Plot iteration against R-hat.
mu_rhat %>%
  ggplot(aes(x = iteration, y = rhat)) +
  geom_line() +
  labs(x = "Iteration", y = expression(hat(R)))

R-hat 对抗迭代 RStan

R-hat against iterations RStan

stan

rstan