在大型稀疏数据帧上拟合随机效应模型时如何加速STAN？

Question

我正在尝试使用RSTAN拟合随机效果模型.我的设计矩阵有198列.它是如此宽,因为我的原始数据帧是一堆因子变量,我转换为二进制指标,以尝试在STAN中拟合模型.我能够使用从一个或两个预测变换器转换的几列来拟合模型,但是花费了10个小时来完成1/2的样本.

这是我用来尝试拟合模型的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);
    }

Answer 1

你能做些什么来加速这个模型的问题与如何更有效地采样的问题交织在一起.关于为什么花费这么长时间的直觉可能与先前a和之间sigma_a(和在较小程度上mu_a)的依赖性有关.

1. 当sigma_a某些迭代很小时,元素a必须接近mu_a
2. 当sigma_a某些迭代很大时,元素a可能很远mu_a.

由于Stan只有一个stepize参数可供操作,因此很难找到适合两种情况的步长.在最好的情况下,你得到一个小到足以保持前一种情况下的准确性的步长,但是壁挂时间会受到不利影响,因为在后一种情况下它必须采用许多跳跃步骤和一小步.

对于像这样的模型,我们通常建议重新参数化 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); }