时间戳

首页 > 解决方案 > 在不运行单独模型的情况下从先前采样

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

问题描述

我想根据这些参数的先验绘制来自 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);
  }
"

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

标签： rstanrstan

解决方案

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

首先，如果程序足够通用，只需传入零大小的数据，这样后验就是先验。例如，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。

推荐阅读