Hello,
I hope this is an appropriate topic for this list. I've run into a very
perplexing problem: I have 2 WinBUGS model codes that do virtually the
same thing. Both codes have virtually the same structure; one was
written by a colleague and the other was copied from a paper by Meyer
and Millar (1999, Canadian Journal of Fisheries and Aquatic Science).
However, one model compiles and runs successfully while the other seems
to send the compiler into an infinite loop! I've played a lot of little
games with the offending model code--e.g., changing variable names,
eliminating tabs in the code, inserting spaces--but nothing has worked.
My guess is that a "reserved" name is being used as a variable, but I've
changed most of them and haven't had any luck. It's probably something
completely obvious, but I've been staring at it too long to see it. Any
suggestions would be greatly appreciated. Both codes and a dataset are
appended below.
Thanks for any and all help,
Buck Stockhausen
************************************************************************
*William T. Stockhausen, Ph.D. |email: [log in to unmask]*
*NRC Postdoctoral Fellow |voice: 508-289-2238 *
*Northeast Fisheries Science Center|fax: 508-457-2169 *
*National Marine Fisheries Service |web: http://www.nefsc.noaa.gov *
* | http://www.usglobec.org *
*__________________________________|___________________________________*
* mailing address: National Marine Fisheries Service *
* 166 Water Street *
* Woods Hole, MA 02543-1026 *
************************************************************************
The model that won't compile:
# version of Meyer and Millar 1999 CJFAS 56:1078-1086 model
model fred
{
# prior distribution of K: lognormal with 10%, 90% quantiles at 80 and 300
# K ~ dlnorm(5.042905,3.7603664)I(10,1000)
K ~ dunif(10,1000)
# prior distribution of r: ln with 10,90% quantiles at 0.13 and 0.48
# r ~ dlnorm(-1.38,3.845)I(0.01,1.2)
r ~ dunif(0.01,1.2)
# prior for q: inverse gamma as proper for p(q) ~ 1/q
invq ~ dgamma(0.001,0.001)I(0.1,1000)
qf <- 1/invq
# prior for sigma2: inv. gamma with 10, 90% q.s at 0.04 and 0.08
isigma2 ~ dgamma(4.0,0.01)
sigma2 <- 1/sigma2
# prior for tau2: inv. gamma with 10, 90% q.s at 0.05 and 0.15
itau2 ~ dgamma(2.0,0.1)
tau2 <- 1/itau2
# conditional prior for Ps
initP ~ dunif(0.01,1.5)
Pmedn[1] <- log(initP)
P[1] ~ dlnorm(Pmedn[1],isigma2)I(0.001,3.0)
for (i in 2:N) {
Pmedn[i] <- max(log(P[i-1]+r*P[i-1]*(1-P[i-1])-C[i-1]/K),0.0001)
P[i] ~ dlnorm(Pmedn[i],isigma2)I(0.0001,3.0)
}
for (i in 1:N) {
Imedn[i] <- log(qf*K*P[i])
Iobs[i] ~ dlnorm(Imedn[i],itau2)
Ires[i] <- Iobs[i]-qf*K*P[i]
}
# management parameters
MSY <- r*K/4
FMSY <- r/2
BMSY <- K/2
for (i in 1:N){
B[i] <- P[i]*K
F[i] <- C[i]/B[i]
Fratio[i] <- F[i]/FMSY
Bratio[i] <- B[i]/BMSY
}
}
# end model
The model that does compile:
# Bayesian Surplus Production Model
model bspvarb1
{
# Prior distributions
K ~ dunif(50.0,200.0)
r ~ dunif(0.01,1.99)
iqUSFall ~ dgamma(0.001,0.001)I(0.1,1000)
qUSFall <- 1/iqUSFall
isigma2 ~ dgamma(4.0,0.01)
sigma2 <- 1/isigma2
itau2USFall ~ dgamma(2.0,0.1)
tau2USFall <- 1/itau2USFall
initP ~ dunif(0.01,1.5)
# compute B as proportions of K each year
Pmean[1] <- log(initP)
P[1] ~ dlnorm(Pmean[1],isigma2)I(0.001,3)
for (i in 2:N){
Pmean[i] <- log(max(P[i-1]+r*P[i-1]*(1-P[i-1])-C[i-1]/K,0.0001))
P[i] ~ dlnorm(Pmean[i],isigma2)I(0.0001,3)
}
# indices
# US Fall - values for all years
for (i in 1:N){
ImeanUSFall[i] <- log(qUSFall*K*P[i])
Iobs[i] ~ dlnorm(ImeanUSFall[i],itau2USFall)
residUSFall[i] <- Iobs[i]-qUSFall*K*P[i]
}
# management parameters
MSY <- r*K/4
FMSY <- r/2
BMSY <- K/2
for (i in 1:N){
B[i] <- P[i]*K
F[i] <- C[i]/B[i]
Fratio[i] <- F[i]/FMSY
Bratio[i] <- B[i]/BMSY
}
}
# end model
A data set:
list(N=23,
C=c(15.9 ,25.7 ,28.5 ,23.7 ,25.0 ,33.3 ,28.2 ,19.7,17.5 ,19.3 ,21.6,23.1 ,22.5 ,22.5 ,23.6 ,29.1 ,14.4 ,13.2 ,28.4 ,34.6 ,37.5 ,25.9,25.3),
Iobs=c(61.89,78.98,55.59,44.61,56.89,38.27,33.84,36.13,41.95,36.63,36.33,38.82,34.32,37.64,34.01,32.16,26.88,36.61,30.07,30.75,23.36,22.36,21.91))
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|