3

我正在尝试使用 RSTAN 拟合随机效应模型。我的设计矩阵有 198 列。它是如此之宽，因为我的原始数据框是一堆因子变量，我将其转换为二进制指标以尝试在 STAN 中拟合模型。我可以使用从一个或两个预测变量转换而来的几列来拟合模型，但完成 1/2 的采样需要 10 个小时。

这是我用来尝试拟合模型的 STAN 代码（基本线性模型）。我试过矢量化，但也许有办法进一步优化？另外，为什么要花这么长时间的直觉是什么？

     data {
      int<lower=0> N;
      int<lower=0> J; 
      int<lower=0> K;
      int<lower=1,upper=J> geo[N];
      matrix[N,K] X;
      vector[N] y;
    }
    parameters {
      vector[J] a;
      vector[K] B;
      real mu_a;
      real<lower=0,upper=100> sigma_a;
      real<lower=0,upper=100> sigma_y;
    } 
    model {
      vector[N] y_hat;
      for (i in 1:N)
        y_hat[i] <- a[geo[i]];
      mu_a ~ normal(0, 1);
      a ~ normal(0, 1);
      y ~ normal(mu_a + sigma_a * y_hat + X * B, sigma_y);
    }
4

1 回答 1

8

您可以做些什么来加速这个模型的问题与如何更有效地对其进行采样的问题交织在一起。关于为什么要花这么长时间的直觉可能与和之间的先前依赖有关a（sigma_a并且在较小程度上mu_a）。

  1. 当sigma_a对于一些迭代来说很小时，元素a必须接近mu_a
  2. 当sigma_a对于某些迭代来说很大时， 的元素a可以远离mu_a。

由于 Stan 只有一个步长参数可供操作，因此很难找到适合这两种情况的步长。充其量，你得到的步长足够小以在前一种情况下保持准确性，但墙时间会受到不利影响，因为在后一种情况下它必须采取许多步长小的越级步骤。

对于这样的模型，我们通常建议重新参数化

data {
  int<lower=0> N;
  int<lower=0> J; // 47 apparently but don't hardcode it
  int<lower=1,upper=J> geo[N];
  vector[N] y;
}  
parameters {
  vector[J] a;
  real mu_a;
  real<lower=0,upper=100> sigma_a;
  real<lower=0,upper=100> sigma_y;
} 
model {
  vector[N] y_hat;
  for (i in 1:N)
    y_hat[i] <- a[geo[i]];
  mu_a ~ normal(0, 1);
  a ~ normal(0, 1);
  y ~ normal(mu_a + sigma_a * y_hat, sigma_y);
}
于 2015-03-23T02:54:28.230 回答