我的任务是使用更新的数据重新运行其他人的模型(我进行了数据更新和格式化)。我了解模型的工作原理但没有编写它,而且它很长而且很详细。最重要的是,我专门使用 JAGS,这是我第一次涉足 WinBUGS,所以界面和错误消息对我来说是新的。该模型需要两天时间才能编译(它确实编译),但随后我收到一个陷阱错误,显示索引超出范围,然后是几十行我无法理解的行。我已经检查了整个模型、for 循环和数据结构,但我看不出哪里有任何索引问题。而且,令人沮丧的是,与 JAGS 不同,WinBUGS 似乎不会告诉您哪个索引超出范围。
我非常感谢解码 Trap 错误消息的任何帮助,以防其中有比我意识到的更多信息来指定索引不匹配的位置。我将在下面发布模型以及陷阱消息,但我看不出发布实际数据的意义,因为没有人会花 2 天时间等待模型编译以重现错误. 相反,我将只发布数据结构,以便每个数据对象的维度都很明显。
一些注意事项:我使用的是以前使用的确切模型(字面意思是复制和粘贴),所以我知道它应该可以工作。之前的数据文件是矩形格式,而我创建和导出的数据文件是 S 格式。我的理解是 WinBUGS 应该能够同时处理这两者,但我想知道这是否可能是问题的一部分。我直接从 WinBUGS 运行模型,手动加载模型和数据文件。我之前尝试使用 r2winbugs 从 R 中执行此操作,但遇到了内存错误和崩溃。
模型:
model{
for(i in 1:n+nzeros){ # loop through the observed and unobserved areas
for(j in 1:J){ #loop through the J fisheries areas
for(t in 1:T){ # loop through the years
y[i,j,t] ~ dpois(mu.y[i,j,t])
mu.y[i,j,t] <- lambda[i,j,t]*x[i,t]
log(lambda[i,j,t]) <- mu.lambda[j]+ theta[z[i],j] + e.lambda[j,t]
}}}
for(i in 1:n+nzeros){
x[i,1] ~ dbern(gamma[i,1])
recruitable[i,1]<-1
for(t in 2:T){
x[i,t] ~ dbern(mu.x[i,t])
mu.x[i,t] <- survived[i,t] + gamma[i,t]*recruitable[i,t]
recruitable[i,t] <- recruitable[i,t-1]*(1-x[i,t-1])
survived[i,t] <- x[i,t-1]*phi[i,t-1]
}}
#priors
for(c in 1:C){
theta[c,1:J] ~ dmnorm(mu.th[], tau.th[,])
for(j in 1:J){
lcapture.clus[c,j]<-mu.lambda[j]+theta[c,j]
pcapture.clus[c,j]<-1-exp(-exp(lcapture.clus[c,j]))
}
}
for(j in 1:J){
mu.th[j] <-0
mu.lambda[j]~dnorm(0, 0.001)
pcapture.overall[j]<-1-exp(-exp(mu.lambda[j]))
for (t in 1:T){
e.lambda[j,t] ~ dnorm(0, tau.e)
lambda.area[j,t]<-mu.lambda[j]+e.lambda[j,t]
p.area[j,t]<-1-exp(-exp(lambda.area[j,t]))
#n.area[j,t]~dbin(p.area[j,t],N.area[j,t])
N.area[j,t]<-n.area[j,t]/p.area[j,t]
#N.area[j,t] ~ dunif(n.area[j,t],500)
#N.area.round[j,t]<-round(N.area[j,t])
}}
tau.e <- 1/(sd.e*sd.e)
sd.e ~ dunif(0,10)
tau.th[1:J, 1:J] ~ dwish(B[,], v)
sigma2[1:J, 1:J] <- inverse(tau.th[1:J,1:J])
v<-J
for(j1 in 1:J){
B[j1,j1]<-1
for(j2 in 1:j1-1){
B[j1,j2]<-0
}
for(j2 in j1+1:J){
B[j1,j2]<-0
}
for(j2 in 1:J){
pd[j1,j2] <- step(sigma2[j1,j2])
}}
for(i in 1:n+nzeros){
for(t in 1:T-1){
logit(phi[i,t]) <- logit(mu.phi) + e.phi[z.cut[i],t]
}
gamma[i,1]~dunif(0,1)
for(t in 2:T){
logit(gamma[i,t]) <-logit(mu.gamma) + e.gamma[z.cut[i],t]
}}
for(c in 1:C){
for (t in 1:T-1){
e.phi[c,t] ~ dnorm(0, tau.phi)
prob.phi[c,t]<-step(e.phi[c,t])
}
for(t in 2:T){
e.gamma[c,t] ~ dnorm(0, tau.gamma)
}}
mu.gamma ~ dunif(0,1)
tau.gamma <- 1/(sd.gamma*sd.gamma)
sd.gamma ~ dunif(0,10)
mu.phi ~ dunif(0,1)
tau.phi <- 1/(sd.phi*sd.phi)
sd.phi ~ dunif(0,10)
# Clustering
for(i in 1:n+nzeros){
z[i] ~ dcat(w[])
z.cut[i]<-cut(z[i])
for(c in 1:C){
post.prob[i,c] <- equals(z[i],c)}}
for(c in 1:C){
w[c] <- ww[c]/sum(ww[])
r[c] ~ dbeta(1,a)
NC[c] <- sum(post.prob[,c])
cl[c] <- step(NC[c]-1)}
ww[1] <- r[1]
for(c in 2:C){
ww[c] <- r[c]*(1-r[c-1])*ww[c-1]/r[c-1]}
CLUSTERS <- sum(cl[])
a<- 1 # controls prior smoothing
#derived parameters
for(i in 1:n+nzeros){
for(c in 1:C){
recruit[i,c,1]<-x[i,1]*equals(z[i],c)
for(t in 2:T){
recruit[i,c,t]<-(1-x[i,t-1])*(x[i,t])*equals(z[i],c)
death[i,c,t] <- x[i,t-1]*(1-x[i,t])*equals(z[i],c)
alive[i,c,t]<-x[i,t]*equals(z[i],c)
}}}
for(c in 1:C){
for(t in 2:T){
R[c,t]<-sum(recruit[1:n+nzeros,c,t])
D[c,t]<-sum(death[1:n+nzeros,c,t])
N[c,t]<-sum(alive[1:n+nzeros,c,t])
Rpc[c,t]<-R[c,t]/N[c,t]
recruitment.clus[c,t]<- exp(logit(mu.gamma) + e.gamma[c,t])/(1+exp(logit(mu.gamma) + e.gamma[c,t]))
}
for(t in 1:T-1){
survival.clus[c,t] <- exp(logit(mu.phi) + e.phi[c,t])/(1+exp(logit(mu.phi) + e.phi[c,t]))
}}
for(t in 2:T){
NN[t]<-sum(N[1:C, t])
DD[t]<-sum(D[1:C, t])
RR[t]<-sum(R[1:C, t])
RRpc[t]<-RR[t]/NN[t]
}
#partitioned p values
for (i in 1:n){
for(j in 1:J){
for (t in 1:T){
y.new[i,j,t]~dpois(mu.y[i,j,t])
# for predictive goodness of fit
p0[i,j,t] <- equals(y[i,j,t],0)
p0new[i,j,t]<-equals(y.new[i,j,t],0)
count[i,j,t]<-y[i,j,t]
countnew[i,j,t]<-y.new[i,j,t]
for (c in 1:C){
d.p0[i,j,t,c]<-p0[i,j,t]*equals(z[i],c)
d.p0new[i,j,t,c]<-p0new[i,j,t]*equals(z[i],c)
d.count[i,j,t,c]<-count[i,j,t]*equals(z[i],c)
d.countnew[i,j,t,c]<-countnew[i,j,t]*equals(z[i],c)
}}}}
for(j in 1:J){
for(t in 1:T){
for(c in 1:C){
fit1.p0[j,t,c] <- sum(d.p0[1:n,j,t,c])
fit.new1.p0[j,t,c] <- sum(d.p0new[1:n,j,t,c])
fit1.count[j,t,c] <- sum(d.count[1:n,j,t,c])
fit.new1.count[j,t,c] <- sum(d.countnew[1:n,j,t,c])
}}}
for(t in 1:T){
for(c in 1:C){
fit2.p0[t,c] <- sum(fit1.p0[1:J,t,c])
fit.new2.p0[t,c] <- sum(fit.new1.p0[1:J,t,c])
fit2.count[t,c] <- sum(fit1.count[1:J,t,c])
fit.new2.count[t,c] <- sum(fit.new1.count[1:J,t,c])
}}
for(c in 1:C){
DF.p0[c] <- sum(fit2.p0[,c])
DF.new.p0[c] <- sum(fit.new2.p0[,c])
pvalue.p0[c] <- step(DF.new.p0[c] - DF.p0[c])
DF.count[c] <- sum(fit2.count[,c])
DF.new.count[c] <- sum(fit.new2.count[,c])
pvalue.count[c] <- step(DF.new.count[c] - DF.count[c])
}
DFoverall.p0<-sum(DF.p0[1:C])
DFoverall.new.p0<-sum(DF.new.p0[1:C])
pvalueoverall.p0<-step(DFoverall.new.p0-DFoverall.p0)
DFoverall.count<-sum(DF.count[1:C])
DFoverall.new.count<-sum(DF.new.count[1:C])
pvalueoverall.count<-step(DFoverall.new.count-DFoverall.count)
}
输入数据(我在 R 中进行了所有数据格式化,因此我将其编译为一个列表,然后使用 r2winbugs 以 S/R 格式导出它,所以它应该可以在 WinBUGS 中工作):
bugs.data <- list(y = y,
x = x.data,
J = dim(y)[2],
T = dim(y)[3],
n = dim(y)[1]-100,
nzeros = 100,
C = 10,
n.area = n.area)
数据结构(我在实际的文本文件中手动确认了这些):
> dim(y)
[1] 727 7 41
> dim(x.data)
[1] 727 41
> dim(n.area)
[1] 7 41
> J
[1] 7
> T
[1] 41
> n
[1] 627
最后是陷阱错误:
index out of range
HostWindows.AppendInt [000002B7H]
.d ARRAY 12 OF CHAR "080'069'784'"
.i INTEGER 12
.j INTEGER 1971081876
.len INTEGER 1971172856
.n INTEGER 1
.s ARRAY 256 OF CHAR " Allocated Memory: " ...
.useSeparators BOOLEAN TRUE
HostWindows.UpdateInfo [000048B8H]
.res INTEGER 2291852
.sstr ARRAY 256 OF SHORTCHAR "l˜"" ...
.str ARRAY 256 OF CHAR " Allocated Memory: " ...
HostWindows.Idle [00004A73H]
.focus BOOLEAN FALSE
.tick Controllers.TickMsg Fields
.w HostWindows.Window NIL
HostMenus.TimerTick [00003422H]
.lParam INTEGER 0
.ops Controllers.PollOpsMsg Fields
.wParam INTEGER 1
.wnd INTEGER 459558
Kernel.Try [00003A61H]
.a INTEGER 459558
.b INTEGER 1
.c INTEGER 0
.h PROCEDURE HostMenus.TimerTick
HostMenus.ApplWinHandler [00003841H]
.Proc PROCEDURE NIL
.hit BOOLEAN Undefined117
.lParam INTEGER 0
.message INTEGER 275
.res INTEGER 2292416
.s ARRAY 256 OF SHORTCHAR 9AX, 4X, 0CX ...
.w INTEGER 0
.wParam INTEGER 1
.wnd INTEGER 459558
<system> (pc=757CC4E6H, fp=0022FB3CH)
<system> (pc=757CC5E6H, fp=0022FBB4H)
<system> (pc=757CCC18H, fp=0022FC14H)
<system> (pc=757C2E40H, fp=0022FC24H)
HostMenus.Loop [00003BDEH]
.done BOOLEAN FALSE
.f SET {0..5}
.n INTEGER 10
.res INTEGER 0
.w HostWindows.Window NIL
Kernel.Start [00002B8CH]
.code PROCEDURE HostMenus.Loop
非常感谢您的帮助!-乔什