如何在 rstan 中一次将发行版适配到多个站点？

Question

可能有一个非常简单的解决方案，但我还不够。我在 rstan 中有代码，其中我将帕累托分布拟合到数据，并且我正在寻找分布的 alpha 参数：

data {
    int<lower=0> N;
    real<lower=0> x_min;
    vector<lower=0>[N] x;
}

parameters {
    // Pareto density
    real<lower=0, upper=5> alpha;
}

model {
    // Prior: Pareto density
    alpha ~ lognormal(1, 1) T[0, 5];
    // x_min uniform on its interval.
    
    // Likelihood: Pareto density
    x ~ pareto(x_min, alpha);
}

这对于单个测量站点非常有效，其中

       N<-length(data$x)
       x<-data$x
       x_min<-x_min
       
       stan_dat<-list(N=N,x=x,x_min=x_min)

但是，我有多个站点（因此，data$site_label 存在）。我可以在 R 的 for 循环中将分布拟合到每个站点并单独获取输出，但是有没有一种方法可以在一个模型中完成所有这些工作？我该如何编码？

我最初使用了 for 循环，其中我分别为每个站点分配了一个分布。它有效，但我想尝试用不同的方式来做

Answer 1

您可以向 Stan 数据添加一个数组，用于索引每个观察的站点。那么

alpha

将是一个参数向量（每个站点一个）而不是单个参数。像

pareto

这样的函数可以轻松矢量化，因此您仍然可以使用单个采样语句将每个观察结果与正确的参数相匹配。

这就是 Stan 代码的样子。我假设

x_min

在所有站点上都是相同的；如果不是，也可以更改为向量。

data {
  int<lower=0> N;
  real<lower=0> x_min;
  vector<lower=0>[N] x;
  int<lower=0,upper=N> N_sites;
  int<lower=1,upper=N_sites> site[N];
}

parameters {
  // Pareto density
  vector<lower=0, upper=5>[N_sites] alpha;
}

model {
  // Prior: Pareto density
  alpha ~ lognormal(1, 1) T[0, 5];

  // Likelihood: Pareto density
  x ~ pareto(x_min, alpha[site]);
}

这里有一些 R 代码，演示了模型为每个站点选取了正确的参数值。我模拟了五个站点的数据。

library(tidyverse)
library(rstan)
library(sads)
theme_set(theme_bw())

# Simulate five sites, each with known alpha
set.seed(24601)
true.values = data.frame(site = 1:5, true.alpha = runif(5, 0, 5))
data = pmap_dfr(
  true.values,
  function(site, true.alpha) {
    data.frame(x = rpareto(100, true.alpha, 1))
  },
  .id = "site"
) %>%
  mutate(site = as.numeric(site))

# Fit a Stan model to the simulated data
stan_dat = list(N = nrow(data),
                x = data$x,
                x_min = 1,
                N_sites = max(data$site),
                site = data$site)
stan_fit = stan("model.stan", data = stan_dat, chains = 3)

# Do the fitted parameters match the true values?  Yes!
rstan::extract(stan_fit, "alpha")$alpha %>%
  data.frame() %>%
  rownames_to_column("n") %>%
  pivot_longer(cols = matches("X"), names_to = "site", names_prefix = "X",
               values_to = "alpha") %>%
  ggplot(aes(x = alpha)) +
  geom_histogram() +
  geom_vline(data = true.values %>% mutate(site = as.character(site)),
             aes(xintercept = true.alpha), color = "red") +
  facet_wrap(~ site) +
  labs(x = "Estimated value of alpha", y = "Number of draws")