使用 rstan 在 R 中绘制高斯过程答案

【问题标题】：Graphing Gaussian Process in R using rstan使用 rstan 在 R 中绘制高斯过程
【发布时间】：2018-01-18 16:56:45
【问题描述】：

我试图了解rstan 哪里出了问题。我已经找到了一种解决方法，但似乎应该有一个比我想出的更好的从后部绘制图形的选择。

我正在尝试学习如何使用rstan 对与我在CV 上打开的另一个问题相关的高斯过程进行建模（无耻的插头，但如果您有可以帮助的想法，我会全力以赴）。

我想作为第一步，我会尝试通过高斯过程的stan documentation 示例。所以我建立了一个模型，简单地设计来绘制随机平方指数协方差函数。

library(rstan)
library(rstanarm)
library(bayesplot)
library(ggplot2)
options(mc.cores=parallel::detectCores())
rstan_options(auto_write = TRUE)

x<-seq(0, 30, by=.01)

model<-'
data{
    int<lower=1> N;
    real x[N];
  }

transformed data {
  matrix[N, N] L;
  matrix[N, N] K;
  vector[N] mu = rep_vector(0, N);
  for (i in 1:(N - 1)) {
    K[i, i] = 1 + 0.1;
    for (j in (i + 1):N) {
      K[i, j] = exp(-0.5 * square(x[i] - x[j]));
      K[j, i] = K[i, j];
    }
  }
  K[N, N] = 1 + 0.1;
  L = cholesky_decompose(K);
}

parameters {
  vector[N] eta;
}

model {
  eta ~ normal(0, 1);
}
generated quantities {
  vector[N] y;
  y = mu + L*eta;
}
'

我遵循文档的建议，在转换后的数据中加入 Cholesky 分解。

使用stan我拟合模型如下：

dat<-list(N=length(x),
          x=x)

fit <- stan(model_code = model,
            data = dat, 
            iter = 1000, 
            chains = 1, 
            pars = c('y', 'eta'),
            control = list(adapt_delta=.99, 
                           max_treedepth=10)
            )

我可以使用以下代码可视化每个抽奖的后验分布：

posterior<-as.matrix(fit)
mcmc_areas(posterior, 
           pars=c('y[1]', 'y[2]'),
           prob = .90
           )

产生：

我真的很想看看每个过程的结果（不是全部 500 个，而是其中一些随机抽取）。

我尝试了多种替代策略，最终得出以下结论：

post.y<-extract(fit, pars='y')

draws<-sample(1:500, size = 10)

DF<-data.frame(Time=x, y=colMeans(post.y$y), Draw=rep('Mu', length(x)))
for(i in 1:length(draws)){
  DF.temp<-data.frame(Time=x, y=post.y$y[i,], Draw=rep(paste0('posterior', i), length(x)))
  DF<-rbind(DF, DF.temp)
}

g1<-ggplot(aes(x=Time, y=y), data=DF)
g2<-g1+geom_line(aes(x=Time, y=y, group=Draw, color=Draw), data=DF[DF$Draw!='Mu',], alpha=.25, show.legend = F)
g3<-g2+geom_line(aes(x=Time, y=y), data=DF[DF$Draw=='Mu',], lwd=1.5)
g3

这段代码产生：

这似乎需要跳过很多额外的障碍。我尝试了使用rstan 系列中的其他函数（例如ppc_dens_overlay）的替代方法，但它们都导致错误或没有返回我想要的。

所以我在这里的问题实际上是关于替代的、更简单的选项，我可以使用它来可视化每个 $y_i$ 值的抽奖的总体平均值以及每个值的所有抽奖的总体平均值（在这种情况，但在其他情况下，当数据以结构方式随时间变化时可能不会）。

我对@987654337@ 比较陌生（使用过rbugs 和rjags），所以我可能根本不知道有一些简单的函数可以使这个过程更容易。

提前感谢您的帮助。

【问题讨论】：

这是一个来自 stan 开发者的关于高斯过程的教程，其中包括一些图：mc-stan.org/events/stancon2017-notebooks/…。现在似乎有必要在 stan 中绘制高斯过程。
感谢您的链接！我会看看。当然，我关心的不是这样一个简单的模型，而是当我扩展到具有多个协方差函数的更复杂的事物时会发生什么。

标签： r bayesian data-visualization rstan

【解决方案1】：

您可以使用matplot 以更少的代码重现您的第二个图形，它可以方便地处理矩阵数据。

post.y <- rstan::extract(fit, 'y')$y
post.y.sub <- post.y[sample(1:nrow(post.y), 10),]

matplot(x, t(post.y.sub), type = 'l', lty = 1, col = adjustcolor(palette(), 0.25))
lines(colMeans(post.y) ~ x, lwd = 2)

如果您更喜欢 ggplot2，困难的部分是将后验样本放入数据框中。我发现 dplyr 和 tidyr 库在这里很有帮助。它看起来像很多代码，但当您的模型变得更复杂时，它会很灵活。

library(dplyr)
library(tidyr)

df.rep <- post.y %>% 
  t() %>%
  as.data.frame() %>% 
  mutate(x = x) %>% 
  gather(rep, post.y, -x)

df.mean <- df.rep %>% 
  group_by(x) %>% 
  summarize(mu = mean(post.y))

df.rep.sub <- df.rep %>% 
  filter(rep %in% sample(unique(rep), 10))

ggplot() +
  geom_line(data = df.rep.sub, aes(x, post.y, col = rep), alpha = 0.25, show.legend = F) +
  geom_line(data = df.mean, aes(x, mu), lwd = 1.5)

【讨论】：