全部版块 我的主页
论坛 计量经济学与统计论坛 五区 计量经济学与统计软件
1783 2
2013-12-14
Could you please take a look at the following WinBUGs code from WinBUGS Manual Volume 1, and explain what  Y, Xa, Xs mean?, Thanks a lot!
=======================================

Dogs: loglinear model for binary data


Lindley (19??) analyses data from Kalbfleisch (1985) on the Solomon-Wynne experiment on
dogs, whereby they learn to avoid an electric shock. A dog is put in a compartment, the lights are
turned out and a barrier is raised, and 10 seconds later an electric shock is applied. The results
are recorded as success (Y = 1 ) if the dog jumps the barrier before the shock occurs, or failure
(Y = 0) otherwise.
Thirty dogs were each subjected to 25 such trials. A plausible model is to suppose that a dog
learns from previous trials, with the probability of success depending on the number of previous
shocks and the number of previous avoidances. Lindley thus uses the following model
pj = Axj Bj-xj
for the probability of a shock (failure) at trial j, where xj = number of success (avoidances) before
trial j and j - xj = number of previous failures (shocks). This is equivalent to the following log linear
model
log pj = axj + b ( j-xj )
Hence we have a generalised linear model for binary data, but with a log-link function rather than
the canonical logit link. This is trivial to implement in BUGS:
model
{
for (i in 1 : Dogs) {
xa[i, 1] <- 0; xs[i, 1] <- 0 p[i, 1] <- 0
for (j in 2 : Trials) {
xa[i, j] <- sum(Y[i, 1 : j - 1])
xs[i, j] <- j - 1 - xa[i, j]
log(p[i, j]) <- alpha * xa[i, j] + beta * xs[i, j]
y[i, j] <- 1 - Y[i, j]
y[i, j] ~ dbern(p[i, j])
}
}
alpha ~ dnorm(0, 0.00001)I(, -0.00001)
beta ~ dnorm(0, 0.00001)I(, -0.00001)
A <- exp(alpha)
B <- exp(beta)
}

二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

全部回复
2013-12-14 01:10:06
Attached here is the stan code by Andrew Gelman for Dogs: loglinear model for binary data. I am  just wondering why he use bernoulli function for loglinear model for binary data? Any suggestion and instruction would be appreciated!

Thanks
==============================================================
# compared with JAGS version in
# the R package BUGSExamples (https://r-forge.r-project.org/R/?group_id=882)
data {
  int<lower=0> Ndogs;
  int<lower=0> Ntrials;
  int Y[Ndogs, Ntrials];
}
transformed data {
  int y[Ndogs, Ntrials];
  int xa[Ndogs, Ntrials];
  int xs[Ndogs, Ntrials];
  for (dog in 1:Ndogs) {
    xa[dog, 1] <- 0;
    xs[dog, 1] <- 0;
    for (trial in 2:Ntrials) {
      for (k in 1:(trial - 1))
        xa[dog, trial] <- xa[dog, trial] + Y[dog, k];
      xs[dog, trial] <- trial - 1 - xa[dog, trial];
    }
  }
  for (dog in 1:Ndogs) {
    for (trial in 1:Ntrials) {
      y[dog, trial] <- 1 - Y[dog, trial];
    }
  }
}
parameters {
  real<upper= -0.00001> alpha;
  real<upper= -0.00001> beta;
}
model {
  alpha ~ normal(0.0, 316.2);
  beta  ~ normal(0.0, 316.2);
  for(dog in 1:Ndogs)  
    for (trial in 2:Ntrials)  
      y[dog, trial] ~ bernoulli(exp(alpha * xa[dog, trial] + beta * xs[dog, trial]));
}
generated quantities {
  real A;
  real B;
  A <- exp(alpha);
  B <- exp(beta);
}
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2018-2-7 12:19:57
楼主你好,我想问问代码里的I(, -0.00001)是什么意思,请问能否指教一下?
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

相关推荐
栏目导航
热门文章
推荐文章

说点什么

分享

扫码加好友,拉您进群
各岗位、行业、专业交流群