如何在没有 for 循环的情况下使用 RStan？

Question

有没有一种方法可以更有效地在 RStan 中执行以下计算？

我只提供了所需的最少代码量：

parameters {
  real beta_0;
  real beta_1;
}     
model {
  vector [n] p_i = exp(beta_0 + beta_1*x)/[1 + exp(beta_0 + beta_1*x)];
  y ~ bernoulli(p_i);
  /* Likelihood:
  for(i in 1:n){
    p_i[i] = exp(beta_0 + beta_1*x[i])/(1 + exp(beta_0 + beta_1*x[i]));
    y[i] ~ bernoulli(p_i[i]);
  */}
// Prior:
  beta_0 ~ normal(m_beta_0, s_beta_0);
  beta_1 ~ normal(m_beta_1, s_beta_1);
}

我收到以下错误消息："Matrix expression elements must be type row_vector and row vector expression elements must be int or real, but found element of type vector"。如果我使用 for 循环（已注释掉），代码工作正常，但我想在我的代码中限制 for 循环的使用。在上面的代码中，x 是一个长度为 n 的向量。

另一个例子：

parameters {
  real gamma1;
  real gamma2;
  real gamma3;
  real gamma4;
}
model {
// Likelihood:
  real lambda;
  real beta;
  real phi;
  for(i in 1:n){
    lambda = exp(gamma1)*x[n_length[i]]^gamma2;
    beta = exp(gamma3)*x[n_length[i]]^gamma4;
    phi = lambda^(-1/beta);
    y[i] ~ weibull(beta, phi);
  }
  //y ~ weibull(exp(gamma1)*x^gamma2, exp(gamma3)*x^gamma4); //cannot raise a vector to a power
// Prior:
  gamma1 ~ normal(m_gamma1, s_gamma1);
  gamma2 ~ normal(m_gamma2, s_gamma2);
  gamma3 ~ normal(m_gamma3, s_gamma3);
  gamma4 ~ normal(m_gamma4, s_gamma4);
}

上面的代码有效，但注释掉的似然计算不起作用，因为我 "cannot raise a vector to a power"（但你可以在 R 中）。再一次，我不想被迫使用 for 循环。上面代码中，n_length，是一个长度为n的向量。

最后一个例子。如果我想从 R 中的正态分布中抽取 10000 个样本，我可以简单地指定

rnorm(10000, mu, sigma)

但是在 RStan 中，我必须使用 for 循环，例如

parameters {
      real mu;
      real sigma;
    }
generated quantities {
  vector[n] x;
  for(i in 1:n) {
    x[i] = normal_rng(mu, sigma);
  }
}

我可以做些什么来加速我的 RStan 示例？

Answer 1

这行代码：

vector [n] p_i = exp(beta_0 + beta_1*x)/[1 + exp(beta_0 + beta_1*x)];

在 Stan 语言中不是有效语法，因为方括号仅用于索引。它可以改为

vector [n] p_i = exp(beta_0 + beta_1*x) ./ (1 + exp(beta_0 + beta_1*x));

它利用元素除法运算符，或者更好

vector [n] p_i = inv_logit(beta_0 + beta_1*x);

在这种情况下，y ~ bernoulli(p_i); 很可能会起作用。更好的是，就做

y ~ bernoulli_logit(beta_0 + beta_1 * x);

它会以数字稳定的方式为您进行转换。您也可以使用 bernoulli_logit_glm，它稍微快一些，尤其是对于大型数据集。

在 Stan 2.19.x 中，我认为您可以从生成数量块中的概率分布中提取 N 个值。但是你太担心 for 循环了。 Stan 程序被转译为 C++，其中循环速度很快，Stan 语言中几乎所有接受向量输入并产生向量输出的函数实际上都涉及 C++ 中的相同循环，就好像您自己完成了循环一样。

如何在没有 for 循环的情况下使用 RStan？

How can I use RStan without for loops?

r

stan

rstan