How to code a matrix in WinBUGS?

问题:

I am trying to code the 2X2 matrix sigma with the 4 elements. Not sure how to code in WINBUGS. My goal is to get the posterior p's, their means and variances and create an ellipse region covered by the two posterior p's. Heres my code below:

model{
#likelihood

for(j in 1 : Nf){
p1[j, 1:2 ] ~ dmnorm(gamma[1:2], T[1:2 ,1:2])  
for (i in 1:2){
    logit(p[j,i]) <- p1[j,i]
    Y[j,i] ~ dbin(p[j,i],n) 
     }

 X_mu[j,1]<-p[j,1]-mean(p[,1])
 X_mu[j,2]<-p[j,2]-mean(p[,2])

v1<-sd(p[,1])*sd(p[,1])
v2<-sd(p[,2])*sd(p[,2])
v12<-(inprod(X_mu[j,1],X_mu[j,2]))/(sd(p[,1])*sd(p[,2]))

sigma[1,1]<-v1
sigma[1,2]<-v12
sigma[2,1]<-v12
sigma[2,2]<-v2
sigmaInv[1:2, 1:2] <- inverse(sigma[,]) 

T1[j,1]<-inprod(sigmaInv[1,],X_mu[j,1])

T1[j,2]<-inprod(sigmaInv[2,],X_mu[j,2])


ell[j,1]<-inprod(X_mu[j,1],T1[j,1])
ell[j,2]<-inprod(X_mu[j,2],T1[j,2])


}   

#priors
gamma[1:2] ~ dmnorm(mn[1:2],prec[1:2 ,1:2])
expit[1] <- exp(gamma[1])/(1+exp(gamma[1]))
expit[2] <- exp(gamma[2])/(1+exp(gamma[2]))
T[1:2 ,1:2] ~ dwish(R[1:2 ,1:2], 2)
sigma2[1:2, 1:2]  <- inverse(T[,])
rho  <-  sigma2[1,2]/sqrt(sigma2[1,1]*sigma2[2,2])
}


  # Data
 list(Nf =20, mn=c(-0.69, -1.06), n=60,
 prec = structure(.Data = c(.001, 0,
            0, .001),.Dim = c(2, 2)),
 R = structure(.Data = c(.001, 0,
         0, .001),.Dim = c(2, 2)),
 Y= structure(.Data=c(32,13,
         32,12,
         10,4,              
        28,11,                  
        10,5,                  
       25,10,
        4,1,
       16,5,
       28,10,
       21,7,
      19,9,
     18,12,
     31,12,
      13,3,
     10,4,
     18,7,
     3,2,
    27,5,
    8,1,
     8,4),.Dim = c(20, 2))

回答1:

You have to specify each element in turn. You can use the inverse function (rather than solve) to invert the matrix.

model{
  sigma[1,1]<-v1
  sigma[1,2]<-v12
  sigma[2,1]<-v21
  sigma[2,2]<-v2
  sigmaInv[1:2, 1:2] <- inverse(sigma[,])
}