0

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

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

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

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$ 值的抽奖的总体平均值以及每个值的所有抽奖的总体平均值（在这种情况下应该是 0，但是在其他情况下，当数据以结构方式随时间变化时可能不会）。

我对rstan（使用过rbugsand rjags）比较陌生，所以我可能根本不知道有一些简单的功能可以使这个过程更容易。

提前感谢您的帮助。

4

1 回答 1

0

您可以使用更少的代码来重现您的第二个图形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)
于 2018-02-19T20:46:59.423 回答