1

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))
       )
user1560215
  • 227
  • 1
  • 13

1 Answers1

3

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.

Chris Jackson
  • 871
  • 5
  • 7
  • Thanks for your suggestions.When I do the first step to sbtract mu's from X's which in my case would be p1[j:1:2]-gamma[1:2] I get an error 'expected multivariate node'. Please see my code above. I need some help with coding it. – user1560215 Sep 19 '13 at 19:28
  • I am really looking for some help with coding. I am trying to perform the operation elementwise but the inprod function doesnt seem to be working and is giving me incorrect results. Please see the modifications in my code – user1560215 Sep 22 '13 at 19:18
  • 1
    The definitions of T1 looks OK in the edited code above. Now you have a vector T1 and you need to take the inner product with the vector X_mu, so this should do it: ell[j]<-inprod(X_mu[j,],T1[j,]) – Chris Jackson Sep 26 '13 at 15:45