Replies to Andrew Millard <[log in to unmask]> please
Ken Rice
BUGS list manager
------------- Begin Forwarded Message -------------
Dear BUGSers
I am attempting to fit a multivariate normal with missing data. According
to the manual this is not possible, except by using a series of
conditional marginal normals. I think the model below using the ones
trick does the same, and I have attempted to model the example mentioned
way-back on this list by Simon Jackman, with severe missingness, and an
analytical solution for (at least) one of the posteriors.
This model seems to produce better results than Jackman's, but still is
someway off approaching the results in Tanner's Tools for Statistical
Inference (2nd ed. 1993 p64). The chains show quite a lot of
autocorrelation. Has anyone suggestions for improving the mixing and
convergence? I suspect that tau could be better parameterised.
Also I wish to actually apply this method to a trivariate dataset. Are
there any simple ways for setting conditional priors on the elements of
tau to ensure that it is positive-definite? Currently I would set the
diagonal elements positive and then condition the other elements using the
symmetry and requirement that the determinant is positive. I'm not a
matrix-algebra person so I'm not even sure if this is sufficient.
Any comments gratefully received,
Andrew
=========================================================================
Dr. Andrew Millard [log in to unmask]
Department of Archaeology, University of Durham, Tel: +44 191 374 4757
South Road, Durham. DH1 3LE. United Kingdom. Fax: +44 191 374 3619
http://www.dur.ac.uk/a.r.millard/
=========================================================================
---- begin WinBUGS code ----
######################################################################
## fit bivariate normal model to the data of Tanner, Table 5.1
## originally implement in WinBUGS by Simon Jackman
## http://tamarama.Stanford.edu/simon/mcmc/tanner51.winbug
## http://www.jiscmail.ac.uk/cgi-bin/wa.exe?A2=ind9811&L=bugs&O=D&F=&S=&P=967
##
## completely rewritten to use multivariate ones trick
## Andrew Millard, Dept Archaeology, Univ. of Durham July 2001
## [log in to unmask] http://www.durham.ac.uk/a.r.millard/
######################################################################
model{
PI <- 3.1415926535897932384626433832795
d <- 2 # dimensionality of MVN
norm.const <- -0.5*d*log(2*PI) # log of normalising constant for MVN
for (i in 1:N){
for (j in 1:2){
#for this problem mu=c(0,0) is given. Priors are
possible here.
mu[i,j] <- 0
}
# use ones trick to give multivariate normal on x
mvn.ones[i] <- 1
mvn.ones[i] ~ dbern( mvn.p[i] )
log(mvn.p[i]) <- norm.const + 0.5*logdettau -
0.5*matprod[i]
# calculate matrix multiplication term in MVN density
# could use loops for generality in higher dimensions
for (j in 1:2){
y[i,j] <- x[i,j] - mu[i,j]
}
matprod[i] <- y[i,1]*y[i,1]*tau[1,1]
+ 2*y[i,1]*y[i,2]*tau[1,2]
+ y[i,2]*y[i,2]*tau[2,2]
# need vague prior-prior on x
for (j in 1:2){
x[i,j] ~ dnorm(0,0.0001)
}
}
# need vague prior on tau
# use nodes aa, ab ac to ensure updates in correct order and
# constraints are obeyed
tau[1,1] <- aa
aa ~ dunif(0,100)
tau[2,2] <- ab
ab ~ dunif(0,100)
limit <- sqrt(tau[1,1]* tau[2,2])
neglimit <- -limit
tau[1,2] <- ac
ac ~ dunif(neglimit,limit) # det(tau) is positive so
ac^2<aa*ab
tau[2,1] <- tau[1,2] # tau is symmetric
logdettau <- logdet(tau[,]) # saves calculating it N times
# nodes for monitoring values of interest
Var1 <- inverse(tau[,],1,1)
Var2 <- inverse(tau[,],2,2)
Cov12 <- inverse(tau[,],1,2)
# Tanner gives plots of posteriors for correlation coefficient
rho <- Cov12/sqrt(Var1*Var2)
# and smaller eigenvalue of the covariance matrix
lambda <- -0.5*(-(Var1+Var2) +
sqrt(pow((Var1+Var2),2) -
4*(Var1*Var2-pow(Cov12,2))))
}
Data
list(x = structure(.Data=c(1, 1,
1, -1,
-1, 1,
-1, -1,
2, NA,
2, NA,
-2, NA,
-2, NA,
NA, 2,
NA, 2,
NA, -2,
NA, -2), .Dim = c(12,2)),
N=12)
Inits
list(x = structure(.Data=c(NA, NA,
NA, NA,
NA, NA,
NA, NA,
NA, 0,
NA, 0,
NA, 0,
NA, 0,
0, NA,
0, NA,
0, NA,
0, NA), .Dim = c(12,2)),
aa = 1, ab = 1, ac=0.5)
------------- End Forwarded Message -------------
-------------------------------------------------------------------
To mail the BUGS list, mail to [log in to unmask]
You can search old messages at www.jiscmail.ac.uk/lists/bugs.html
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
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]
|