我是贝叶斯分析的新手,正在尝试使用rstan估算后验密度分布。练习试图重新创建我的大学使用stan给我们的示例,但是我对如何正确地转换变量感到有些困惑。我当前的代码运行没有错误,但结果与大学给我们的结果不符(尽管很接近),为清楚起见,下面的图表以黑色标明了稳定的估算值。通过参考手册并将随机位拼凑在一起,可以使代码正常工作,但是特别是我不太确定为什么需要target以及伽马的转换是否确实正确。任何指导将不胜感激!

模型

r - 转换rstan中的变量(贝叶斯分析)-LMLPHP

斯坦码

data {
  int<lower=0> I;
  int<lower=0> n[I];
  int<lower=0> x[I];
  real<lower=0> a;
  real<lower=0> b;
  real m;
  real<lower=0> p;
}

parameters {
  real<lower=0> lambda;
  real mu;
  real<lower=0, upper=1> theta[I];
}

transformed parameters {
  real gam[I];
  for( j in 1:I)
   gam[j] = log(theta[j] / (1-theta[j])) ;
}


model {
  target +=  gamma_lpdf( lambda |  a, b);
  target +=  normal_lpdf( mu | m , 1/sqrt(p));
  target +=  normal_lpdf( gam | mu, 1/sqrt(lambda));
  target +=  binomial_lpmf( x | n , theta);
}

R代码
library(rstan)
fit <- stan(
  file = "hospital.stan" ,
  data = dat ,
  iter = 20000,
  warmup = 2000,
  chains = 1
)

dat
structure(
  list(
      I = 12L,
      n = c(47, 211, 810, 148, 196, 360, 119,  207, 97, 256, 148, 215),
      x = c(0, 8, 46, 9, 13, 24, 8, 14, 8, 29, 18, 31),
      a = 2,
      b = 2,
      m = 0,
      p = 0.01),
  .Names = c("I", "n", "x", "a", "b", "m", "p")
)

---用解决方案更新---

本·古德里奇(Ben Goodrich)指出的问题是,我是从theta导出伽马,因为伽马是我的随机变量,所以它应该是相反的。正确的stan代码如下。
data {
    int<lower=0> I;
    int<lower=0> n[I];
    int<lower=0> x[I];
    real<lower=0> a;
    real<lower=0> b;
    real m;
    real<lower=0> p;
}

parameters {
    real<lower=0> lambda;
    real mu;
    real gam[I];
}

transformed parameters {
    real<lower=0 , upper=1> theta[I];
    // theta = inv_logit(gam);  // Alternatively can use the in-built inv_logit function which is vectorised
    for(j in 1:I){
        theta[j] = 1 / ( 1 + exp(-gam[j]));
    }
}

model {
    target +=  gamma_lpdf( lambda |  a, b);
    target +=  normal_lpdf( mu | m , 1/sqrt(p));
    target +=  normal_lpdf( gam | mu, 1/sqrt(lambda));
    target +=  binomial_lpmf( x | n ,  theta );
}

最佳答案

提示一下,尝试将gam(ma)放入parameters块中,然后根据您一开始就给出的分布在theta块中声明并定义transformed parameters

Stan的初学者通常认为Stan具有逻辑上解决Stan程序的含义的能力,当确实将Stan程序真正按字面意义转换为C++时,就会反复执行transformed parametersmodel块中的代码行。

无论是gam(ma)还是theta是原始参数,它之所以与众不同,原因与变量更改原理有关。如果确实需要,如果在target中添加了Jacobian行列式项(以对数为单位),则可以通过原始参数化获得相同的结果,但是通过将gam(ma)移至parameters块和theta,可以更轻松地避免这种情况到transformed parameters块。有关变量更改原理的详细信息,请参见case study

关于r - 转换rstan中的变量(贝叶斯分析),我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/41839180/

10-09 13:53