Dear all
Thanks to all that responded.
Here is a summary of the responses I got to the question:
>I wonder whether anyone has come up with a trick to get around the
>limitation in BUGS that dmulti(p[], N) can't handle the situation where N is
>a stochastic node.
Three approaches were suggested by people responding :
1. use the zeroes trick for an arbitratry distribution (Andrew Millard)
2. use a constructive definition of a multinomial in terms of binomials
(Andrew Millard, David Elston, Finn Krogstad)
3. saying that the dmulti(p[],N) vector is the sum of N independent
vectors all of which are distributed as dmulti(p[],1) means that the problem is one of summing a set where the size is a random quantity
(Mark Brewer, Andy Lynch)
Before going into the more detailed suggestions for 2 and 3 it should be noted that David Spiegelhalters initial respons to my question whether it was any plans to allow N in multinomial models to be a stochastic node was:
Sorry, it's not allowed due to the single-site updating method -all the
n's would have to be simultaneously updated. Not sure if there is some
trick by expressing the multinomial as a normalised set of Poissons
- I think there might be. It might be worth asking on the list for
some clever person (although it has been asked before, if you check
the archives).
David
Ranta Jukka elaborated a bit on the problem:
The problem is the single site updating
algorithm in WinBUGS. If you have both N and the multinomial variable x_1,...,
x_k as random, there's no way of updating correctly without doing it
simultaneously at least for two components, since the x_i must always sum up to
N. Prospects look bad as long as we are not allowed to write more detailed
specification about which variables to update at each step and from which
proposal densities (as in the general Metropolis-Hastings algorithm). I think
the options are now: either to change your model to something else that is
reasonably sensible in your application, or to write your own sampler.
Ranta Jukka
I am not sure how this problem affects the suggested solutions using approach 2 and 3, (any comments will be of interest,) but here are the detailed suggestions.
Approach 2.
Finn Krogstad suggested the following (with a minor change in the code by me).
One trick is to represent the multinomial as a set of binomials which
don't need a constant order parameter.
model{
for (i in 1:I){
p[i]<-e[i]/sum(e[])
e[i]~dgamma(.001,.001)
for (j in 1:J){
y[j,i] ~ dbin(p[i], N)
}}
N~dpois(lambda)
lambda~dgamma(0.001,0.001)
}
#data
list(y=structure(.Data = c(53, NA, 7, 53, 22, NA), .Dim = c(2, 3)),I=3, J=2)
#list
list(e=c(1,1,1), lambda=82)
Approach 3.
Andy Lynch described it this way:
I don't know how efficient it is, as I have only used it for small problems,
but by saying that the dmulti(p[],N) vector is the sum of N independent
vectors all of which are distributed as dmulti(p[],1) means that the problem
is one of summing a set where the size is a random quantity.
This can be done using the indicator function trick given in the bugs user
manual. This was possible for me because my N had an upper bound, If yours
doesn't then I guess that this is not an option.
In any case, my BUGS knowledge is severely limited, so I am sure that there
will be better solutions offered by the great and the good
Andy Lynch
Mark Brewer gave a code example:
Just something off the top of my head - instead of having, say,
model{ y[1:K] ~ dmulti(p[],N)}
list(N=2,K=2,p=c(0.3,0.7))
which clearly has N as fixed, you might do something
like
model{
for(i in 1:bigN){
ytt[i,1:K] ~ dmulti(p[],1)
oneorzero[i] <- step(N-i)
for(k in 1:K){
yt[i,k] <- ytt[i,k] * oneorzero[i]
}
}
for(k in 1:K){
y[k] <- sum(yt[1:bigN,k])
}
N ~ dcat(alpha[])
alpha[1] <- 0.5
alpha[2] <- 0.3
alpha[3] <- 0.2
}
list(K=2,p=c(0.3,0.7),bigN=10)
which is probably hideously inefficient, but seems
to work for this trivial example. The idea is to make
bigN larger than N will ever be, and to create a
vector to sample from dmulti(p[],1) bigN times. You
only want N of these at any one time, and since
the sum(yt[1:N,k]) construct didn't seem to like a
variable N either, we use the step function to
obtain the value 1 for when i <= N and 0 otherwise,
form the product of this 0-1 vector and the sampled
bigN y's, and finally sum the yt's to get a sample
from a dmulti(p[],N).
A long way around, but it might help.
Best wishes,
Mark
Regards
Audun Stien
Audun Stien
CEH Banchory,
Hill of Brathens, Glassel, Banchory,
Kincardineshire,
AB31 4BW,
Scotland, U.K.
Telephone : (0044) (0) 1330 826336
Fax : (0044) (0) 1330 823303
E.mail: [log in to unmask]
-------------------------------------------------------------------
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
|