我把这个程序做成了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来做的

语法或数据结构错误的地方
请先依据下列修改看看
model{       
for (i in 1:Subjects){
for (j in 1:Replications){
x[i,j,1:3] ~ dmnorm(mu[i,1:3],P[,])}
mu[i,1:3] ~ dmnorm(theta[1:3],T[,])}
for (k in 1:3){
theta[k] ~ dnorm(0,0.0001)
ICC[k] <- V[k,k]/(V[k,k]+W[k,k])}         # Intra-class correlation
P[1:3,1:3] ~ dwish(U[1:3,1:3],3)
T[1:3,1:3] ~ dwish(R[1:3,1:3],3)
W[1:3,1:3] <- inverse(P[1:3,1:3])
V[1:3,1:3] <- inverse(T[1:3,1:3])
p.mu[1:3] ~ dmnorm(theta[1:3],T[,])      # Predicted means
p.x[1:3] ~ dmnorm(p.mu[1:3],P[,])         # Predicted measurements
for (y in 1:2){
        for (z in (y+1):3){
        bias[y,z] <- theta[y]-theta[z]         # Bias estimation
        diff[y,z] <- p.x[y]-p.x[z]}}}               # Limits of agreement

INITIALS
list(theta=c(127,143,143),
P=structure(.Data=c(1,0,0,0,1,0,0,0,1),.Dim=c(3,3)),
T=structure(.Data=c(1,0,0,0,1,0,0,0,1),.Dim=c(3,3)))

DATA
list(Subjects=85, Replications=3,
x=structure(.Data=c(100,98,122,106,98,...(data omitted)...,128),.Dim=c(85,3,3)),
U=structure(.Data=c(0.100,0.005,0.005,
          0.005,0.100,0.005,
          0.005,0.005,0.100), .Dim=c(3,3)),
R=structure(.Data=c(0.100,0.005,0.005,
          0.005,0.100,0.005,
          0.005,0.005,0.1), .Dim=c(3,3)))

非常感谢epoh大师！
Dim=c(85,3,3)),这个85，3，3是什么意思？
Subjects=85, Replications=3,又是什么意思？
附件中程序运行时提示：node dimension doesnt match,如何修改？

附件中有odc文件
model{        
for (i in 1:Subjects){
for (j in 1:Replications){
x[i,j,1:3] ~ dmnorm(mu[i,1:3],P[,])}
mu[i,1:3] ~ dmnorm(theta[1:3],T[,])}
for (k in 1:3){
theta[k] ~ dnorm(0,0.0001)
ICC[k] <- V[k,k]/(V[k,k]+W[k,k])}         # Intra-class correlation
P[1:3,1:3] ~ dwish(U[1:3,1:3],3)
T[1:3,1:3] ~ dwish(R[1:3,1:3],3)
W[1:3,1:3] <- inverse(P[1:3,1:3])
V[1:3,1:3] <- inverse(T[1:3,1:3])
p.mu[1:3] ~ dmnorm(theta[1:3],T[,])      # Predicted means
p.x[1:3] ~ dmnorm(p.mu[1:3],P[,])         # Predicted measurements
for (y in 1:2){
        for (z in (y+1):3){
        bias[y,z] <- theta[y]-theta[z]         # Bias estimation
        diff[y,z] <- p.x[y]-p.x[z]}}}               # Limits of agreement

DATA
list(Subjects=85, Replications=3,
x=structure(.Data=c(100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98, 100,98,122,106,98),.Dim=c(85,3,3)),
U=structure(.Data=c(0.100,0.005,0.005,
          0.005,0.100,0.005,
          0.005,0.005,0.100), .Dim=c(3,3)),
R=structure(.Data=c(0.100,0.005,0.005,
          0.005,0.100,0.005,
          0.005,0.005,0.1), .Dim=c(3,3)))
INITIALS
list(theta=c(127,143,143),
P=structure(.Data=c(1,0,0,0,1,0,0,0,1),.Dim=c(3,3)),
T=structure(.Data=c(1,0,0,0,1,0,0,0,1),.Dim=c(3,3)))

.Dim=c(85,3,3)
请看Bayesian Modeling Using WinBUGS.pdf
CHAPTER 3,page 126/518

3.4.6.2 List data format
  Arrays: using the same syntax as in matrices but the
          .Dim argument will have at least three values
          specifying the length of each corresponding dimension

假设数据结构如下:  
dimension 5 x 3 x 2

  , , 1
        [,1]  [,2]  [,3]
  [1,]    1    2     3
  [2,]    4    5     6
  [3,]    7    8     9
  [4,]   10   11    12
  [5,]   13   14    15

  , , 2
        [,1]  [,2]  [,3]
  [1,]   16    17    18
  [2,]   19    20    21
  [3,]   22    23    24
  [4,]   25    26    27
  [5,]   28    29    30

first define(x11j;j=1,2),then elements(x12j;j=1,2),(x13j;j=1,2),
followed by (x21j;j=1,2),then elements(x22j;j=1,2),(x23j;j=1,2),
........

语法如下:
  A=structurec(.Data=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,
                16,17,18,19,20,21,22,23,24,25,26,27,28,29,30),
               .Dim=c(5,3,2))

########
请注意,相同语法在R/S-PLUS,
结果是不同的 !!   

#########
这篇文献是我两年前看的
当时也是因为关注multivariate hierarchical Bayesian
想到winbugs 可能对楼主有帮助,所以贴上
其实现在也忘光了!!
没想到zhangtao兄对此也有兴趣
这篇文献及data file,winbugs code 可由此下载
  http://www.biomedsearch.com/nih/ ... ch-to/19161599.html

非常感谢epoh大师的指导！

node dimension doesnt match
是因为 .Dim=c(85,3,3)
需要数据 85 x 3 x 3 = 765,你给的不够多

我把数据附上
burn in   :  95,000
iteration : 200,000
chain     :       1

这个结果跟作者文献.

page 7/13,table 1 相近.



node    mean  sd       MC error 2.5%   median 97.5% start sample
ICC[1] 0.9609 0.007197 4.188E-5 0.9452 0.9615 0.9734 95001 105000
ICC[2] 0.9596 0.007439 4.37E-5  0.9434 0.9602 0.9725 95001 105000
ICC[3] 0.9206 0.01423  6.766E-5 0.8899 0.9216 0.9453 95001 105000
node   mean  sd     MC error       2.5% median 97.5% start sample
theta[1]125.2 3.082  0.1709         120.2 124.9 131.3 95001 105000
theta[2]125.2 3.052  0.1692         120.2 124.8 131.2 95001 105000
theta[3]141.1 3.276  0.1449         135.1 141.0 147.6 95001 105000

######################
作者文献

θm     126.9 (120.6, 133.5)   126.9 (120.5, 133.3)   142.6 (135.9, 149.3)
ICCm     0.96 (0.95, 0.97)      0.96 (0.94, 0.97)      0.92 (0.89, 0.95)

multimodel.rar
大小:(1.69 KB)

 马上下载

本附件包括：

• multimodel.bug

T[1:3,1:3] ~ dwish(R[1:3,1:3],3)
epoh老师，dwish是个什么函数？我在winbugs中还没有见过。

Continuous Multivariate
1.Dirichlet
    p[] ~ ddirch(alpha[])

2.Multivariate Normal
    x[] ~ dmnorm(mu[], T[,])

3.Multivariate Student-t
    x[] ~ dmt(mu[], T[,], k)

4.Wishart
    x[,] ~ dwish(R[,], k)

WinBUGS User Manual Version 1.4
  page 59/60
  http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/manual14.pdf

非常感谢epoh大师！

zhangtao兄
mcmc.src  原本就可运行.
garch.src   需要用"EAFE"当index的数据
                   所以我拿garch.asc  (265 x 1)
                   转换成garch.dht
                   请输入
                    load x[265,1]=garch.asc
                   call saved(x,"garch.dht","EAFE")  /*create garch.dht*/;

                    and revised garch.src
                    /*dsn = "indx";*/
                    dsn = "garch.dht";
                    就可运行.
                    运行后会存贮为文件tgarch.out
garch_p.src 需要自己定义参数,才能运行.

garch.asc
     http://faculty.washington.edu/rons/garch.html

epoh老师，您上面的问题的，原始数据应该在这里面有，请您查看

http://faculty.washington.edu/rons/software/garch.asc

zhangtao兄
FGARCHPROC.prg 修改如下即可运行
1. startvalues={0.01,0.01,0.01,0.01};

   /*{b,f,g,vcv,ret,hhat}=FIGARCH11(x,0.2 1 2 3 1 1)*/

   {b,f,g,vcv,ret,hhat}=figarch11(x,startvalues);

2. /* __output=0;*/

   output file = figarch.out reset;

3.你的数据长度是459

   所以proc figarch11_loglike(x,e);

   proc figarch11_hhat(x,e);

   t必须配合修改为459(原作者数据长度1000)

4.修改后FGARCHPROC.prg,结果figarch.out如下:


  

FIGARCH11_proc.rar
大小:(2.19 KB)

 马上下载

非常感谢，深切感受到相互帮助的感动！
也感受到Gauss确实是非常好的语言！

刚看了楼主的方程
应该可以考虑用
R Package ｀bayesm＇
function rhierLinearModel()
   implements a Gibbs Sampler for hierarchical linear models
   with a normal prior.

# Model:
# nreg regression equations
#        y_i = X_ibeta_i + epsilon_i  
#        epsilon_i ~ N(0,tau_i)
#             nvar X vars in each equation
#
# Priors:
#        tau_i ~ nu.e*ssq_i/chisq(nu.e)  tau_i is the variance of epsilon_i
#        beta_i ~ N(ZDelta[i,],V_beta)
#               Note:  ZDelta is the matrix Z * Delta; [i,] refers to ith row of this product!
#
#          vec(Delta) | V_beta ~ N(vec(Deltabar),Vbeta (x) A^-1)
#          V_beta ~ IW(nu,V)  or V_beta^-1 ~ W(nu,V^-1)
#              Delta, Deltabar are nz x nvar
#              A is nz x nz
#              Vbeta is nvar x nvar
#        
#          NOTE: if you don't have any z vars, set Z=iota (nreg x 1)


详细数学式子请看
   http://hosho.ees.hokudai.ac.jp/~ ... ierLinearModel.html

23# epoh
epoh老师，您的建议很好！但是在您提供的链接上没有
bayesm包？

Package :‘bayesm’
Title :Bayesian Inference for Marketing/Micro-econometrics
Author: Peter Rossi.
下载bayesm_2.2-4.zip 由本机安装
  http://cran.r-project.org/web/packages/bayesm/index.html

要效果加倍,当然是参考Peter Rossi.的书籍
Bayesian Statistics and Marketing
  https://bbs.pinggu.org/thread-824880-1-1.html

25# epoh
epoh老师，我这几年对几个概念理解不好，如果您对这方面了解的话，想请教您讲一讲：
1、金融和时间序列计量中的“谱分析”和“谱密度”函数的思想是什么？有什么用？有什么优势？
2、什么是“鞅差分序列”？
3、griddy-gibbs抽样和传统的gibbs抽样有什么区别和联系？
非常感谢！

依题意及数据
应该是nreg=4; nobs=10; nvar=8
不过楼主给的x数据不足
所以我修改数据并改为
nreg=3; nobs=10; nvar=3

#############
library(bayesm)
R=10
nreg=3; nobs=10; nvar=3
Z=cbind(c(rep(1,nreg)),3*runif(nreg)); Z[,2]=Z[,2]-mean(Z[,2])
x = structure(.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,2,7,4,3,2,7,2,3,4,2),.Dim=c(10,3,3))

y=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),.Dim = c(10,3))

iota=c(rep(1,nobs))
regdata=NULL
for (reg in 1:3) {
X=cbind(iota,x[,,reg])
Y=y[,reg]
regdata[[reg]]=list(y=Y,X=X) }
Data1=list(regdata=regdata,Z=Z)
Mcmc1=list(R=R,keep=1)
out=rhierLinearModel(Data=Data1,Mcmc=Mcmc1)

##########
out
Starting Gibbs Sampler for Linear Hierarchical Model
    3  Regressions
    2  Variables in Z (if 1, then only intercept)

Prior Parms:
Deltabar
     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
A
     [,1] [,2]
[1,] 0.01 0.00
[2,] 0.00 0.01
nu.e (d.f. parm for regression error variances)=  3
Vbeta ~ IW(nu,V)
nu =  7
V
     [,1] [,2] [,3] [,4]
[1,]    7    0    0    0
[2,]    0    7    0    0
[3,]    0    0    7    0
[4,]    0    0    0    7

MCMC parms:
R=  10  keep=  1

MCMC Iteration (est time to end - min)
  Total Time Elapsed:  0

> names(out)
[1] "Vbetadraw" "Deltadraw" "betadraw"  "taudraw"

zhangtao 兄
griddy-gibbs抽样和传统的gibbs抽样有什么区别
这个答案可在Analysis of Financial Time Series_2ed找到

Analysis of Financial Time Series_2ed
chap12 Markov Chain Monte Carlo Methods with Applications
page 577/638
12.4.3 Griddy Gibbs
1. In financial applications, model may contain some nonlinear parameters
      (e.g., the moving-average parameters in an ARMA model or
       the GARCH parameters in a volatility model)
2. nonlinear parameters do not have a closed-form expression,
     implementing a Gibbs sampler in this situation may become complicated
3. Tanner (1996) describes a simple procedure to obtain random
    draws in a Gibbs sampling when the conditional posterior distribution
    is univariate.
4. The method is called the Griddy Gibbs sampler and is widely applicable
5 .However, the method could be inefficient in a real application

详细的 Griddy Gibbs Algorithm
请参考page 553/638
范例请参考 Example 12.3.GARCH(1,1)
ht is drawn by using the Griddy Gibbs with 500 grid points ,
and the range of ht is (0, 1.5s^2),where s^2 is the
sample variance of the log return rt
可惜没有S-PLIU SOURCE CODE

我曾经学习北大何书元教授的《应用时间序列分析》和
华东师大茆师松教授的〈高等数理统计〉等书，就是实在弄不明白
“谱分析”、“谱密度函数”、谱密度核函数、谱密度极大似然估计、最大熵谱估计、还有
Fisher信息阵是怎么来的？有什么用？思想是什么？
在baidu and google上也做过搜索，在人大论坛上也发过帖子，
但是也没人知道，但是这些东西不清楚。

Fisher information matrix
可以利用matlab ecmnfish.m算出

%%%%%%%%%%
ecmnfish
  Fisher information matrix
  Fisher = ecmnfish(Data, Covariance, InvCovariance, MatrixFormat)

Arguments:

  Data: NUMSAMPLES-by-NUMSERIES matrix of observed multivariate normal data
  Covariance :NUMSERIES-by-NUMSERIES matrix with covariance estimate of Data
  InvCovariance:(Optional) Inverse of covariance matrix: inv(Covariance)

%%%%%%%%
matlab Command Window
>> ecmguidemo

这是matlab 自动产生

3 series
60 observations
0% missing data

比较用Hessian 及 Fisher information matrix的结果

   Mean(Hessian)      Std.Dev(Hessian)        Mean(Fisher)       Std.Dev(Fisher)
0.023355856212147877 0.005975609497120777 0.02335585621214785  0.005975609497120777
0.03739210439061549  0.013830706723727054 0.033608644398810866 0.012373501469109024
0.053544968623126306 0.031212803971525954 0.05267961862804869  0.03040015987039729


%%%%%%%%
不知zhangtao兄所谓的谱分析
谈的是否底下的内容
我的理解是假设你的数据属cycle的话
通常会使用Spectral Analysis
这在Time series analysis and its applications with R examples _2ed
chap 4 Spectral Analysis and Filtering
都有详细讨论到
4.3 The Spectral Density
4.4 Periodogramand Discrete Fourier Transform
4.5 Nonparametric Spectral Estimation
4.8 Parametric Spectral Estimation
请先行参考.

############
martingale difference sequence:

个人认为看完底下这份教材
page 1,2 的定义及举例garch(1,1)
应该不难理解.
   http://www.ssc.wisc.edu/~bhansen/718/NonParametrics13.pdf


Simple time series models are built using fundamental "“shocks”et";
typically normalized so that Eet = 0:
(1) The most basic shock et is iid.
(2) Martingale Difference Sequence (MDS). et is a MDS if E (et |et-1,et-2,..) = 0
    it is a serially uncorrelated (but not necessarily an independent) process
(3) White noise. If E(et)(et-h) = 0 for all h>=1
   
     iid > MDS > White noise

An iid shock is a MDS,a MDS is white noise, but the reverse is not true.
An iid shock is unpredictable. A MDS is unpredictable in the mean.