6

我的任务是使用更新的数据重新运行其他人的模型(我进行了数据更新和格式化)。我了解模型的工作原理但没有编写它,而且它很长而且很详细。最重要的是,我专门使用 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

非常感谢您的帮助!-乔什

4

0 回答 0