Как закодировать матрицу в WinBUGS?

Я пытаюсь закодировать сигму матрицы 2X2 с 4 элементами. Не уверен, как кодировать в WINBUGS. Моя цель — получить апостериорные p, их средние и отклонения и создать область эллипса, покрытую двумя апостериорными p. Вот мой код ниже:

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

winbugs winbugs14
person user1560215    schedule 31.03.2015    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Вы должны указать каждый элемент по очереди. Вы можете использовать функцию inverse (а не solve) для инвертирования матрицы.

model{
  sigma[1,1]<-v1
  sigma[1,2]<-v12
  sigma[2,1]<-v21
  sigma[2,2]<-v2
  sigmaInv[1:2, 1:2] <- inverse(sigma[,])
}
person guyabel    schedule 01.04.2015
comment
Я изменил свой код выше, но продолжаю получать несколько определений ошибки node ell[1,2] - person user1560215; 01.04.2015
comment
Может быть inprod2... возможно, вы хотите inprod, или может быть, у вас есть ell в ваших данных? Если это последнее, вы определяете его дважды, отсюда и сообщение об ошибке. - person guyabel; 02.04.2015
comment
Пожалуйста, посмотрите мой код выше, я его отредактировал, но теперь я получаю сообщение об ошибке множественных определений узла sigmaInv[1,1] - person user1560215; 03.04.2015
comment
Перепробовал несколько решений, ничего не работает. Его несколько определений узла для ell или sigmaInv. Пожалуйста помоги. - person user1560215; 04.04.2015