JY Barnagaud
2012-03-30
Hi forum;
I heard about that JAGS can't work out a Poisson with a parameter = 0
which is quite annoying in zip models (like the one partially quoted below)
N ~ dpois(esp)
esp <- w * lambda
w ~ dbern(psi)
I've had a "runtime error" while trying to update that model with Jags 3.2.0
under R -> could that poisson(0) be the reason? Or should I look at an issue
somewhere else? The 1000 adaptation interations run well, and then the first
8% of MCMC iterations, before R fails.
many thanks for advice!
bugbird
Martyn Plummer
2012-03-30
JAGS 3.2.0 can handle a zero mean parameter. If that were the problem it would
happen almost immediately.
You haven't given enough information to help anyone solve your problem. At the
very least you need to give the model in full and also the full error message
(It doesn't just say "RUNTIME ERROR" does it?).
JY Barnagaud
2012-03-30
Hi Martyn, thanks for your answer
Well you have answered my specific question about the P(0) issue - so it's
something else.
I did not keep the error message (would I have wanted I guess than even a
screen print wouldn't have been possible as the whole computer was frozen), it
was something like that the code ended in an inadequate way - in an
independent window, not produced inside the R device (which subsequently
failed).
I guess it might be due either to the model code itself or bad inits choice.
some slight simplification of the model is enough to fix it but here it is
anyway in case something is obvious.
Thanks again;
bugbird
model{
for (g in 1:numori){
mlb.mean~dunif(0,1)
mu.lam.beta<-log(mlb.mean)-log(1-mlb.mean)
mlb2.mean~dunif(0,1)
mu.lam.beta2<-log(mlb2.mean)-log(1-mlb2.mean)
} # g
for (h in 1:numori){ # loop over numori origins
tau.lam.beta~dgamma(0.1,0.1)
tau.lam.beta2~dgamma(0.1,0.1)
} # h
for (g in 1:numori){
mpb.mean~dunif(0,1)
mu.psi.beta<-log(mpb.mean)-log(1-mpb.mean)
mpb2.mean~dunif(0,1)
mu.psi.beta2<-log(mpb2.mean)-log(1-mpb2.mean)
} # g
for (h in 1:numori){
tau.psi.beta~dgamma(0.1,0.1)
tau.psi.beta2~dgamma(0.1,0.1)
} # h
for (i in 1:S){ # loop over S species
for (f in 1:nzone){
psi.alpha~dnorm(mu.psi.int_,tau.psi.int_)
} # f
mu.psi.int_~dnorm(0,0.001)
tau.psi.int_<- 1/(sigma.psi.int_*sigma.psi.int_)
sigma.psi.int_~dunif(0,100)
psi.beta_ ~ dnorm(mu.psi.beta[ST_],tau.psi.beta[ST_])
psi.beta2_ ~ dnorm(mu.psi.beta2[ST_],tau.psi.beta2[ST_])
for (f in 1:nzone){
lam.alpha~dnorm(mu.lam.int_,tau.lam.int_)
}
mu.lam.int_~dnorm(0,0.001)
tau.lam.int_<- 1/(sigma.lam.int_*sigma.lam.int_)
sigma.lam.int_~dunif(0,100)
lam.beta_ ~ dnorm(mu.lam.beta[ST_],tau.lam.beta[ST_])
lam.beta2_ ~ dnorm(mu.lam.beta2[ST_],tau.lam.beta2[ST_])
det.alpha_ ~ dnorm(0,0.001)
det.beta1_ ~ dnorm(0,0.001)
det.beta2_ ~ dnorm(0,0.001)
det.beta3_ ~ dnorm(0,0.001)
for (j in 1:R){ # loop over R sites
w ~ dbern(psi)
logit(psi)<-psi.alpha[i,ZON]+psi.beta_NATFOR500+psi.beta2_SHA500
N ~ dpois(esp)
esp <- w * lambda
log(lambda)<-lam.alpha[i,ZON]+lam.beta_NATFOR500+lam.beta2_SHA500
for (k in 1:K){
obs ~ dbin(p,N)
p <- theta*w
logit(theta) <- det.alpha_ + det.beta1_ * HOUR + det.beta2_ * VH + det.beta3_
* VH2
} # k
} # j
} # i
for (i in 1:S){
for (j in 1:R){
for (k in 1:K){
eval<-N*p
sd.resi<-sqrt(eval*(1-p))+0.5
E<-(obs-eval)/sd.resi
E2<-pow(E,2)
obs.new~dbin(p,N)
E.new<-(obs.new-eval)/sd.resi
E2.new<-pow(E.new,2)
} # k
} # j
} # i
fit<-sum(E2)
fit.new<-sum(E2.new)
} # model
Martyn Plummer
2012-03-30
Are you saying that JAGS locked up your whole computer? Because that sounds
bad.
JY Barnagaud
2012-03-30
Basically JAGS made R fail but it seems that it was more the failure of R than
JAGS itself which caused this temporary freeze. Might be due to the large
amount of memory used by R at the time of the failure due to the size of the
model & of the data (?)
(BTW my computer runs under Windows XP, 3.20 Ghz, dual core, 2.25Go RAM)
Martyn Plummer
2012-03-30
You might have just run out of memory. Try running the same model on a more
recent computer. If you have a really big problem, you need a 64-bit operating
system. On 32-bit there is a hard memory limit of 2Gb.
JY Barnagaud
2012-03-30
So hopefully I'll find a recent computer with 64b system somewhere and follow
your advice !
One very last question: I'm new to Jags so sorry if that's something trivial:
should I consider the n.adapt as equivalent to a burnin period (like
specifying the burn in period in OpenBUGS)? As far as I understood the answer
is "no" but not sure about it neither why.
Thanks a lot for your rapid help anyway.
Jack Tanner
2012-03-30
There's a discussion of adaptation in the JAGS manual. Please have a look and
see if that answers your question.
JY Barnagaud
2012-04-02
About adaptation: not too clear to me, it seems that it's not a MCMC run, but
I can't figure out what it is then.About running model: crashes also with 64 bits (update of the mode for 10000
iterations works fine, then monitoring parameters for 10000 aditional
iterations makes R fail with a "runtime error"). I'm monitoring
Martyn Plummer
2012-04-02
Let's do the math. You need 8 bytes to store a double-precision float, so the
minimum memory needed for your 10000 iterations is
> (917 * 15 * 3 + 917 * 15 + 6 * 15 ) * 10000 * 8 [1] 4408800000
or 4 gigabytes. Due to the way memory allocation works, you may need up to
twice that amount. So yes, you are going to run out of memory quite quickly.
Tim Handley
2012-04-02
Regarding adaptation:
There are several different ways of generating samples with MCMC. JAGS tries
to be smart about this, and choose the method that is most efficient for your
particular model. Each of these methods has semi-hidden parameters that you
won't see unless you're a programmer or a statistician. These are not model
parameters, they do not fit your data. Rather, they are meta-sampling-
parameters which can be tweaked to optimize the efficiency of sample
generation. During the adaptation period, JAGS is drawing samples, but it is
also tweaking the meta-sampling-parameters to help you get the best possible
performance.
In general, you want to run JAGS in two steps. First, call jags.model() with
n.adapt=100 to 1000. This allows JAGS to optimize the meta-parameters. You
cannot store the samples generated here, but that's fine. Once the sampler is
optimized, you then call jags.samples(). This second function uses the
optimized sampler to collect samples of your model parameters.
As an example, run jags.model with n.adapt=1. Then run jags.samples() with
n.iter=10000 and n.thin=100. Time the call to jags.samples(), and look at the
within-chain autocorrelation. Then do the same thing but with n.adapt=1000.
The second call to jags.camples() will likely be faster and more efficient
(lower autocorrelation) because you've given JAGS a lot of time (1000
iterations) in which to optimize the efficiency of the sampling.
JY Barnagaud
2012-04-03
that's exactly the kind of explanation I'm able to understand, thanks a lot!