首页 文章

当从二项式随机变量建模值时，Stan / RStan中if（）条件的问题

提问于
浏览 次
2

我正在尝试使用Stan和R来拟合一个模型，呃，模拟观察到的实现y_i = 16,9,10,13,19,20,18,17,35,55，它们来自二项式分布随机变量，比方说，Y_i，参数m_i（试验次数）和p_i（每次试验中的成功概率） .

yi = c(16, 9, 10, 13, 19, 20, 18, 17, 35, 55)

出于本实验的目的，我将假设所有m_i都是固定的并由m_i = 74,99,58,70,122,77,104,129,308,119给出 .

mi = c(74, 99, 58, 70, 122, 77, 104, 129, 308, 119)

我将使用Jeffrey的先前：\ alpha = 0.5和\ beta = 0.5 .

alpha = 0.5, beta = 0.5

我试着

  • 查找p_i的贝叶斯估计 .

  • 找到p_i的范围（即参数k如下：
    enter image description here

我在2.的尝试是这部分代码：

real k;
        real mx = 0;
        real mn = 0;
        if (p > mx) 
          mx = p;
        if (mn > p) {
          mn = p;
        }
        k = mx - mn;

我的Stan代码如下：

```{stan output.var="BinModBeta"}
  data {
    int <lower = 1> mi[10];
    int <lower = 0> yi[10];
    real <lower = 0> alpha;
    real <lower = 0> beta;
  }

  parameters {
    real <lower = 0, upper = 1> p[10];
  }

  transformed parameters {
    real k;
    real mx = 0;
    real mn = 0;
    if (p > mx) 
      mx = p;
    if (mn > p) {
      mn = p;
    }
    k = mx - mn;
  }

  model {
    yi ~ binomial(mi, p);
    p ~ beta(alpha, beta);
  }


我的R代码如下：

```java
```{r}
library(rstan)

data.in <- list(mi = c(74, 99, 58, 70, 122, 77, 104, 129, 308, 119), yi = c(16, 9, 10, 13, 19, 20, 18, 17, 35, 55), alpha = 0.5, beta = 0.5)
model.fit1 <- sampling(BinModBeta, data=data.in)

print(model.fit1, pars = c("p"), probs=c(0.1,0.5,0.9), digits = 5)


现在，我刚刚开始学习斯坦，所以我真的不确定这是否正确 . 但是，似乎这段代码适用于我的第一个目标（至少，无论我编写的代码似乎都有用......） .   **But my trouble starts when attempting to code my second aim.** 

当我尝试编译上面的Stan代码时，我收到以下错误：


![enter image description here](https://i.stack.imgur.com/15b3X.png)


现在，基于这个错误信息，似乎我的问题是由于p是10个实数值的向量而不是单个实数的事实 . 但是，由于我对斯坦缺乏经验，我完全不确定如何解决这个问题 . 

如果有人愿意花时间帮我解决这个问题，我将不胜感激 .
r statistics stan rstan

1 回答

  • 2

    这是我要做的：

    model <- "
data {
  int <lower = 1> mi[10];
  int <lower = 0> yi[10];
  real <lower = 0> alpha;
  real <lower = 0> beta;
}

parameters {
  real <lower = 0, upper = 1> p[10];
}

model {
    p ~ beta(alpha, beta);        // Prior
    yi ~ binomial(mi, p);         // Likelihood
}

generated quantities {
    real k;
    k = max(p) - min(p);
}
"

library(rstan);
yi = c(16, 9, 10, 13, 19, 20, 18, 17, 35, 55);
mi = c(74, 99, 58, 70, 122, 77, 104, 129, 308, 119);
fit <- stan(
    model_code = model,
    data = list(mi = mi, yi = yi, alpha = 0.5, beta = 0.5))
fit;
#Inference for Stan model: 6a01a3b25656e1b18183baf19183abf7.
#4 chains, each with iter=2000; warmup=1000; thin=1;
#post-warmup draws per chain=1000, total post-warmup draws=4000.
#
#         mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
#p[1]     0.22    0.00 0.05    0.13    0.19    0.22    0.25    0.32  4000    1
#p[2]     0.10    0.00 0.03    0.05    0.07    0.09    0.11    0.16  4000    1
#p[3]     0.18    0.00 0.05    0.09    0.14    0.17    0.21    0.28  4000    1
#p[4]     0.19    0.00 0.05    0.11    0.16    0.19    0.22    0.29  4000    1
#p[5]     0.16    0.00 0.03    0.10    0.14    0.16    0.18    0.22  4000    1
#p[6]     0.26    0.00 0.05    0.17    0.23    0.26    0.30    0.37  4000    1
#p[7]     0.18    0.00 0.04    0.11    0.15    0.17    0.20    0.25  4000    1
#p[8]     0.13    0.00 0.03    0.08    0.11    0.13    0.15    0.20  4000    1
#p[9]     0.11    0.00 0.02    0.08    0.10    0.11    0.13    0.15  4000    1
#p[10]    0.46    0.00 0.04    0.38    0.43    0.46    0.49    0.55  4000    1
#k        0.38    0.00 0.05    0.28    0.34    0.38    0.41    0.47  4000    1
#lp__  -530.01    0.05 2.26 -535.38 -531.33 -529.65 -528.37 -526.69  1782    1
#
#Samples were drawn using NUTS(diag_e) at Tue Apr 24 22:02:07 2018.
#For each parameter, n_eff is a crude measure of effective sample size,
#and Rhat is the potential scale reduction factor on split chains (at
#convergence, Rhat=1).

    评论：

    • 我会将涉及计算 k 的部分移到 generated quantities 块中;这与在不同时间执行的不同程序块有关 . 虽然 transformed parameters 块在每个跳跃步骤执行一次，但 generated quantities 块仅在每个样本绘制时执行一次 . 因此重新计算 k 的开销会减少 . 参见例如here了解详情 . 请注意， pi 后部密度的不确定性会正确地传播到 k .

    • 在计算 k 时，可以使用 Stan s内部 max ， min 函数 . 这比使用 if 条件确定 pi 的最小值/最大值更快，并且还无需定义 mn 和 mx .

    回复于

相关问题