Расшифровка ловушки - сообщение "Индекс вне диапазона" в WinBUGS
Мне было поручено перезапустить чужую модель с обновленными данными (я сделал обновление и форматирование данных). Я понимаю, как работает модель, но не написал ее, и она довольно длинная и подробная. Кроме того, я использую исключительно JAGS, и я впервые захожу в WinBUGS, поэтому интерфейс и сообщения об ошибках являются новыми для меня. Для компиляции модели требуется два дня (и она компилируется), но затем я получаю ошибку Trap, которая говорит, что Index Out of Range выходит за несколькими десятками строк, которые мне непонятны. Я просмотрел всю модель, циклы for и структуру данных и не вижу, где возникают проблемы с индексацией. И, к сожалению, в отличие от JAGS, WinBUGS не показывает, какой индекс находится вне диапазона.
Я был бы очень признателен за любую помощь в расшифровке сообщения об ошибке Trap, на случай, если там будет больше информации, чем я понимаю, чтобы указать местоположение несоответствия индекса. Я собираюсь опубликовать модель ниже, вместе с сообщением о ловушке, но я не вижу смысла в публикации фактических данных, поскольку никто не собирается тратить 2 дня в ожидании компиляции модели, чтобы воспроизвести ошибку, Вместо этого я просто опубликую структуру данных, чтобы измерения каждого объекта данных были очевидны.
Несколько соображений: я использую точную модель, которая была ранее использована (буквально скопирована и вставлена), поэтому я знаю, что она должна работать. Предыдущие файлы данных представлены в прямоугольном формате, а созданные и экспортированные файлы - в формате S. Насколько я понимаю, WinBUGS должен быть в состоянии справиться с обоими, но мне интересно, может ли это быть частью проблемы. Я запускаю модель напрямую из WinBUGS, загружаю модель и файлы данных вручную. Я пытался сделать это из R с r2winbugs ранее, но получал ошибки памяти и сбои.
Модель:
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
И, наконец, ошибка Trap:
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
Заранее большое спасибо за вашу помощь! -Josh