急！！求解关于winbugs软件怎么运用矩阵正态分布的问题

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

大家好，我这有一个关于如何用winbugs来构建多层贝叶斯模型的问题，是有关矩阵正态分布的构建的，该模型如下所示：
Yij=∑βik*Xjk+εij=xj’βi+εij
其中:
xj=（x1，x2，…xk）j’ j=1,2…J
βi=(β1,β2,…βk）i’ i=1,2…I
εij～N(0,δ2)
βi=Γ*Ζi+ζi
i=1,2…I
ζi～Nk(0,Σ)
Ζi=（z1，z2，…zm）i’
Γ为k﹡m 矩阵
Γ～Nk*m(w,W)(服从矩阵正态分布）也就是vec（Γ）～Nkm*1(vec（w）,W)
Σ-1～Wishart(v,V)


目前我是这样建的：
model;
{
   for( j in 1 : J ) {
      for( i in 1 : I ) {
         y[i , j] ~ dnorm(mu[i , j],tau)
      }
   }
   for( j in 1 : J ) {
      for( i in 1 : I ) {
         mu[i , j] <- inprod(x[j , 1:4],beta[1:4 , i])
      }
   }
   for( i in 1 : I ) {
      beta[1:4 , i] ~ dmnorm(mu[1:4 , i],R[1:4 , 1:4])
   }
   tau ~ dgamma(0.001,0.001)
   sigama <- 1 / sqrt(tau)
   for( i in 1 : I ) {
      mu[1:4 , i] <- inprod(gaba[1:4 , 1:3],z[1:3 , i])
   }
   R[1:4 , 1:4] ~ dwish(omega,4)
gaba[1:4 , 1:3]~ dmnorm(mean,pre)
}




这里面有个问题就是矩阵gaba[1:4,1:3]我在文中给它的先验分布是服从矩阵正态分布，可是没有相关的code，

很多资料中他们先将这个矩阵进行向量化算子，就是vec（gaba），然后再让他服从vec（gaba）~dmnorm（mean，pre），不知道这一步怎么处理，而且这边弄不好的话，下面的data load老是不对，

不知道哪位高手能帮我解决一下这个问题啊，谢谢啊，急求！

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Winbugs软件 winbugs WINBUG 正态分布 BUGS 软件 矩阵 正态分布 求解 winbugs

不知道将R软件调用winbugs的话是否可以做出来呢，怎么没有人回复呢？？？？哎。。。

1# azheforever
1、你要把你的模型的思想说清楚？2、你的模型要做什么？
然后其他朋友在读懂你的模型和程序的基础上才有可能帮你调试程序和解决问题？

我对这个多层bayes model没有研究过，所以清楚的程序和说明非常重要！

谢谢楼上，我来解释一下，我这个多层贝叶斯模型的构建是用于消费者选择决策的研究，
(1)首先构建一个回归：
Yij=∑βik*Xjk+εij=xj’βi+εij      εij～N(0,δ2)
其中βi=(β1,β2,…βk）i’ i=1,2…I    和 δ2是我们要求的参数，
(2)但是所能获得消费者个体数据较少，无法进行最小二乘估计，所以给βi和δ2加个先验分布，就是：
多重线性回归：βi=Γ*Ζi+ζi   i=1,2…I,   ζi～Nk(0,Σ)  Zi=(Ζi=（Z1，Z2，…Zm）i’
以及共轭先验：δ2~IG(a,b)
其中Zi，a，b都是已知
(3)由于上述先验中仍有未知参数Γ和Σ
所以对这些新的参数再加个先验分布，也就是
Γ～Nk*m(w,W)     
Σ-1～Wishart(v,V)
总而言之βi=(β1,β2,…βk）i’ i=1,2…I ，Γ，Σ-1，δ2
都是我们需要求的参数，
不知道有没有说清楚啊，我看winbugs好像对Γ～Nk*m(w,W)无法表示，外文文献是通过R软件的baysm package做的，不知道有没有哪位朋友帮忙编出来啊
，万分感谢啊

想法非常好，有一定的难度！

我把这个程序做成了odc文件，可以解压后在winbugs中调试，希望
懂的朋友献爱心帮朋友解决这个问题，共同进步！

model;
{
   for( j in 1 : J ) {
      for( i in 1 : I ) {
         y[i , j] ~ dnorm(mu[i , j],tau)
      }
   }
   for( j in 1 : J ) {
      for( i in 1 : I ) {
         mu[i , j] <- inprod(x[j , 1:4],beta[1:4 , i])
      }
   }
   for( i in 1 : I ) {
      beta[1:4 , i] ~ dmnorm(mu[1:4 , i],R[1:4 , 1:4])
   }
   tau ~ dgamma(0.001,0.001)
   sigama <- 1 / sqrt(tau)
   for( i in 1 : I ) {
      mu[1:4 , i] <- inprod(gaba[1:4 , 1:3],z[1:3 , i])
   }
   R[1:4 , 1:4] ~ dwish(omega,4)
gaba[1:4 , 1:3]~ dmnorm(mean,pre)
}


vec&pound;¨gaba&pound;&copy;~dmnorm&pound;¨mean&pound;&not;pre&pound;&copy;


list(J=8,I=10,
x=structure(.Data=c(1,1,1,1,1,2,1,2,2,0,0,1,1,0,0,0,0,2,1,1,2,1,1,0,1,0,1,1,1,1,0,0),.Dim = c(8,4)),
z = structure(.Data = c(1,1,1,1,0,0,1,0,0,1,10,36,17,21,70,54,64,40,31,25,600,3000,3000,2500,2000,4000,800,1300,1500,1000),
.Dim = c(3,10)),
mean=structure(.Data=c(0,0,0,
0,0,0,
0,0,0,
0,0,0),Dim=c(4,3))
pre=structure(.Data=c(0.001,0,0,0,0,0,0,0,0,0,0,0,
0,0.001,0,0,0,0,0,0,0,0,0,0,
0,0,0.001,0,0,0,0,0,0,0,0,0,
0,0,0,0.001,0,0,0,0,0,0,0,0,
0,0,0,0,0.001,0,0,0,0,0,0,0,
0,0,0,0,0,0.001,0,0,0,0,0,0,
0,0,0,0,0,0,0.001,0,0,0,0,0,
0,0,0,0,0,0,0,0.001,0,0,0,0,
0,0,0,0,0,0,0,0,0.001,0,0,0,
0,0,0,0,0,0,0,0,0,0.001,0,0,
0,0,0,0,0,0,0,0,0,0,0.001,0,
0,0,0,0,0,0,0,0,0,0,0,0.001),.Dim=c(12,12)),
omega=structure(.Data = c(0.1,0,0,0,0,0.1,0,0,0,0,0.1,0,0,0,0,0.1),.Dim = c(4,4)),
y = srtucture(.Data = c(1,2,3,4,5,6,7,8,
3,5,4,3,6,7,1,2,
2,2,5,2,4,3,2,3,
2,7,4,3,2,7,5,2,
3,5,2,5,6,5,3,2,
1,2,3,4,7,5,4,1,
2,3,4,2,2,3,7,1,
1,4,2,7,5,4,2,6,
2,3,7,5,4,3,2,4,
2,1,3,4,2,4,5,2),
.Dim(10,8)),

Inits
    list(theta = c(0.0, 0.0), X = c(0.0, 0.0, 0.0))
楼主朋友，程序是你编的，你最了解，希望你把初值补上。

这个模型好像winbugs确实做不了，好像只能用R做了，哎。。不知道哪位R高手会用MCMC来做的