如何在WinBUGS中編碼矩陣？

我想用4個元素編碼2X2矩陣西格馬。不知道如何在WINBUGS中編碼。我的目標是獲得後驗p，它們的均值和方差，並創建一個由兩個後驗p覆蓋的橢圓區域。下面是我的代碼：如何在WinBUGS中編碼矩陣？

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

來源

2015-03-31 user1560215

您必須依次指定每個元素。您可以使用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-01 09:39:21 gjabel

我已經修改了我上面的代碼，但不斷收到節點ELL的錯誤多個定義[1,2] – user1560215 2015-04-01 19:58:11

可能是'inprod2' ......也許你想'inprod'，或者它可能是你有你的數據中有「ell」嗎？如果它是後者，則您將其定義兩次，因此會顯示錯誤消息。 – gjabel 2015-04-02 12:24:32

請參閱我上面的代碼，我已經編輯它，但現在我得到錯誤節點sigmaInv [1,1] – user1560215 2015-04-03 00:18:24

如何在WinBUGS中編碼矩陣？

回答

相關問題