如何代码WinBUGS软件矩阵?(How to code a matrix in WinBUGS?)

2019-10-22 11:22发布

我想代码的4种元素的2X2矩阵西格玛。 不知道如何在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))

Answer 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[,])
}


文章来源: How to code a matrix in WinBUGS?