如何在 rstan 中绘制参数时间序列

Question

我有一个已经 rstan 名为“fit”的合适对象，类型为：

Inference for Stan model: MySecondModel.
1 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=1000.

                                 mean se_mean    sd        2.5%         25%         50%         75%       97.5% n_eff Rhat
step_size                        0.01    0.00  0.00        0.01        0.01        0.01        0.01        0.01   493 1.00
theta_init                       0.18    0.00  0.02        0.14        0.17        0.18        0.19        0.22  1167 1.00
theta_steps_unscaled[1]         -0.02    0.02  0.91       -1.81       -0.60        0.00        0.58        1.73  1453 1.00
theta_steps_unscaled[2]         -0.02    0.03  0.94       -1.95       -0.65        0.00        0.57        1.81  1382 1.00
theta_steps_unscaled[3]          0.14    0.03  0.95       -1.73       -0.51        0.15        0.80        2.01  1224 1.00
theta_steps_unscaled[4]         -0.34    0.02  0.91       -2.05       -0.95       -0.36        0.26        1.44  1499 1.00
theta_steps_unscaled[5]         -0.32    0.02  0.88       -2.02       -0.94       -0.31        0.31        1.32  2079 1.00
theta_steps_unscaled[6]         -0.99    0.02  0.86       -2.69       -1.56       -0.98       -0.41        0.66  1759 1.00
theta_steps_unscaled[7]          0.23    0.03  0.94       -1.59       -0.40        0.24        0.88        2.06  1411 1.00
theta_steps_unscaled[8]         -0.34    0.02  0.89       -2.13       -0.94       -0.38        0.24        1.34  1774 1.00
theta_steps_unscaled[9]         -0.78    0.02  0.89       -2.49       -1.39       -0.75       -0.17        0.86  1509 1.00
theta_steps_unscaled[10]        -0.66    0.02  0.87       -2.34       -1.28       -0.71       -0.05        1.05  1928 1.00

我想访问整个向量“theta_steps_unscaled”，首先是最终调试链，最后将其 mean/mode 绘制为时间序列。我一直在寻找类似 regex_pars 的东西，但似乎没有“stan”对象。

我能达到的最好结果是下面那个，但这不是我真正想要的。

我也试过这个，分析链：

traceplot(fit, "theta_steps_unscaled[6]")

Error in mcmc.list(x) : Arguments must be mcmc objects

最后将其转换为 mcmc：

traceplot(as.mcmc(fit, "theta_steps_unscaled[6]"))

虽然这引发了：

Error in xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 'list' object cannot be coerced to type 'double'

请注意我的最后一个目标是绘制整个矢量 theta_steps_unscaled[1:some_index]，获得类似的。即具有可信区间、mode/mean 和其他贝叶斯参数。

Answer 1

tidybayes 是一种以方便的格式提取样本的简便方法。（它也有方便的绘图功能；this vignette gives an overview。）我会使用 tidybayes 来提取您的样本，并使用 ggplot2 以您想要的方式绘制它们。

首先，提取采样值：

library(tidyverse)
theme_set(theme_bw())
library(tidybayes)

samples.df = fit %>%
  spread_draws(theta_steps_unscaled[t]) %>%
  ungroup()

接下来，用 x 轴上的步长和 y 轴上的可信区间绘制它们：

samples.df %>%
  group_by(t) %>%
  summarise(median = median(theta_steps_unscaled),
            upper.95 = quantile(theta_steps_unscaled, 0.975),
            lower.95 = quantile(theta_steps_unscaled, 0.025)) %>%
  ggplot(aes(x = t, y = median)) +
  geom_line(size = 1, color = "blue") +
  geom_ribbon(aes(ymin = lower.95, ymax = upper.95),
              fill = "blue", alpha = 0.1) +
  scale_x_continuous("step", breaks = 1:10) +
  scale_y_continuous(expression(theta))

情节看起来像这样：

如何在 rstan 中绘制参数时间序列

How to plot parameters timeseries in rstan

r

stan

rstan