在 rstan 中变换变量(贝叶斯分析)

Transforming variables in rstan (bayesian analysis)

提问人:gowerc 提问时间:1/25/2017 最后编辑:gowerc 更新时间:1/25/2017 访问量:1048

问:

我是贝叶斯分析的新手,正在尝试使用 rstan 来估计后验密度分布。该练习试图重新创建我的大学使用 stan 提供给我们的示例,但我对如何正确转换变量有点困惑。我当前的代码运行没有错误,但结果与大学给我们的结果不匹配(尽管很接近),下图清晰可见,黑色的斯坦估计值。我通过查阅手册并将随机位拼凑在一起来使代码工作,但特别是我不太确定为什么需要以及伽玛的转换是否确实正确。任何指导将不胜感激!target

enter image description here

斯坦码

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
)

达特

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)所指出的,问题是我从θ中推导出伽玛,而伽玛应该是相反的,因为伽玛是我的随机变量。正确的 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 );
}
R 贝叶斯 斯坦

评论

0赞 Ben Goodrich 1/25/2017
该区块仅执行一次;该块在马尔可夫链的每次迭代中都会被执行多次。transformed datatransformed parameters

答:

3赞 Ben Goodrich 1/25/2017 #1

作为提示,尝试将 (马) 放入块中,然后根据您一开始给出的分布在块中声明和定义。gamparametersthetatransformed parameters

Stan 的初学者通常认为 Stan 具有逻辑地计算 Stan 程序含义的能力,而实际上它被相当字面地转译为 C++,并且代码行和块被一遍又一遍地执行。transformed parametersmodel

原始参数是(马)还是原始参数有区别的原因与变量变化原理有关。如果您真的愿意,如果将雅可比行列式项(以对数单位表示)添加到 ,则可以获得与原始参数化相同的结果,但是通过将 (马) 移动到块和块来更容易避免这种情况。有关变量变化原理的详细信息,请参阅此案例研究gamthetatargetgamparametersthetatransformed parameters

评论

0赞 gowerc 1/25/2017
这修好了!在我标记为答案之前,您能否帮助解释为什么这会修复它?我本来以为将 θ 作为 γ 的函数或将 gamma 作为 θ 的函数应该是可以互换的?
0赞 gowerc 1/25/2017
话虽如此,我现在在想这是因为 gamma 是我们从正态分布中采样的 rv,而 theta 是该 rv 的函数。