初学WinBUGS，请教一个问题。这是我的模型，check model 时候说“expected camma”，不知道什么原因。能帮忙看看嘛？此外我想通过语句算DIC，能帮忙写具体代码吗？谢谢啦！
model {
# trick using ones for specifying the CTB and estimating parameters
C <- 10000 # this just has to be large enough to ensure all pp's < 1
for (i in 1:N) {
ones <- 1
ones ~ dbern(pp)
pp <- (pi*a/b*pow(y/b,a*p-1)/pow(1+pow(y/b,a),p+q)+(1-pi)*a/k/b*pow(y/k/b,a*p-1)/pow(1+pow(y/b/k,a),p+q))/(exp(loggam(p))*exp(loggam(q))/exp(loggam(p+q)))/ C
#likelihood for observaion y using for goodness of fit measure logl
l<-log((pi*a/b*pow(y/b,a*p-1)/pow(1+pow(y/b,a),p+q)+(1-pi)*a/k/b*pow(y/k/b,a*p-1)/pow(1+pow(y/b/k,a),p+q))/(exp(loggam(p))*exp(loggam(q))/exp(loggam(p+q))))
#frequency for y use for A^2 measure
v<- rank(y,i)/N
}
#prior
a~dgamma(0.01,0.01)
b~dgamma(0.01,0.01)
p~dgamma(0.01,0.01)
q~dgamma(0.01,0.01)
pi~dgamma(0.01,0.01)
k~dgamma(0.01,0.01)
#logl
logl<-sum(l)
#DIC
#A^2
for(i in 1:N){
vv<-(2*i-1)*(log(v)+log(1-v[N-i+1]))
}
A2<- -N-1/N*sum(vv)
}
#data
list(y=c(290.4, 537.19, 756.8, 769.19, 787.69, 796.18, 933.62, 967.97, 1010.56, 1017.4, 1033.49, 1034.33, 1056.93, 1124.09, 1165.73, 1248.49, 1268.24, 1284.56, 1363.85, 1436.2, 1445.96, 1469.48, 1507.47, 1662.36, 1674.58, 1690.91, 1739.96, 1776.56, 1932.09, 1975.89, 2099.79, 1217.64, 2202.96, 2222.8, 2255.72, 2274.61, 2328.64, 2384.37, 2847.83, 2947.04, 2948.35, 3036.51, 3287.68, 3331.62, 3416.67, 3604.66, 3671.16, 3739.3, 3941.3, 4017.01, 4100, 4166.98, 4355.02, 5117.93, 5335.96, 5453.02, 5568.96, 5761.83, 6161.81, 6348.69, 6859.37, 7972.2, 8028.32, 10047.22, 10560.1, 11179.54, 11461.39, 14538.13, 14789.81, 17186.09, 18582.57, 22857.33, 23177.85, 23446.13, 28409.82, 57612.82, 59582.78, 113164.7, 123228.9, 626402.8),N=80)
#initial values
list(a=1,b=2,p=2,q=5,pi=0.5,k=2)
扫码加好友，拉您进群
全部版块
我的主页
收藏
