我试图了解我哪里出错了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
(使用过rbugs
and rjags
)比较陌生,所以我可能根本不知道有一些简单的功能可以使这个过程更容易。
提前感谢您的帮助。