r - 在不运行单独模型的情况下从先前采样

Question

我想根据这些参数的先验来绘制来自 stan 模型的参数估计的直方图。我尝试通过在 stan 中运行模型，用 ggplot2 绘制它，然后使用 R 的随机生成器函数(例如 rnorm() ， rbinom() )覆盖先验分布的近似值来做到这一点，但我遇到了许多缩放问题，这些问题使图表不可能看起来正确。

我在想一个更好的方法是直接从先验分布中采样，然后根据参数估计绘制这些样本，但是运行一个完整的单独模型来从先验中采样似乎非常耗时。我想知道是否有一种方法可以在现有模型内或与之并行执行此操作。

这是一个示例脚本。

# 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);
  }
"

但是有没有办法从第一个模型中的先验中获取样本？也许通过生成的数量块？

Answer 1

有两种方法可以做到这一点。

首先，如果程序足够通用，只需传入零大小的数据，以便后验是先验。例如，N = 0在您给出的回归示例中(以及正确的零大小的 x 和 y)。

其次，您可以在生成量块中编写纯蒙特卡罗生成器(不使用 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。