1

我正在寻找一些帮助来编写 WinBUGS 中的椭圆方程。我需要在我的数据中使用 p1 形成一个双变量椭圆。我尝试将方程用作 (X-mu)'sigmainverse(X-mu)，其中 X 是二元正态变量，mu 是均值向量，sigmainverse 是 var-covariance 矩阵的逆矩阵。在我的示例中，p1 是具有平均 gamma 和逆 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))
       )
4

1 回答 1

3

* 运算符不会乘以矩阵和向量，只会乘以标量。不幸的是，WinBUGS 中没有通用的矩阵乘积函数。相反，您可以使用对“inprod”函数（或更快的“inprod2”）的两次调用来获取 T 的每一行与 X - mu 的内积，从而给出一个新的（临时）向量节点。然后使用另一个inprod 将该向量与转置后的 (X - mu) 相乘，得到 ell[j]。或者如果速度是一个问题，只需手动写出内部产品，根据一些报道，这可能会更快。

于 2013-09-19T10:42:56.083 回答