时间戳

首页 > 解决方案 > 如何在没有 for 循环的情况下使用 RStan？

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

问题描述

有没有办法在 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);
}

我收到以下错误消息：“矩阵表达式元素必须是 row_vector 类型，并且行向量表达式元素必须是 int 或 real，但找到了 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);
}

上面的代码有效，但注释掉的可能性计算不起作用，因为我“不能将向量提升到幂”（但你可以在 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 示例的速度吗？

标签： rstanrstan

解决方案

这行代码：

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++ 中的相同循环，就好像您自己完成了循环一样。

推荐阅读