Code Equation of Ellipse in WinBUGS

Question

I was looking for some help to code an equation of ellipse within WinBUGS. I need to form a Bivariate ellipse using p1's in my data. I tried to use the equation as (X-mu)'sigmainverse(X-mu), where X is the Bivariate Normal Variable, mu is the vector of means and sigmainverse is th e inverse of the var-covariance matrix. In my example p1's are Bivariate Normal Variables with mean gamma and inverse sigma2 matrix. Within double quotes is what I did but it doesnot work. Heres the WinBUGS code:

    model
    {      
     for (j in 1 : Nf)

  { 
  p1[j, 1:2 ] ~ dmnorm(gamma[1:2 ], T[1:2 ,1:2 ]) 

# gamma is the MVN mean or mean of logit (p)
#T is the precision matrix inverse sigma of MVN or logit(p)
# precision equals reciprocal of variance
# precision matrix is the matrix inverse of the covariance matrix

  for (i in 1:2)
  {
 logit(p[j,i])<-p1[j,i]

 Y[j,i] ~ dbin(p[j,i],n)

 wp[j,i] <- p[j,i]*dbw[j,i] 
 }

 sumwp[j] <- sum(wp[j, ])

 #X_mu[j,1:2]<-p1[j,1:2]-gamma[1:2]**

 #ell[j]<-((t(p1[j,1:2]-gamma[1:2]))*T[1:2,1:2]*(p1[j.1:2]-gamma[1:2]))**

    X_mu[j,1]<-p1[j,1]-gamma[1]
    X_mu[j,2]<-p1[j,2]-gamma[2]

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

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


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

       #ell[j]<-((t(p1[j,1:2]-gamma[1:2]))*T  
  }  

 # Hyper-priors:  

 gamma[1:2] ~ dmnorm(mn[1:2],prec[1:2 ,1:2])

  T[1:2 ,1:2] ~ dwish(R[1:2 ,1:2], 2)

 sigma2[1:2, 1:2] <-inverse(T[,])
  #sigma2 is the covariance matrix

 rho <- sigma2[1,2]/sqrt(sigma2[1,1]*sigma2[2,2])
#rho is the correlation matrix



 }

  expit[i]<-exp(gamma[i])/(1+exp(gamma[i]))
 }


 # 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)),

       dbw=structure(.Data=c(0.25,0.25,
       0.25,0.25,
      0.25,0.25,
      0.25,0.25,
       0.25,0.25,
       0.25,0.25,
      0.25,0.25,
      0.25,0.25,
       0.25,0.25,
     0.25,0.25,
      0.25,0.25,
      0.25,0.25,
       0.25,0.25,
    0.25,0.25,
      0.25,0.25,
      0.25,0.25,
       0.25,0.25,
      0.25,0.25,
      0.25,0.25,
       0.25,0.25
        ),.Dim=c(20,2))
       )

Answer 1

The * operator won't multiply matrices and vectors, just scalars. Unfortunately there's no general matrix product function in WinBUGS. Instead you could use two calls to the "inprod" function (or the faster "inprod2") to take the inner product of each row of T with the X - mu, giving a new (temporary) vector node. Then use another inprod to take the inner product of that vector with the transposed (X - mu), giving your ell[j]. Or if speed is a concern, just write the inner product out by hand, according to some reports this may be faster.