Dear List,
There is a (well-known?) example of using data augmentation to deal with
an odd pattern of missing data in a bivariate normal setup; this example
appears at pp93ff of Tanner's _Tools for Statistical Inference_ (3rd ed)
and also appears in Tanner & Wong's 1987 JASA piece. In the JASA article
this example is attributed to comments by Murray on Dempster, Laird and
Rubin's EM article from 1977. So this is an example with quite a
history...
I was wondering if anyone had ever coded this example using Gibbs sampling
and in BUGS in particular? I'm having a lot of difficulty getting this
right in BUGS, though do a better job with my own Gibbs sampling routine
in S (sampling from conditional univariate Normals, and Wishart sampling).
In particular, given the pattern of missing data (see attached file), we
expect a bimodal posterior on rho (see Fig 5.4 Tanner 3rd ed, p96).
I'm a little unfamiliar with multivariate nodes and Wishart priors, and
I'd greatly appreciate any insights into how to implement this (simple?!)
example in BUGS.
Thanks -- Simon Jackman
==========================================================================
Simon Jackman, Asst Professor, ph: +1 (650) 723-4760
Dept of Political Science, fax: +1 (650) 723-1808
455 Serra Mall, Bldg 160, [log in to unmask]
Stanford University, http://tamarama.stanford.edu/simon
Stanford CA, 94305-2044, USA.
1 1
1 -1
-1 1
-1 -1
2 NA
2 NA
-2 NA
-2 NA
NA 2
NA 2
NA -2
NA -2
###############################################################################
##
## fit bivariate normal model to the data of Tanner, Table 5.1
##
## simon jackman, dept of political science, stanford univ, nov 1998
###############################################################################
model tanner51;
const
N=12,
K=2;
var
x[N,K], # data
mu[K],
Omega[K,K], # precision matrix, data
Sigma[K,K], # covariance matrix, data
R[K,K], # prior on Sigma (Wishart)
tau1, tau2,
rho;
data x in "tanner51.dat";
{
for (i in 1:N){
#x[i,] ~ dmnorm(mu[],Omega[,]);
mu[1] <- rho*sqrt(Sigma[1,1]/Sigma[2,2])*x[i,2];
mu[2] <- rho*sqrt(Sigma[2,2]/Sigma[1,1])*x[i,1];
x[i,1] ~ dnorm(mu[1],tau1);
x[i,2] ~ dnorm(mu[2],tau2);
}
tau1 <- 1.0/(Sigma[1,1]*(1-(rho*rho)));
tau2 <- 1.0/(Sigma[2,2]*(1-(rho*rho)));
Omega[,] ~ dwish(R[,],2);
Sigma[,] <- inverse(Omega[,]);
rho <- Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2]);
R[1,1] <- 1.0; R[1,2] <- 0.0;
R[2,1] <- 0.0; R[2,2] <- 1.0;
}
|