Кодовое уравнение эллипса в WinBUGS

Question

Я искал некоторую помощь для кодирования уравнения эллипса в WinBUGS. Мне нужно сформировать двумерный эллипс, используя p1 в моих данных. Я попытался использовать уравнение как (X-mu)'sigmainverse(X-mu), где X - это двумерная нормальная переменная, mu - вектор средних значений, а sigmainverse - это инверсия матрицы вар-ковариации. В моем примере p1 - это двумерные нормальные переменные со средней гаммой и обратной матрицей sigma2. В двойных кавычках это то, что я сделал, но это не работает. Вот код WinBUGS:

    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

Оператор * не будет умножать матрицы и векторы, только скаляры. К сожалению, в WinBUGS нет общей функции продукта матрицы. Вместо этого вы могли бы использовать два вызова функции "inprod" (или более быстрый "inprod2"), чтобы получить внутреннее произведение каждой строки T с X - mu, давая новый (временный) векторный узел. Затем используйте другой inprod, чтобы взять внутреннее произведение этого вектора с транспонированным (X - mu), давая свой ell[j]. Или, если скорость имеет значение, просто напишите внутренний продукт вручную, согласно некоторым сообщениям, это может быть быстрее.