Sorry it's taken me so long to post a summary of replies. I plead
weight of other projects + illness + losing use of a Windows computer
(I normally run Linux).
It turned out that I had 2 problems. Firstly my code was wrong, and
secondly there is a bug in WinBUGS version 2.01beta (Nov 2004), which
I was using.
The error message from WinBUGS2.01 was "UpdaterGamma.Updater", while
WinBUGS1.4 gives the error "Cannot sample from interval censored gamma
full conditional w[6]". This is much more helpful, because there were
n=5 observed taxis, so clearly, as soon as it tries to sample a
failure time from one of the unobserved ones, it has difficulty. I
still don't understand why.
But Mike Thompson fixed the problem, with a very clever piece of
coding (lines #1, #2 and #3 in the code below). It retains the upper
limit M on the population size, with mythical taxis N+1,...,M, but
uses a cut function to ensure that these taxis do not influence the
likelihood.
Mike's code executes without problems under WinBUGS1.4, but fails
under 2.01 with the same "UpdaterGamma.Updater" error as before. I've
reported this as a bug (but had no acknowledgement: is this normal?)
Thanks also to Craig Aumann, Finn Krogstad, Mary Kynn, Charles
Linveld, and Bob O'Hara. Bob pointed out that the results are likely
to be very sensitive to the assumed form of the distribution of failure
times. I'll investigate that further.
Ted Catchpole.
# -------------
model heterotaxi
{
# Likelihood
# lam.uo = lambda.unobserved
# lam.ne = lambda.nonexistent
for (i in 1:n){ # loop over observed taxis
for (j in offset[i]:(offset[i+1]-1)){
t[j] ~ dexp(lambda[i]) # observed times for taxi i
}
u[i] ~ dexp(lambda[i])I(x[i],) # censored time for taxi i
}
for (i in 1:M){
lam.ne[i] <- cut(lambda[i+n]) #1
lam.uo[i] <- step(N-i)*lambda[i+n] + (1-step(N-i))*lam.ne[i] #2
w[i] ~ dexp(lam.uo[i])I(T,1.0E3) #3
}
for (i in 1:(n+M)){
lambda[i] ~ dgamma(alpha,beta)
}
# Prior
alpha ~ dgamma(0.01,0.01)
beta ~ dgamma(0.1,0.1)
N ~ dpois(mu)I(1,M)
mu ~ dunif(0,M)
}
# Initial values
list(alpha=1, beta=10, N=10, mu=10,
u = c(10,10,10,10,10),
w = c(
40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40))
# Data
list(n=5, T=30, offset=c(1,7,14,21,25,32), M=200,
t=c(0.07, 4.60, 3.76, 7.88, 1.16, 5.43, 2.13, 5.58, 6.58, 2.07, 3.38,
8.16, 0.36, 3.23, 1.73, 3.98, 6.98, 6.99, 5.55, 0.85, 0.05, 13.97,
11.77, 3.34, 0.26, 15.15, 0.46, 5.92, 0.72, 0.77, 1.24),
x=c(7.10, 1.74, 0.69, 0.87, 5.49))
# -------------
On 14/04/05 16:26, ecatchpole wrote,:
> The following code compiles fine, but then appears to halt almost
> immediately with the message "UpdaterGamma.Updater"appearing.
> Can anyone suggest what is going wrong, please?
>
> The problem is this. Taxis fail at a heterogeneous Poisson rate lambda
> --- i.e., a separate lambda for each taxi, from some population, say
> lambda ~ gamma(alpha,beta). They are repaired instantly and sent back
> out on the road. After time T we have one or more failure times, plus
> one censored failure time, for each of n taxis seen by the repair shop,
> plus an unknown number N that have not failed by time T. We wish to
> estimate N.
>
> Ted Catchpole
>
> # -------------
> model heterotaxi
> {
> # Likelihood
> for (i in 1:n){ # loop over observed taxis
> for (j in offset[i]:(offset[i+1]-1)){
> t[j] ~ dexp(lambda[i]) # observed times for taxi i
> }
> u[i] ~ dexp(lambda[i])I(x[i],) # censored time for taxi i
> }
> for (i in 1:M){ # M is a large dummy value, bigger
> than N
> c[i] <- T*step(N-i) # N is the population size
> w[i] ~ dexp(lambda[i+n])I(c[i],) # w[1],...,w[N] are the unobserved
> values
> # all censored at time T
> }
> for (i in 1:(n+M)){
> lambda[i] ~ dgamma(alpha,beta)
> }
>
> # Prior
> alpha ~ dgamma(1,5) # I've tried vague priors for alpha and beta too
> beta ~ dgamma(1,1) # --- it doesn't seem to make any difference
> N ~ dpois(mu)I(1,M) # N has a vague prior
> mu ~ dunif(0,M)
> }
>
> # Initial values
> list(alpha=0.2, beta=1, N=10, mu=10,
> u = c(10,10,10,10,10),
> w = c(
> 40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
> 40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
> 40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
> 40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
> 40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
> 40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
> 40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
> 40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
> 40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
> 40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,20))
>
> # Data
> list(n=5, T=30, offset=c(1,7,14,21,25,32), M=200,
> t=c(0.07, 4.60, 3.76, 7.88, 1.16, 5.43, 2.13, 5.58, 6.58, 2.07, 3.38,
> 8.16, 0.36, 3.23, 1.73, 3.98, 6.98, 6.99, 5.55, 0.85, 0.05, 13.97,
> 11.77, 3.34, 0.26, 15.15, 0.46, 5.92, 0.72, 0.77, 1.24),
> x=c(7.10, 1.74, 0.69, 0.87, 5.49))
>
> # -------------
>
--
Dr E.A. Catchpole
Visiting Fellow
Univ of New South Wales at ADFA, Canberra, Australia
and University of Kent, Canterbury, England
- www.ma.adfa.edu.au/~eac
- fax: +61 2 6268 8786
- ph: +61 2 6268 8895
-------------------------------------------------------------------
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
|
