在没有 运行 单独模型的情况下从先验抽样
Sampling from prior without running a separate model
我想根据这些参数的先验值绘制 stan 模型的参数估计直方图。我尝试通过 运行 在 stan 中构建一个模型，用 ggplot2 绘制它，然后使用 R 的随机生成函数（例如 rnorm()、rbinom()）覆盖先验分布的近似值，但是我 运行 遇到许多缩放问题，这些问题使图表无法正确显示。
我在想一个更好的方法是直接从先验分布中抽样，然后根据参数估计值绘制这些样本，但是 运行 建立一个完整的单独模型 just 从先验中抽样似乎非常耗时。我想知道是否有一种方法可以在现有模型中或与之并行地执行此操作。
这是一个示例脚本。
# simulate linear model
a <- 3 # intercept
b <- 2 # slope
# data
x <- rnorm(28, 0, 1)
eps <- rnorm(28, 0, 2)
y <- a + b*x + eps
# put data into list
data_reg <- list(N = 28, x = x, y = y)
# create the model string
ms <- "
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
}
model {
vector[N] mu;
sigma ~ cauchy(0, 2);
beta ~ normal(0,10);
alpha ~ normal(0,100);
for ( i in 1:N ) {
mu[i] = alpha + beta * x[i];
}
y ~ normal(mu, sigma);
}
"
# now fit the model in stan
fit1 <- stan(model_code = ms, # model string
data = data_reg, # named list of data
chains = 1, # number of Markov chains
warmup = 1e3, # number of warmup iterations per chain
iter = 2e3) # show progress every 'refresh' iterations
# extract the sample estimates
post <- extract(fit1, pars = c("alpha", "beta", "sigma"))
# now for the density plots. Write a plotting function
densFunct <- function (parName) {
g <- ggplot(postDF, aes_string(x = parName)) +
geom_histogram(aes(y=..density..), fill = "white", colour = "black", bins = 50) +
geom_density(fill = "skyblue", alpha = 0.3)
return(g)
}
# plot
gridExtra::grid.arrange(grobs = lapply(names(postDF), function (i) densFunct(i)), ncol = 1)
现在我明白我可以通过简单地从模型字符串中省略可能性来从先验中采样，就像这样
ms <- "
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
}
model {
sigma ~ cauchy(0, 2);
beta ~ normal(0,10);
alpha ~ normal(0,100);
}
"
但是有什么方法可以从第一个模型中的先验中获取样本吗？也许通过生成的数量块？
有两种方法可以做到这一点。
首先，如果程序足够通用，只需传入零大小的数据，这样后验就是先验。例如，您给出的回归示例中的 N = 0 将起作用（以及正确的零大小 x 和 y）。
其次，您可以在生成的数量块中编写一个纯 Monte Carlo 生成器（不使用 MCMC）。类似于：
generated quantities {
real<lower = 0> sigma_sim = cauchy_rng(0, 2); // wide tail warning!
real beta_sim = normal_rng(0, 10);
real alpha_sim = normal_rng(0, 20);
}
第二种方法效率更高，因为它可以方便地抽取独立样本，而不必执行任何 MCMC。
今天早上我在公共汽车上想到了如何做到这一点的答案。当然，当我写完它时，@Bob Carpenter 发布了我正在寻找的解决方案。相比之下，我的方法相当麻烦和笨拙，但它 确实 有效。
我们需要做的就是指定反映实际先验但永远不会向下游传递到似然函数的先验。
因此在上面的示例中，我们需要做的就是在模型字符串中创建这些镜像变量。我们称它们为 p_alpha、p_beta 和 p_sigma。这些将是 alpha、beta 和 sigma 的类似物，但不会出现在任何似然函数中。
请注意，我们必须在 parameters{} 块和 model{} 块中创建这些变量。
ms <- "
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
// priors to sample from
real p_alpha;
real p_beta;
real p_sigma;
// real priors
real alpha;
real beta;
real<lower=0> sigma;
}
model {
vector[N] mu;
// priors to sample from
p_sigma ~ cauchy(0, 2);
p_beta ~ normal(3,1); // for didactic purposes
p_alpha ~ normal(0,100);
// actual priors
sigma ~ cauchy(0, 2);
beta ~ normal(0,10);
alpha ~ normal(0,100);
// likelihood
for ( i in 1:N ) {
mu[i] = alpha + beta * x[i];
}
y ~ normal(mu, sigma);
}
"
请注意，镜像参数的分布规范应与我在 p_alpha/alpha 和 p_sigma/[=20 中所做的实际先验相匹配=].出于教学目的，我故意使 p_beta 的中心和散布与 beta 不同，因为我将在下面的同一张图上绘制它们。
现在运行再次模型
fit1 <- stan(model_code = ms,
data = data_reg,
chains = 1,
warmup = 1e3,
iter = 2e3)
并提取样本
post <- as.data.frame(extract(fit1, pars = c("p_alpha", "p_beta", "p_sigma", "alpha", "beta", "sigma")))
head(post)
# output
p_alpha p_beta p_sigma alpha beta sigma
1 -81.44259 3.275672 -1.1416369 3.121382 2.499459 2.354001
2 161.03740 3.694711 0.2989131 3.648288 2.335520 2.140973
3 126.58106 3.495947 -2.0027929 3.846835 2.266247 3.037055
4 18.55785 3.283425 -0.4045153 2.903958 1.854639 1.807591
5 103.02826 5.213568 -18.3721863 3.980290 1.725396 2.178264
6 49.50477 1.737679 6.5971377 4.209471 2.535044 2.941958
这是作为单独图的先验和后验
所以现在我们有相同数据帧中相同参数的原始先验和后验。
现在如果我们想把先验和后验放在同一个图上怎么办？
首先将两个参数 p_beta 和 beta 放入数据框，使其成为长格式，以便估计值在一列中，分布（先验与后验）在另一列中。
library(dplyr)
betaDF <- post %>% dplyr::select(grep("^.*beta$", names(.))) %>%
gather(key = source, value = estimate) %>%
transform(source = factor(ifelse(source == "p_beta", "prior", "posterior"), levels = c("prior", "posterior")))
现在绘制它
ggplot(betaDF, aes(x = estimate, fill = source)) +
geom_density(alpha = 0.3) +
coord_cartesian(xlim = c(-5,10)) +
labs(x = "beta")
我想根据这些参数的先验值绘制 stan 模型的参数估计直方图。我尝试通过 运行 在 stan 中构建一个模型，用 ggplot2 绘制它，然后使用 R 的随机生成函数（例如 rnorm()、rbinom()）覆盖先验分布的近似值，但是我 运行 遇到许多缩放问题，这些问题使图表无法正确显示。
我在想一个更好的方法是直接从先验分布中抽样，然后根据参数估计值绘制这些样本，但是 运行 建立一个完整的单独模型 just 从先验中抽样似乎非常耗时。我想知道是否有一种方法可以在现有模型中或与之并行地执行此操作。
这是一个示例脚本。
# simulate linear model
a <- 3 # intercept
b <- 2 # slope
# data
x <- rnorm(28, 0, 1)
eps <- rnorm(28, 0, 2)
y <- a + b*x + eps
# put data into list
data_reg <- list(N = 28, x = x, y = y)
# create the model string
ms <- "
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
}
model {
vector[N] mu;
sigma ~ cauchy(0, 2);
beta ~ normal(0,10);
alpha ~ normal(0,100);
for ( i in 1:N ) {
mu[i] = alpha + beta * x[i];
}
y ~ normal(mu, sigma);
}
"
# now fit the model in stan
fit1 <- stan(model_code = ms, # model string
data = data_reg, # named list of data
chains = 1, # number of Markov chains
warmup = 1e3, # number of warmup iterations per chain
iter = 2e3) # show progress every 'refresh' iterations
# extract the sample estimates
post <- extract(fit1, pars = c("alpha", "beta", "sigma"))
# now for the density plots. Write a plotting function
densFunct <- function (parName) {
g <- ggplot(postDF, aes_string(x = parName)) +
geom_histogram(aes(y=..density..), fill = "white", colour = "black", bins = 50) +
geom_density(fill = "skyblue", alpha = 0.3)
return(g)
}
# plot
gridExtra::grid.arrange(grobs = lapply(names(postDF), function (i) densFunct(i)), ncol = 1)
现在我明白我可以通过简单地从模型字符串中省略可能性来从先验中采样，就像这样
ms <- "
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
}
model {
sigma ~ cauchy(0, 2);
beta ~ normal(0,10);
alpha ~ normal(0,100);
}
"
但是有什么方法可以从第一个模型中的先验中获取样本吗？也许通过生成的数量块？
有两种方法可以做到这一点。
首先，如果程序足够通用，只需传入零大小的数据，这样后验就是先验。例如，您给出的回归示例中的 N = 0 将起作用（以及正确的零大小 x 和 y）。
其次，您可以在生成的数量块中编写一个纯 Monte Carlo 生成器（不使用 MCMC）。类似于：
generated quantities {
real<lower = 0> sigma_sim = cauchy_rng(0, 2); // wide tail warning!
real beta_sim = normal_rng(0, 10);
real alpha_sim = normal_rng(0, 20);
}
第二种方法效率更高，因为它可以方便地抽取独立样本，而不必执行任何 MCMC。
今天早上我在公共汽车上想到了如何做到这一点的答案。当然，当我写完它时，@Bob Carpenter 发布了我正在寻找的解决方案。相比之下，我的方法相当麻烦和笨拙，但它 确实 有效。
我们需要做的就是指定反映实际先验但永远不会向下游传递到似然函数的先验。
因此在上面的示例中，我们需要做的就是在模型字符串中创建这些镜像变量。我们称它们为 p_alpha、p_beta 和 p_sigma。这些将是 alpha、beta 和 sigma 的类似物，但不会出现在任何似然函数中。
请注意，我们必须在 parameters{} 块和 model{} 块中创建这些变量。
ms <- "
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
}
parameters {
// priors to sample from
real p_alpha;
real p_beta;
real p_sigma;
// real priors
real alpha;
real beta;
real<lower=0> sigma;
}
model {
vector[N] mu;
// priors to sample from
p_sigma ~ cauchy(0, 2);
p_beta ~ normal(3,1); // for didactic purposes
p_alpha ~ normal(0,100);
// actual priors
sigma ~ cauchy(0, 2);
beta ~ normal(0,10);
alpha ~ normal(0,100);
// likelihood
for ( i in 1:N ) {
mu[i] = alpha + beta * x[i];
}
y ~ normal(mu, sigma);
}
"
请注意，镜像参数的分布规范应与我在 p_alpha/alpha 和 p_sigma/[=20 中所做的实际先验相匹配=].出于教学目的，我故意使 p_beta 的中心和散布与 beta 不同，因为我将在下面的同一张图上绘制它们。
现在运行再次模型
fit1 <- stan(model_code = ms,
data = data_reg,
chains = 1,
warmup = 1e3,
iter = 2e3)
并提取样本
post <- as.data.frame(extract(fit1, pars = c("p_alpha", "p_beta", "p_sigma", "alpha", "beta", "sigma")))
head(post)
# output
p_alpha p_beta p_sigma alpha beta sigma
1 -81.44259 3.275672 -1.1416369 3.121382 2.499459 2.354001
2 161.03740 3.694711 0.2989131 3.648288 2.335520 2.140973
3 126.58106 3.495947 -2.0027929 3.846835 2.266247 3.037055
4 18.55785 3.283425 -0.4045153 2.903958 1.854639 1.807591
5 103.02826 5.213568 -18.3721863 3.980290 1.725396 2.178264
6 49.50477 1.737679 6.5971377 4.209471 2.535044 2.941958
这是作为单独图的先验和后验
所以现在我们有相同数据帧中相同参数的原始先验和后验。
现在如果我们想把先验和后验放在同一个图上怎么办？
首先将两个参数 p_beta 和 beta 放入数据框，使其成为长格式，以便估计值在一列中，分布（先验与后验）在另一列中。
library(dplyr)
betaDF <- post %>% dplyr::select(grep("^.*beta$", names(.))) %>%
gather(key = source, value = estimate) %>%
transform(source = factor(ifelse(source == "p_beta", "prior", "posterior"), levels = c("prior", "posterior")))
现在绘制它
ggplot(betaDF, aes(x = estimate, fill = source)) +
geom_density(alpha = 0.3) +
coord_cartesian(xlim = c(-5,10)) +
labs(x = "beta")