1

我正在尝试用 4 个元素对 2X2 矩阵 sigma 进行编码。不知道如何在 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))
4

1 回答 1

1

您必须依次指定每个元素。您可以使用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[,])
}
于 2015-04-01T09:39:21.377 回答