R2jags 中的先验权

问题描述 投票:0回答:0

在固定功率先验模型中,模型设置为:

$$ \pi(p_i \mid lpha,\mathcal{D}_0) \propto L(p_i\mid \mathcal{D}_0)^{w} \pi(p_i) $$

假设事件服从概率为$p_i$的伯努利分布,我想要的是将伯努利似然提高到jags的$w$次方:

mod1<-function(){
  for(i in 1 : N) {
    y[i] ~ dbern(p[i])^w #<-------- 
    logit(p[i]) <- a[group[i]] + t[group[i]]*d[i]
  }
  for (j in 1:3) {
    a[j] ~ dnorm(mu.a, pow(sd.a,-2))
    t[j] ~ dnorm(mu.t, pow(sd.t,-2))
  }
  mu.a~dnorm(0,0.01)
  mu.t~dnorm(0,0.01)
  sd.a ~ dunif(0,10)
  sd.t ~ dunif(0,10)
}

但是 JAGS 不识别这样的代码。

作为替代方案,我编写了以下 STAN 代码:

mod1 <- '
data {
  int<lower = 1> N;
  real w;
  int y[N];
  vector[N] d;
  int group[N];
  }

parameters {
  vector[3] a;
  vector<lower=0>[3] t;
  real mua;
  real mut;
  real sigmaa;
  real sigmat;
}

model {
  real p[N];
  //prior
  mua~normal(0,10);
  mut~normal(0,10);
  sigmaa ~ uniform(0,10);
  sigmat ~ uniform(0,10);
  for (j in 1:3) {
  a[j] ~ normal(mua, sigmaa); 
  t[j] ~ normal(mut, sigmat);
  }
  
  //likelihood

  for (i in 1:N) {
    p[i]=inv_logit(a[group[i]] + t[group[i]]*d[i]);
    target += bernoulli_lpmf(y[i]|(p[i]))*w;
  }
}
'

这有效,但我想在这个项目中使用 JAGS (R2jags),并且想知道是否有与上述 STAN 代码等效的 JAGS。谢谢!

如果你想测试代码,这里是示例数据:

y <- rbinom(160, 1, 0.3)
d <- runif(160, 0, 41)
group <- sample(1:3, 160, replace = TRUE)
w <- rep(0.5, 160)
N=160
mydata <- data.frame(y, d, group, x4,N)
jags stan rjags r2jags
© www.soinside.com 2019 - 2024. All rights reserved.