Hi All,
Many thanks to those who responded to my question about the multinomial
probit model.
Below is the summary of the reponses.
1. It seems that the MNP is best implemented outside of WinBUGS.
Alternative packages/programs include:
a) Kosuke Imai and David van Dyk at Harvard have a (free) distributed
program that does Bayesian MNP. It's at
http://www.people.fas.harvard.edu/~kimai/research/BAMDA.html
b) Simon Jackman at Stanford implemented McCullock and Rossi(1994) in
gauss. Please see http://jackman.stanford.edu/mcmc
c) The new version of SAS (version 8.2) estimates MNP using GHK algorithm.
2. Andrew Millard provided the following code on Cholesky decomposition in
a MVN with missing data.
######################################################################
## fit bivariate normal model to the data of Tanner, Table 5.1
## originally 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 prior 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 Cholesky decomposition to ensure positive-definite and
remove
# need for conditional priors
# NB not clear (yet) what are the implied priors on tau, Var-Cov
matrix
# and correlations
# on and below the diagonal
for (i in 1:d) {
for (j in 1:i) {
A[i,j] ~ dunif(-100, 100) # could this be
improved?
for (k in 1:i){
sum.term[i,j,k] <- A[i,k]*A[j,k]
}
tau[i,j] <- sum(sum.term[i,j,1:i])
}
}
#above the diagonal
for (i in 1:d) {
for (j in i+1:d) {
A[i,j] <- 0
tau[i,j] <- tau[j,i]
}
}
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, d=2)
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)),
A=structure(.Data=c(1,NA,1,1),.Dim=c(2,2)))
--Yuanping
-------------------------------------------------------------------
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
|