改进用于分层logit / probit模型的Stan程序的速度

问题描述

我使用以下代码来估计分层Logit模型。

data {
  int<lower=1> D;   // number of covariates
  int<lower=0> N;   // number of observations
  int<lower=1> L;   // number of groups (levels)
  int<lower=0,upper=1> y[N];  // binary response
  int<lower=1,upper=L> ll[N];    // group ids
  row_vector[D] x[N];   //covariates
}
parameters {
  real mu[D];
  real<lower=0> sigma[D];
  vector[D] beta[L];
}
model {
  
  mu ~ normal(0,100);  // prior for mu
//  sigma ~ uniform(0,100) ; // prior for sigma
  for (l in 1:L) beta[l] ~ normal(mu,sigma);
  
  {
    vector[N] x_beta_ll;
    for (n in 1:N) x_beta_ll[n] = x[n] * beta[ll[n]];
    y ~ bernoulli_logit(x_beta_ll);
  }
}

对于我的数据，N = 264,713（总观察数），L = 10,000（组数），D = 26（协变量数）。然后，我收到了此消息，并且花了很长时间才能采样：

Chain 1: Gradient evaluation took 0.131 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1310 seconds.
Chain 1: Adjust your expectations accordingly!

我想知道是否有提高速度的好方法。顺便说一句，这是我使用的R代码：

fit <- stan(
  file = "hierarchical_logit.stan",# Stan program
  data = data,# named list of data
  chains = 1,# number of Markov chains
  warmup = 1000,# number of warmup iterations per chain
  iter = 2000,# total number of iterations per chain
  cores = 5,# number of cores (could use one per chain)
  verbose = TRUE
)

我希望得到专家的提示。

谢谢。

解决方法

第一件事是，在进行logit转换时，您的先验情况看起来很糟糕。这取决于您要提供的变量，但是我猜想它将把这些变量的几乎所有密度都置于0和1。对于连续变量，如果您一开始就提供标准变量，则将暗示这些概率在完全没有空间的情况下从0步进到1，并且类似地，对于拦截项，大多数密度也将为0和1。值得考虑这些-并尝试模拟它，例如在R中，任何拦截项的先验是hist(boot::inv.logit(rnorm(1000,100)),breaks = 100))。例如，在拦截项上更健康，更没有先验的先验将是normal(0,1.5)之类的东西，在概率标度的[0,1]之间大致相同（如McElreath所建议）。类似的想法和模拟将帮助您确定其他变量的先验条件（我不知道它们是什么）。根据变量/数据的外观，提供一些常识性的先验可能会很快，并且无论如何都很重要。

除非有某些特殊原因，否则我要做的第一件事就是使用brms或rstanarm运行模型。这是一个标准模型，它们将自动运行结构良好的模型，您可以只使用lmer型语法。如果那仍然太慢了，那么您将需要做一些额外的工作。

如果您想自己做（/ brms太慢了），请确保传递标准化变量，或者在模型中标准化它们（https://mc-stan.org/docs/2_25/stan-users-guide/standardizing-predictors-and-outputs.html）。

目前，您有一个中心化的参数化-非中心化的参数化可能会快得多，因此，如果速度较慢，则值得尝试一下（Stan手册中有一节）。看起来像这样：

parameters {
  real mu[D];
  real<lower=0> sigma[D];
  vector[D] beta_z[L];
}
transformed parameters {
  vector[D] beta[L];
  for(l in 1:L){
    beta[l] = mu + sigma * beta_z[l];
  }
}
model {
...as before,but replace beta[l] ~ normal(mu,sigma) with...
for (l in 1:L) beta_z[l] ~ normal(0,1);
}

最后，您可以编写一个使用链内并行化来提高速度的模型。但是，这需要一些努力，在您发现更简单的步骤仍然太慢之前，我不会这样做。有关工作示例，请参见https://mc-stan.org/users/documentation/case-studies/reduce_sum_tutorial.html。

此外，对于有关Stan的问题，Stan discourse非常好-我可能还没有做过其他事情。