如何在没有for循环的情况下使用RStan？
提问于 2019-06-24 17:16:25
有没有一种方法可以在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，但找到了向量类型的元素”。如果我使用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示例吗？
回答 1
Stack Overflow用户
回答已采纳
发布于 2019-06-25 05:41:55
下面这行代码：
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中，我认为您可以从generated quantities块中的概率分布中提取N个值。但是您太担心for循环了。Stan程序被转换到C++，其中循环速度很快，并且Stan语言中几乎所有接受向量输入和生成向量输出的函数实际上都涉及C++中的相同循环，就好像您自己完成了循环一样。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/56741210
复制相关文章
点击加载更多
腾讯云开发者
