在固定功率先验模型中,模型设置为:
$$ \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)