2

我正在尝试使用 RStan 代码模拟指数随机变量。

模拟指数随机变量的 R 代码如下所示:

A <- rexp(1000, 4)
B <- rexp(1000, lambda2)
C <- rexp(1000, lambda3)
D <- rexp(1000, lambda4)
E <- rexp(1000, lambda5)
F <- rexp(1000, lambda6)

A + B = AB
C + D = CD

mean(AB)
mean(CD)

如您所见,在 R 中模拟指数随机变量对我来说非常简单。A <- rexp(1000, 4)举个例子:它生成 1000 个结果的随机样本,其中 lambda = 4。然后我可以对这些模拟值进行统计分析(发现手段等)。

现在我想在使用 RStan 时做同样的事情。使用 RStudio 的 R 笔记本功能,可以“插入”不同类型的代码,其中之一是 Stan:

在此处输入图像描述

我有以下 Stan 和 R 代码:

{stan output.var="exponential"}
generated quantities{
  real total;
  real number;

  number = 0.0;
  total = 0.0;
  while(total < 1000) {
    total += exponential_rng(4);
    number += 1.0;
  }
  number -= 1.0;
}

{r}
simul2 <- sampling(exponential, algorithm="Fixed_param")

{r}
print(simul2, pars=c("number"), digits = 5)

但是,当我执行此代码(特别是print命令)时,我得到以下输出:

在此处输入图像描述

我不明白为什么我会得到如此荒谬的统计数据(4000 的意思!?)。如上所示,我的目标是能够获得与在 R 中进行分析时相同的统计数据;换句话说,我的目标是获得与 If I had done 时相同的值A <- rexp(1000, 4),在这种特定情况下,以及X <- rexp(1000, lambda)更普遍的情况下。

显然,我的 Stan 代码是不正确的,所以我非常感谢人们可以花时间解释正确的使用方法。

4

1 回答 1

2

这告诉你 的后验分布number,这只是你的计数器。您的 Stan 代码执行的操作与您的 R 代码不同,但如果您将Nlambdas 作为数据传递,则执行 R 代码的模拟并不难:

data {
  int<lower=1> N;
  real<lower=0> lambda2;
  real<lower=0> lambda3;
  real<lower=0> lambda4;
}
generated quantities {
  real mean_AB;
  real mean_CD;
  {
    vector[N] A;
    vector[N] B;
    vector[N] C;
    vector[N] D;
    for (n in 1:N) {
      A[n] = exponential_rng(4);
      B[n] = exponential_rng(lambda2);
      C[n] = exponential_rng(lambda3);
      D[n] = exponential_rng(lambda4);
    }
    mean_AB = mean(A + B);
    mean_CD = mean(C + D);
  }
}
于 2018-03-20T17:19:57.007 回答