WinBUGS在统计分析中的应用(第一部分)

开篇词

首先非常感谢COS论坛提供了这样一个良好的平台,敝人心存感激之余,也打算把一些学习心得拿出来供大家分享,文中纰漏之处还请各位老师指正。下面我将以WinBUGS的统计应用为题,分几次来谈一谈WinBUGS这个软件。其中会涉及到空间数据的分析、GeoBUGS的使用、面向R及SPLUS的接口包R2WinBUGS的使用、GIS与统计分析等等衍生出的话题。如有问题,请大家留下评论,我会调整内容,择机给予回答。

第一节 什么是WinBUGS?

WinBUGS对于研究Bayesian统计分析的人来说,应该不会陌生。至少对于MCMC方法是不陌生的。WinBUGS (Bayesian inference Using Gibbs Sampling)就是一款通过MCMC方法来分析复杂统计模型的软件。其基本原理就是通过Gibbs sampling和Metropolis算法,从完全条件概率分布中抽样,从而生成马尔科夫链,通过迭代,最终估计出模型参数。引入Gibbs抽样与MCMC的好处是不言而喻的,就是想避免计算一个具有高维积分形式的完全联合后验概率公布,而代之以计算每个估计参数的单变量条件概率分布。具体的算法思想,在讲到具体问题的时候再加以叙述,在此不过多论述。就不拿公式出来吓人了(毕竟打公式也挺费劲啊)。

第二节 为什么要用WinBUGS?

第一、因为同类分析软件中它做得最好。同类的软件:OpenBUGS、JAGS等在成熟度、灵活性以及兼容性方面和它相比还有一定距离。在处理spatial data的方面,它采用了Gibbs抽样和MCMC的方法,在模型支持方面又具有极大的灵活性,较之名声大噪的GeoR包,虽然也实现了bayesian的手法,但是灵活性还是不及WinBUGS。

第二、因为它免费。免费的东西总有吸引人之处。

第三、有各色的R包为WinBUGS实现了针对R的、SPLUS的、Matlab的软件接口。只要你喜欢,就直接在R下调用WinBUGS吧,无非是装个R2WinBUGS包,简简单单。

第四、详细的文档、帮助、指导、范例。当然没有中文版的,小小一点遗憾。

第三节 如何得到WinBUGS?

WinBUGS目前是一款免费的软件,去http://www.mrc-bsu.cam.ac.uk/bugs/下载就好了。不过要用高级功能(如GeoBUGS)的话,还是去http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/contents.shtml注册一下好了,挺方便的。系统会立即把注册码发到你邮箱(真是好人啊)。不过只可以用一个月。这倒无妨,到时在注册一下就好了。

第四节 初试WinBUGS

WinBUGS-GUI
WinBUGS-GUI

我们先找一个例子来实际地运行一下WinBUGS,感受一下基本的操作流程,然后再考虑高级的操作。

第一步,打开WinBUGS,通过菜单File -> New新建一个空白的窗口

第二步,在第一步中新建的空白窗口中输入三部分内容:模型定义、数据定义、初始值定义(代码见附录)

第三步,点击菜单Model -> Specification,弹出一个Specification Tool面板。

第四步,在第二步中的提到的那个窗口中,将model这个关键字高亮起来,点击check model。你会看到WinBUGS的左下角状态栏上显示”model is syntactically correct.”

第五步,把定义的data前的关键字list也高亮起来,点Specification Tool面板上的load data

第六步,改Specification Tool面板上的马尔科夫链的数目,默认为1就好了

第七步,点击Specification Tool面板上的compile

第八步,把定义的初始值中的list关键字也高亮起来,再点击Specification Tool面板上的load inits

第九步,关了Specification Tool面板

第十步,点击菜单Inference -> Samples,弹出一个Sample Monitor Tool面板。

第十一步,在Sample Monitor Tool面板的node中填要估计的参数名,这里可以是tau, alpha0, alpha1, b, 把它们一个一个填在node中,逐一点set。

第十二步,关了Sample Monitor Tool面板

第十三步,点击菜单Model -> Update,弹出一个Update Tool面板。

第十四步,将Update Tool面板中的updates改大点,比如50000,点update按钮。

第十五步,运行完后,关了Update Tool面板

第十六步,点击菜单Inference -> Samples

第十七步,在弹出的Sample Monitor Tool面板上选一个node

第十八步,点history看所有迭代的时间序列图,点trace看最后一次迭代的时间序列图,点auto cor看correlogram时间序列图,点stat看参数估计的结果

Estimation results by WinBUGS 1.4
Estimation results by WinBUGS 1.4

附第二步中的代码如下:

#MODEL
model
{
    for (i in 1:N) {
        O[i] ~ dpois(mu[i])
        log(mu[i]) <- log(E[i]) + alpha0 + alpha1 * X[i]/10 +
            b[i]
        # Area-specific relative risk (for maps)
        RR[i] <- exp(alpha0 + alpha1 * X[i]/10 + b[i])
    }
    # CAR prior distribution for random effects:
    b[1:N] ~ car.normal(adj[], weights[], num[], tau)
    for (k in 1:sumNumNeigh) {
        weights[k] <- 1
    }
    # Other priors:
    alpha0 ~ dflat()
    alpha1 ~ dnorm(0, 1e-05)
    tau ~ dgamma(0.5, 5e-04)
    # prior on precision
    sigma <- sqrt(1/tau)
    # standard deviation
}
#DATA

list(N = 56, O = c(9, 39, 11, 9, 15, 8, 26, 7, 6,
    20, 13, 5, 3, 8, 17, 9, 2, 7, 9, 7, 16, 31, 11, 7, 19, 15,
    7, 10, 16, 11, 5, 3, 7, 8, 11, 9, 11, 8, 6, 4, 10, 8, 2,
    6, 19, 3, 2, 3, 28, 6, 1, 1, 1, 1, 0, 0), E = c(1.4, 8.7,
    3, 2.5, 4.3, 2.4, 8.1, 2.3, 2, 6.6, 4.4, 1.8, 1.1, 3.3, 7.8,
    4.6, 1.1, 4.2, 5.5, 4.4, 10.5, 22.7, 8.8, 5.6, 15.5, 12.5,
    6, 9, 14.4, 10.2, 4.8, 2.9, 7, 8.5, 12.3, 10.1, 12.7, 9.4,
    7.2, 5.3, 18.8, 15.8, 4.3, 14.6, 50.7, 8.2, 5.6, 9.3, 88.7,
    19.6, 3.4, 3.6, 5.7, 7, 4.2, 1.8), X = c(16, 16, 10, 24,
    10, 24, 10, 7, 7, 16, 7, 16, 10, 24, 7, 16, 10, 7, 7, 10,
    7, 16, 10, 7, 1, 1, 7, 7, 10, 10, 7, 24, 10, 7, 7, 0, 10,
    1, 16, 0, 1, 16, 16, 0, 1, 7, 1, 1, 0, 1, 1, 0, 1, 1, 16,
    10), num = c(3, 2, 1, 3, 3, 0, 5, 0, 5, 4, 0, 2, 3, 3, 2,
    6, 6, 6, 5, 3, 3, 2, 4, 8, 3, 3, 4, 4, 11, 6, 7, 3, 4, 9,
    4, 2, 4, 6, 3, 4, 5, 5, 4, 5, 4, 6, 6, 4, 9, 2, 4, 4, 4,
    5, 6, 5), adj = c(19, 9, 5, 10, 7, 12, 28, 20, 18, 19, 12,
    1, 17, 16, 13, 10, 2, 29, 23, 19, 17, 1, 22, 16, 7, 2, 5,
    3, 19, 17, 7, 35, 32, 31, 29, 25, 29, 22, 21, 17, 10, 7,
    29, 19, 16, 13, 9, 7, 56, 55, 33, 28, 20, 4, 17, 13, 9, 5,
    1, 56, 18, 4, 50, 29, 16, 16, 10, 39, 34, 29, 9, 56, 55,
    48, 47, 44, 31, 30, 27, 29, 26, 15, 43, 29, 25, 56, 32, 31,
    24, 45, 33, 18, 4, 50, 43, 34, 26, 25, 23, 21, 17, 16, 15,
    9, 55, 45, 44, 42, 38, 24, 47, 46, 35, 32, 27, 24, 14, 31,
    27, 14, 55, 45, 28, 18, 54, 52, 51, 43, 42, 40, 39, 29, 23,
    46, 37, 31, 14, 41, 37, 46, 41, 36, 35, 54, 51, 49, 44, 42,
    30, 40, 34, 23, 52, 49, 39, 34, 53, 49, 46, 37, 36, 51, 43,
    38, 34, 30, 42, 34, 29, 26, 49, 48, 38, 30, 24, 55, 33, 30,
    28, 53, 47, 41, 37, 35, 31, 53, 49, 48, 46, 31, 24, 49, 47,
    44, 24, 54, 53, 52, 48, 47, 44, 41, 40, 38, 29, 21, 54, 42,
    38, 34, 54, 49, 40, 34, 49, 47, 46, 41, 52, 51, 49, 38, 34,
    56, 45, 33, 30, 24, 18, 55, 27, 24, 20, 18), sumNumNeigh = 234)
#INITIAL VALUES
list(tau = 1, alpha0 = 0, alpha1 = 0, b = c(0, 0,
    0, 0, 0, NA, 0, NA, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

WinBUGS在统计分析中的应用 第一部分完

WinBUGS在统计分析中的应用(第一部分)》有80个想法

  1. 我会在第二部分中将相关的理论部分加上,看来得要学习一下怎么嵌入LaTex了.

    1. 请问如何输入幂的形式的,比如a^b,但是中间那个符号不能识别

      1. ”请问如何输入幂的形式的,比如a^b,但是中间那个符号不能识别“ 请问这个问题您解决了没,求教

  2. 嗯,楼主辛苦了!跟着一点一点学,期待着空间数据分析。。。

  3. 贴一首诗吧,刚BUGS team发过来的,非常有意思:

    Each year you wait with bated breath,
    The old WinBUGS key nearing death,
    And will the brand new key appear
    In time to join the festive cheer?
    The waiting’s over – raise your glass,
    And drink to rituals that pass.
    Relax, sit back and have a chortle;
    This time your WinBUGS key’s immortal.

    1. 可以。Google Openbugs 的主页,上面又说明。我更喜欢linux下的版本。
      另外,R里面有 R2bugs 的package,在数据不大的情况下,很方便。

  4. 前几天上课时听说了这个软件,真是及时雨!谢谢了!

  5. 小弟现正在剑桥mrc-bsu做postdoc,具体的项目就是BUGS的开发以及在生物及医学方面应用。WinBUGS这个软件是我两个老板David Spiegehalter,Dave Lunn 和其它一些牛人共同开发的。我们现在正在从WinBUGS 转向openBUGS,目的是将它做成open source的软件以应用在更广的领域。 我现在正在开发BUGS中的WBDiff部分并将它应用在二型糖尿病的动态系统的数据分析中。有兴趣或有问题的同学可以和我联系:

    [email protected]

    还有我们这里每年会举办3-4次BUGS的培训,2天的课程,在英国或能到英国出差的同学有兴趣的话可以参加,主讲人是David speigehalter 和 Dave Lunn。

    1. 请问在
      winbugs中我编译FRAILTY模型时提示
      educational version cannot do this model 难道这个软件还有别的版本???我用的版本是1.4
      我编译的模型是帮助文件Examples Volume 1中的最后一个文件Cox regression with random effects

  6. 希望能和您联系,我这边看了gibbs sampling有不少问题,不知道您可否提供帮助?

  7. 弱弱的问一句,执行到第四步check model的时候,winbugs坐下角显示的不是“model is syntactically correct”,而是在alpha1 ~ dnorm(0, 1e-05)一句1e处显示e处应为”expected right parenthesis,在代码最后0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))指示”invalid or unexpected token scanned”,
    这要怎么弄呢? 第一次接触这个软件,还不会使,请高手指点。

    1. 数字格式不对。参见下面的回复。

      To齐韬:你可能需要修改正文中的代码,我在OpenBUGS 3.0.7中运行你的代码会出错,错误信息和这里一样。我检查了一下,是数字格式的问题。

      1. 请问就是把您给的代码复制进WINBUGS,到compile时,出现educational version cannot do this model.什么原因阿?对这个软件不甚掌握,论文却急用请指教,谢谢

  8. 把 alpha1 ~ dnorm(0, 1e-05)
    tau ~ dgamma(0.5, 5e-04)
    替换成

    alpha1 ~ dnorm(0, 0.00001)
    tau ~ dgamma(0.5, 0.0005)

    这样就可以了

    1. 出错的原因是数字的格式不对,应该写作1.0E-05(1是整数,1.0才是浮点数;e要大写)

      这里设置的正态分布实际上应该是想要一个flat prior,所以用dflat()就可以了。

      1. 您好,请问进行MH算法需要另外编辑程序吗?

      2. 不需要,要是MH都要编程,那WinBUGS就没有存在必要了……呃,确切地说,我说的是OpenBUGS……我从没用过WinBUGS,因为它已经确定没有继续的技术支持了,它的开发已经停止。OpenBUGS还在更新中。

  9. 请问winbugs能计算递归吗?想用之解决一个随机变量自我n重卷积的问题。用递归来算,比如一个正态分布s~(mu,tau),自我卷积n次求结果。该怎么弄呢?谢谢大虾们指教。

  10. 请教一个问题:我运行第五步,load data的时候,软件的左下角显示expected variable name, 然后点击compile的时候,软件界面上没有任何变化,load inits和gen Init按钮是灰色的,不能点击运行。是怎么回事啊,请各位朋友帮帮忙,谢谢了。

    1. 你要选中所有的数据,然后点load data,或者至少你要选中list这个词,这样BUGS才知道数据的起止。如果选择范围不恰当,BUGS就不能正确载入数据。

      另一个可能性是你的模型中使用了数据没有提供的变量名,这种情况下,你要自己检查模型中的变量,如果有遗漏,就在数据的list()中加上。

  11. model
    {
    for (i in 1 : I) {
    for(j in 1 : J) {
    for(k in 1:K) {
    epsilon[i,j,k] ~ dnorm(0.0,tau[i])
    y[i,j,k] <-alpha[i,k] + beta[i,k] * theta[j,k] + epsilon[i,j,k]
    }
    }
    }

    list(
    I = 5,J = 5,K = 4,
    y = structure(
    .Data = c(
    6.8,7.2,7.35,7.63,

    .Dim = c(5,5,4)
    )
    )

    编译出现multiple definitions of node y[1,1,1]
    如果去掉数据部分的y定义就可以编译过去。
    希望大家帮我解决下这个问题,万分感谢!

    1. 去掉y即使模型语法正确,程序肯定也无法运行。

      你再仔细看看你这里y的定义,c的回括号打到哪里去了?

      还有,一共4个数字,请问该如何定义为一个5x5x4的数组?

    2. 请问,上面这位仁兄,你的问题解决了吗?我现在也遇到与你同样的问题。到底问题出在哪里呢?

  12. model                                                               
    {  
    #distribution of Ys                                                             
    ###################                                                             
    for (i in 1:N) {
        ysigmadet[i]<-exp(th[i,1]+th[i,2])*(1-rhoep[i]*rhoep[i]);        
        Yisigma2[i,1,1] <- exp(th[i,2])/ysigmadet[i];
        Yisigma2[i,2,2] <- exp(th[i,1])/ysigmadet[i];
        Yisigma2[i,1,2] <- -rhoep[i]*exp(0.5*th[i,1]+0.5*th[i,2])/ysigmadet[i];             Yisigma2[i,2,1] <- Yisigma2[i,1,2];
        Y[i,1:2]~ dmnorm(muy[],Yisigma2[i,,]);  
        }
    muy[1]<-0;
    muy[2]<-0;
    thmean[1,1] <- mu1;   
    thmean[1,2] <- mu2;
    th[1,1]~dnorm(thmean[1,1],itaua2);
    th[1,2]~dnorm(thmean[1,2],itaub2);
    sig1[1]<-exp(0.5*th[1,1]);
    sig2[1]<-exp(0.5*th[1,2]);
    q[1]~dnorm(psi0,itau2);
    rhoep[1]<-(exp(q[1])-1)/(exp(q[1])+1);
    for (i in 2:N) {
    		thmean[i,1] <- mu1 + phi1*(th[i-1,1]-mu1);   
    		thmean[i,2] <- mu2 + phi2*(th[i-1,2]-mu2);
    		th[i,1]~dnorm(thmean[i,1],itaua2);
    		th[i,2]~dnorm(thmean[i,2],itaub2);
    		sig1[i]<-exp(0.5*th[i,1]);
    		sig2[i]<-exp(0.5*th[i,2]);
    		qmean[i]<-psi0+psi*(q[i-1]-psi0);
    		q[i]~dnorm(qmean[i],itau2);
    		rhoep[i]<-(exp(q[i])-1)/(exp(q[i])+1);
       }                                                                            
    #distribution of phi, mu, rhoep                                                    
    ###########################                                                     
    phi1star ~ dbeta(20,1.5);   
    phi1 <- 2*phi1star -1;  
    phi2star ~ dbeta(20,1.5);   
    phi2 <- 2*phi2star -1;
    psistar ~ dbeta(20,1.5);   
    psi <- 2*psistar -1;
    
    itaua2 ~ dgamma(2.5,0.025);                                                      
    taua <- sqrt(1/itaua2);
    itaub2 ~ dgamma(2.5,0.025);                                                      
    taub <- sqrt(1/itaub2);                                             
    itau2 ~ dgamma(2.5,0.025);                                                      
    tau <- sqrt(1/itau2);                                             
    mu1 ~ dnorm(0,0.04);
    mu2 ~ dnorm(0,0.04);
    psi0~dnorm(0.7,0.1);
    } 
    list(phi1star=0.99, phi2star=0.99,mu1=0,mu2=0,itaua2=100,itaub2=100,psistar=0.99,psi0=1.9,itau2=100)
  13. 这是一个DC-MSV模型,大家能帮我看一下,到底是哪个初始值没有定义吗?
    在编译完成之后,提示说有未定义的初始值,是怎么回事?
    急急急,先谢过各位了!!!!!

    1. 这是个体力活儿,自个儿把模型中所有的参数都列出来,然后对比一下是不是所有的参数都有初始值吧

      1. 请教谢老师,我用winbugs做SVCJ(带跳的随机波动率)模型,MODEL检验和上传数据都行,就是到了load inits那里总是提示有uninitialized variables,可我该设的参数的初始值差不多都设了,仔细检查了不知道还有什么变量没设。如果采用软件自动赋值初始化(gen inits)的话,只要数据稍多一点,比如500条记录,软件就半天没反应了,而如果只是几十条数据的话还行。不知道是不是编的程序的问题?另还有什么参数没设初值呢?
        model;
        {
        mu ~ dnorm(1,0.4)
        for( i in 2 : n ) {
        v[i] ~ dnorm(vmean[i],ivd[i])I(0,50)
        }
        kt ~ dnorm( 0.0, 1.0)
        k ~ dnorm( 0.0, 1.0)
        # theta <- kt / k
        for( i in 2 : n ) {
        ksy[i] ~ dnorm(ksymean[i],tauy)
        }
        for( i in 2 : n ) {
        ksv[i] ~ dexp(muv)
        }
        rhoj ~ dnorm( 0.0,4)
        muv ~ dgamma(1,0.1)
        for( i in 2 : n ) {
        j[i] ~ dbern(lambda)
        }
        lambda ~ dbeta(1,20)
        #ev ~ dnorm( 0.0, 1.0)
        #ey ~ dnorm( 0.0, 1.0)
        #e ~ dnorm( 0.0,1.0E-6)
        sigy2 <- 1 / tauy
        tauy ~ dgamma( 1,0.1)
        muy ~ dnorm( 0.0,0.1)
        for( i in 2 : n ) {
        #rmean[i] <- mu + ksymean[i] * j[i]
        rmean[i] <- mu + ksy[i] * j[i]
        }
        tauv ~ dgamma( 1, 0.1)
        sigv2 <- 1 / tauv
        for( i in 2 : n ) {
        #vmean[i] <- v[i – 1] +( kt-k* v[i – 1]) + 1/muv * j[i] + rho * (y[i] – y[i – 1] – mu – ksymean[i] * j[i])
        vmean[i] <- v[i – 1] +( kt-k* v[i – 1]) +ksv[i] * j[i] +sqrt(sigv2)* rho * (y[i] – y[i – 1] – mu – ksy[i] * j[i])
        }
        for( i in 2 : n ) {
        #ivd[i] <- 1 / (rho*rho*sigy2*j[i]+(1-rho*rho)*sigv2*j[i]+v[i – 1]+1/pow(muv,2)*j[i])
        ivd[i] <- tauv/((1-rho*rho)*v[i-1])
        }
        for( i in 2 : n ) {
        #ird[i] <- 1 / (sigy2*j[i]+v[i – 1])
        ird[i] <- 1/v[i-1]
        }
        for( i in 2 : n ) {
        r[i] <- y[i]-y[i-1]
        r[i] ~ dnorm(rmean[i],ird[i])
        }
        rho ~ dunif(-1,1)
        for( i in 2 : n ) {
        ksymean[i] <- muy + rhoj * ksv[i]
        }
        v[1] <- 1
        ksv[1] <- 0.0
        ksy[1] <- 0.0
        j[1]<-1
        }
        list(y=c(3.274324143,3.287163251,3.284207325,3.273737212,3.279863062),
        n=5)
        list(mu=1,k=1,kt=1,tauv=1,rho=1,lambda=1,muy=1,tauy=1,muv=1,rhoj=1)

    2. 您好,我也在做这个dcc-msv模型。代码部分也是这个。可是我到了compile这一步以后就一直提示“vairable N is not defined”.请问你遇到了这个问题吗?谢谢

      1. 请问您的问题解决了吗?我也遇到这个问题

  14. 每次启动winBUGS1.4.3都要弹出Licence Agreement,这是不是没安装好呢?

  15. 第七步时,它说edicational version cannot do this model,是不是必须注册了才行呀

  16. 学《高级统计学》的时候,在介绍MCMC及EM的时候,老师介绍过这个软件,
    非常好!就是我自己编不出程序来。我是统计的爱好老师与初学者。
    另外,问一下谢老师,openbugs在那儿有下载的呢?
    非常感谢!

  17. 请问谢老师,我的问题也是
    它说edicational version cannot do this model,是什么原因呢?谢谢了

      1. 你好,我是WINBUGS的初学者,现在想用来做SV-T模型的MCMC估计,但总提示数据载入不了,还有初始值设置也不明白,请帮忙指导一下,真的很急,谢谢!!!

        model{
        for(i in 1:n)
        {y[i]~dt(0,p[i],omega)
        p[i]<-exp(-theta[i])
        }
        theta[1]~dnorm(mu,itau2)
        for(j in 2:n)
        {theta[j]~dnorm(theta2[j],itau2)
        theta2[j]<-mu+phi*(theta[j-1]-mu)}
        phi<-2*phi1-1
        tau<-sqrt(1/itau2)
        mu~dnorm(0,0.01)
        itau2~dgamma(2.5,0.025)
        phi1~dbeta(20,1.5)
        omega~dchisqr(8)
        }
        
        list(N=236, y=c(-2.863137825,0.236862576,0.708743605,-0.410171084,-0.197900711,-0.017711569,-1.099961877,
        -0.131942478,-0.628070303,0.329370517,-0.425979456,0.003865601,-1.562805432,-1.029900337,0.30855106,
        -1.76405028,0.536937996,0.88903611,-0.478923481,-2.231280334,0.95179952,-0.240613751,-0.444776657,
        0.238243297,2.343495649,-1.345814335,0.409521478,0.938909523,-0.541400821,-1.120087252,-0.465308369,
        0.22970175,-0.442665133,0.208833427,-0.866786092,0.294963792,1.018515417,-0.296166322,-0.055762223,
        -0.336635628,-0.431943542,1.38828862,-0.107762343,-0.502041508,0.094863359,-0.266547417,-0.188964816,
        -2.292527804,-0.355015921,-0.569077156,0.017514077,-0.25962396,0.583479024,0.068167933,0.014009904,
        1.352042569,0.177236882,-0.564620898,0.349787556,-0.707430539,-0.129763575,1.671152083,0.776163597,
        -0.021446623,0.632287045,-0.008614125,0.209323915,-0.289051354,1.681336675,-0.155862607,-0.077731799,
        0.82786481,-1.041519493,0.261535927,-0.412195369,1.112557672,-0.100371081,-1.144249282,0.040589822,
        -0.189283934,0.434968595,1.094428192,0.290705338,-0.0138222,0.283080016,-0.782503525,-0.085722785,
        0.129237661,-0.863186702,0.146550995,0.106793381,0.866076389,-0.290592628,-0.298612901,0.585978478,
        0.060169353,0.871721164,0.067681259,-0.016429211,-0.691862319,-0.099393894,0.397879547,-0.0049166,
        -0.686676288,-0.782562372,-0.017342403,0.106779034,-0.020327393,0.4941572,-0.184977933,0.045901852,
        0.883238737,1.461385829,1.173875952,0.238537867,0.716623196,0.32768553,1.063529252,0.24979205,
        1.218106841,0.438037488,-0.11434406,0.212401723,1.731100734,-0.575671103,-0.861522307,0.323921825,
        -1.171227643,1.687658722,-0.670770608,-0.347428641,0.656840307,0.654451244,0.023281919,-0.413060303,
        -0.236891092,-0.290287567,-2.747778165,0.143692637,-1.71000562,-0.935081269,0.692751369,
        -0.411491889,0.928809898,-0.926643719,1.04862006,0.412139218,-0.193032523,-0.009653293,-0.79853627,
        -0.100936572,0.309291433,0.102463117,-0.050883159,0.694099117,-0.842781601,
        -0.269146402,0.441508007,1.405110928,0.033024447,-0.328643359,-0.323136664,-0.09309922,
        -0.158225088,1.260303693,-0.789521669,-0.337560523,-0.291325249,-1.066279143,
        -0.618890031,0.248983641,0.034759251,0.882100279,0.48269568,-0.066941817,
        -0.184608038,-0.023517406,-0.672468729,0.253540536,0.191561488,-0.087120801,-0.686079288,
        -1.793905963,0.036004183,1.022793487,-1.434730455,0.589690201,-0.217817872,
        -0.145532668,0.401124556,0.627431273,0.224008103,0.485978832,0.082761656,
        -0.662787178,0.939936177,0.179767043,1.528823088,-0.165954634,0.393731707,-0.033624367,
        -0.370122603,0.936833337,-1.404577085,-0.053542798,0.249369071,0.293947579,0.548939249,-0.006808959,
        -0.11640434,-0.288842871,0.870856137,0.706003843,0.278474099,-0.229181126,-0.72134986,
        -0.524715596,0.113296712,-0.787005295,0.586313224,
        -0.47130353,0.155244387,0.154718569,0.025373314,0.50625468,-0.025369976,0.419284333,-0.019841149,
        -0.297545391,-0.131656389,-0.25774229,0.435472774,0.420723963,0.210803862,0.551788416,-0.27312702))
        
        list(mu=1, tau=2)
      2. 首先我想知道你了解这个模型吗?软件不是用来往里面乱塞东西的,尽管BUGS软件做得很简单(通常给定先验、给定数据,剩下的什么都不用管了),但你也得稍微明白一下这些代码的来龙去脉吧。比如你的模型中有个变量n,在数据中有吗?初始值中提供了一个变量tau,可是模型中它是参数吗?

  18. 你好,我刚学WINBUGS软件,现在想用来做SV-N模型,但是在模型的编辑(compile)这一步时,软件总是显示array index is greater than array bound for y 不知为什么,麻烦请帮忙指导一下,真的很急,谢谢!!!

    model
    {
    for(i in 1:N)
    {y[i]~dnorm(0,p[i])
    p[i]<-exp(-theta[i])
    }
    theta[1]~dnorm(mu,itau2)
    for(j in 2:N)
    {theta[j]~dnorm(theta2[j],itau2)
    theta2[j]<-mu+phi*(theta[j-1]-mu)}
    phi<-2*phi1-1
    mu~dnorm(0,0.01)
    itau2~dgamma(2.5,0.025)
    phi1~dbeta(20,1.5)
    }
    list(N=1000,
    y=c(-0.01332,0.215247,0.165196,0.268095,-0.08932,0.784046,0.194879,-0.55396,-0.06357,-0.75685,
    -0.80325,-0.24164,-0.14004,-0.06396,-0.03859,0.03763,-0.41906,-0.2917,-0.01348,-0.01342,
    -0.21579,0.189124,0.215324,0.778916,-0.60096,-0.97844,0.646815,-0.92657,0.341113,-0.08932,
    -0.06394,-0.19079,-0.16536,0.393219,0.242184,0.064052,-0.08932,0.03867,-0.01245,-0.06368,
    -0.06368,-0.01235,0.270668,0.764389,0.249,-0.42764,-0.01135,-0.45268,0.09219,0.300745,
    -0.14142,-0.21945,0.275452,0.434102,-0.03683,-0.0368,0.015803,0.200342,0.016219,0.759039,
    0.123898,-0.14267,-0.19551,-0.19539,-0.56526,-0.45794,0.252919,-0.74644,0.357049,-0.32588,
    -0.0368,0.173698,0.306509,-0.27424,0.14849,0.361428,0.176777,0.204207,0.823436,0.072621,
    -0.11633,-1.51045,-0.77916,-0.14219,0.414057,-0.61912,0.148742,0.255552,0.711109,-0.83659,
    0.150263,-1.33421,-0.71908,0.093948,-0.97636,0.014628,-0.08932,-0.7889,-0.45015,-0.29492,
    0.06484,0.348767,0.454509,0.248824,-1.28074,0.013718,0.323902,-0.03755,0.195911,0.248824,
    -0.29755,0.197099,-0.08932,-0.11539,0.067213,0.38176,-0.06309,0.015692,0.252739,-0.03659,
    -0.19475,0.042482,0.280657,0.149247,0.28292,-0.14258,0.123898,-0.06264,-0.03593,0.124525,
    -0.14283,-0.16953,-0.08932,0.151481,0.286416,-0.08932,0.126021,-0.35843,-0.2236,0.018089,
    -0.06245,-0.16991,-0.08932,-0.08932,0.233432,-0.11626,0.315482,0.154349,-0.17061,-0.00803,
    -0.14352,-0.30583,-0.2244,0.045758,0.208499,-0.00794,-0.08932,-0.00788,0.346171,-0.55197,
    -0.00783,-0.11649,0.209957,0.047011,0.40302,0.185253,0.351572,0.021206,-0.47564,-0.17191,
    0.048362,-0.72111,-0.19879,0.266892,-0.28129,-0.11672,0.10259,0.130458,0.241257,0.408604,
    0.188379,-0.58863,-0.03397,-0.14468,0.049125,0.104826,0.077391,-0.14492,-0.03372,0.35661,
    -0.42396,0.049973,0.30174,-0.11731,0.247001,-0.03316,0.361129,0.136667,-0.1176,-0.14585,
    -0.20228,0.08017,-0.54066,-0.20184,-0.03308,0.729776,-0.08932,0.280049,-0.57207,-0.96368,
    -0.06124,-0.53775,0.359109,0.276515,0.193005,-0.65318,-0.65002,0.162603,0.388274,-0.03298,
    -0.08932,0.136348,0.30683,0.393839,-0.37382,0.166687,0.195904,0.196719,-0.17522,0.082549,
    -0.11799,-0.03198,-0.06064,-0.03193,-0.08932,0.313325,-0.49197,0.025555,0.198448,-0.29085,
    -0.11808,0.227453,0.141692,-0.11823,-0.06042,-0.06041,-0.0604,-0.06039,-0.205,-0.69442,
    -0.06059,0.025687,0.169932,0.112786,-0.0315,0.113312,-0.58072,-0.37726,0.314021,-0.37759,
    0.343392,0.055333,-0.00243,0.317183,0.027126,0.290074,-0.08932,-0.00156,-0.14784,-0.08932,
    -0.49796,-0.08932,-0.14756,0.582508,-0.84845,0.318723,0.320395,0.410459,0.028637,-0.20728,
    0.590852,0.059156,-0.26747,0.237524,-0.2083,-0.14876,-0.23776,0.029408,-0.26736,0.118424,
    -0.6522,-0.26642,-0.05983,-0.35446,-0.14815,0.14618,-0.05985,-0.53056,-0.11867,-0.26522,
    0.115927,-0.47016,-0.23541,-0.38087,-0.14753,-0.69846,-0.08932,0.287326,-0.03125,-0.61076,
    0.403076,0.844174,-0.26502,-0.26471,-0.78781,-0.08932,-0.35,-0.11824,-0.0604,-0.40701,-0.40601,
    -0.06057,-0.4051,-0.2039,-0.3752,0.196556,-0.11795,-0.2609,-0.43159,-0.06084,-0.2885,
    -0.62795,0.705459,0.110364,-0.71557,-0.06094,-0.1177,-0.06094,-0.08932,-0.74005,-0.62371,
    -0.11737,-0.14539,-0.48093,-0.39594,0.049934,-0.42321,-0.20037,-0.36641,-0.06165,-0.22762,
    0.27064,-0.25562,-0.06163,0.299244,0.300759,-0.14514,-0.5903,-0.56017,0.104289,-0.117,
    -0.3105,0.353524,-0.00607,0.523323,-0.28466,0.357731,-0.00528,-0.42508,-0.14517,-0.86797,
    0.021543,0.494743,0.217984,-0.31291,0.386403,-1.03852,0.691714,-0.64783,-0.39518,-0.61542,
    -0.42017,-0.33674,-1.26308,-0.1707,-0.22481,-0.5486,-0.19708,-0.00851,-0.06237,-0.51966,
    0.314067,-0.00845,-0.70769,-0.33025,-0.32967,0.204521,-1.10078,-0.1952,-0.32713,-0.32657,
    0.438662,-0.27443,0.069324-0.14223,0.680627,-0.91215,-0.2478,0.572669,0.791271,0.071622,
    -0.43771,-0.11607,1.311186,0.673001,-0.19858,-0.17119,-0.55197,-0.44167,-0.33253,0.316358,
    -1.35499,0.151803,0.367729,-1.02804,0.28511,1.015364,0.45401,-0.84916,-1.45865,-1.01833,
    0.26577,-1.961,-0.27197,0.171707,-1.85103,0.270853,1.783782,2.16811,-3.62531,-4.12899,
    7.486337,-10.1731,-1.89421,0.125552,-2.59122,-1.0174,-0.02001,-0.57347,-0.73123,1.732474,
    0.775472,0.145696,1.763023,-0.56756,-1.7923,-2.29282,0.324472,-1.12061,1.241889,1.823496,
    -0.86352,0.944279,-0.84209,-0.27663,4.87394,-2.0368,-2.11779,-0.77208,0.216165,-3.44495,
    -0.0438,1.401702,3.390926,-4.83265,-2.81318,0.378125,-1.24284,-1.4475,-0.5452,0.75897,
    -0.3511,0.281736,1.120681,-1.08041,0.724853,0.420149,0.601489,1.148178,-0.08932,-0.08932,
    -0.38321,0.772166,1.055641,0.072024,1.327823,0.992787,1.100318,1.308628,-0.62188,-0.90683,
    1.06676,-2.10378,-0.82262,0.004995,1.264457,0.678805,0.00711,1.539893,0.64887,0.753937,
    -0.98196,2.23317,-0.79426,-1.26155,-0.68265,0.826866,-2.78887,-0.83712,0.223592,-0.47431,
    -0.85488,-0.42242,2.023095,-0.74221,-0.73797,-0.9002,-2.96907,0.095478,-0.6427,-1.95726,
    -1.78997,-1.01689,1.193881,0.670628,0.789557,0.114598,-1.97659,0.804938,0.496265,-2.01332,
    0.332573,-2.05039,-0.95833,0.344236,-0.58775,-1.24988,0.747503,-0.32606,-0.23968,1.033073,
    -0.08932,2.394049,1.571026,-1.86087,-0.46647,0.132358,-0.5101,0.887812,-0.37902,-0.84304,
    1.491123,0.834836,1.785019,1.164467,-2.89992,1.260041,-1.4614,0.754565,0.139983,-0.8896,
    1.33272,1.282952,0.780965,-1.28694,0.825211,0.264661,1.411242,1.312344,1.20892,0.976496,
    4.031302,-1.66091,-1.46025,0.492296,1.366828,-1.34238,-1.42722,1.121591,-1.75057,0.185644,
    0.714905,3.219945,3.360211,-2.09393,1.591838,1.483805,0.513915,0.048286,1.409041,0.190633,
    2.466183,2.828889,0.029197,0.148137,1.80064,-1.89009,0.05951,0.448314,0.210618,-1.19464,
    2.133659,2.526614,-3.58232,-1.34607,0.626885,-1.30976,-0.08932,3.556209,-2.93283,2.020461,
    1.074907,-2.43449,-2.32189,0.73835,-1.47475,0.057155,1.446897,0.388148,-0.80468,1.135854,
    1.76174,1.765445,0.505266,2.196492,0.943612,2.474078,-4.93029,-1.84982,2.084334,-4.45138,
    1.632738,1.505209,0.924307,-0.78729,0.863671,0.358535,0.489458,-0.34697,1.4993,1.69114,
    -2.3913,0.628462,-0.15479,0.074425,-0.1876,2.162429,1.225378,-0.83307,0.451037,-0.19086,
    -0.22455,2.096557,2.074811,-2.42596,0.014144,-0.26171,1.823447,0.121574,1.043033,0.696677,
    0.378053,0.70663,0.092463,0.567855,0.756738,0.838501,1.375493,-0.08932,-0.20277,0.289323,
    0.138555,-0.05129,0.981434,-0.24299,0.179754,0.373642,-0.43675,1.503254,-1.2958,1.588112,
    -1.65063,0.376517,-0.16711,-0.67083,0.881748,2.560625,1.935039,-1.55098,1.372336,-2.31388,
    1.808575,-8.60344,-0.98385,-0.86552,2.031688,2.039227,1.186551,3.55609,-0.89318,0.472703,
    1.615957,-0.00739,-1.99723,-4.0315,0.453315,-6.23177,-0.78137,0.895548,2.134639,-1.90924,
    1.655664,-0.94706,2.430082,2.300002,-0.24521,-2.20845,2.772436,3.83118,-1.42673,-0.77137,
    1.807809,2.7627,0.79508,-1.76721,0.620345,-0.08932,2.327586,-0.04639,1.993183,0.394089,
    2.137942,-0.71797,-1.24643,0.220893,-0.26671,2.924173,1.661068,1.408833,-0.04214,0.668616,
    3.690198,-2.52847,-0.18566,-0.66542,0.390525,1.608512,0.745858,2.966263,1.345129,-0.14091,
    -0.34685,0.013612,0.065277,0.947494,0.223831,1.703975,1.304055,-0.08932,2.479827,1.19239,
    -1.37103,-0.80657,1.573688,-0.70231,-0.03375,0.188997,0.413613,0.078887,0.022974,0.756945,
    -1.16006,-1.42564,0.632302,0.862266,1.612364,-1.67846,-0.033,0.079835,1.33156,-0.03206,
    -1.90549,0.361637,-1.21293,-0.53525,-0.42247,1.756884,0.704333,0.882948,0.371508,-0.89439,
    -0.31816,-0.37463,0.539431,0.082844,0.428963,0.60593,0.435308,0.203333,-0.90661,0.143507,
    0.027296,0.085861,-0.5558,0.027092,0.085554,1.025386,-1.14567,0.0275,0.261966,0.912749,
    2.06645,-0.33116,0.273655,-0.33145,0.943808,-0.02822,-0.15043,-0.57683,0.94949,0.89834,
    -2.1158,0.459293,0.955234,-0.45924,1.398616,-0.5256,-1.81563,0.400875,-0.76273,-0.75822,
    -2.60292,-0.73734,-0.90804,-0.55417,0.901074,0.556789,0.62032,0.148349,-2.26652,-0.26379,
    0.20163,0.553767,-3.49121,-0.08932,-0.03262,-1.04897,-0.53775,1.432263,-0.65557,0.023672,
    -0.25877,-0.14574,0.874131,-0.48716,0.024185,0.708858,0.139904,-0.43296,-0.48875,-0.26002,
    1.399824,1.305239,0.262378,-0.90804,0.085554,1.202534,-0.56102,0.382377,-4.08638,-1.7227,
    -0.31254,-1.58313,1.85141,-0.42471,2.453906,2.637856,1.272855,0.388862,0.210708,-2.99151,
    -2.22628,2.457112,-0.38198,-1.42461,-2.36992,4.171541,3.564873,0.645525,0.527203,1.846424,
    2.659333,1.412502,-0.08932,0.108241,0.43943,1.851187,1.820695,0.393937,3.002983,0.555149,
    0.559329,1.73488,-2.70575,-0.51881,0.411932,0.414457,-0.01715,0.999618,1.455394,0.505475,
    1.716596,1.133,-0.39631,1.14431,0.689498,0.61684,-0.32526,0.14662,0.860052,0.548639,
    0.391834,-0.33019,0.071192,0.636195,-0.00838,0.723691,0.074077,-0.08932,-0.49732,-1.14239,
    -0.16987,-0.00877,-0.57164,0.473607,0.395723,0.316687,0.809744,-0.00719,-0.25353,0.734406,
    -0.00657,-0.08932,-1.64993,0.07381,0.565881,0.487523,1.074448,-1.17041,-1.89265,-0.97894,
    1.861959,1.566345,0.412771,-0.00539,0.923345,0.676965,1.374329,0.084289,-1.29832,-0.60302,
    0.338578,-0.08932,0.426586,0.515954,1.39607,1.329141,1.078181,-0.08932,-1.07809,1.715781,
    0.092993,-0.08932,0.367926,0.831495,0.374286,-0.46038,0.281736,-0.18222,0.562753,-0.27606,
    0.847892,-0.08932,-0.37141,1.043835,-0.18424,0.38619,-0.18461,-1.41384,-0.37088,-0.08932))

    list(mu=0,itau1=0.02,phi1=0.975)

    1. 您好,请问你这个问题解决了吗?我最近在做用SV预测波动率的实证,遇到了同样的问题,难道theta[j]需要初始化吗?

  19. 当我运行compile这一步是,软件界面上说educational version cannot do this model,是什么意思啊!有知道的麻烦告诉一声,有急用!谢谢

  20. 请问用openbugs时,前面模型都没有问题,但是在load inits这步,总显示unable to generate initial values for node [011483COH] of type Grapht.Mixing,是什么意思?难道初值设定有问题吗?谢谢

  21. 请教大家一下,如果炫耀载入的数据非常多,如何载入那么多数据啊,有没有简便方法?

  22. 等WinBUGS的update时,顺便把评论都看了一遍,觉得很多朋友对Bayes分析还是很有兴趣的,但椰丝儿们留错了地方,cos论坛的椰子板块会更靠谱。看到后来发现很多评论挂着一大长串代码的都在问一件事情:SV模型怎么做?这让我觉得很惊讶,不过想来也合理。

    国内的金融方面最近几年一直很热,学生极多。很多学校会对金融的学生要求论文,而高质量的也是学生保研,各种评优的重要依据。最快的方法,就是用各种统计和计量方法来让自己的文章升级(记得有个神文写的是如何拿计量经济五大杀器力克CSSCI)。我接触的一些研究生,甚至是本科生,都在搞一些很高级的处理方法和模型(高级的意思是读一两本书都搞不清楚)。金融方面比较好用的,无非就是波动模型,因此在GARCH模型已经搞烂掉的现在,SV模型凭借它”拒人于千里之外“的难理解的特点越来越受到学生的青睐。这里的MCMC之所以被提如此之多,就是因为目前SV估计比较好的方法,就是它了。

    有时候一想起这些就会有些郁闷,但觉得孩子们本身没错儿,靠自己的努力和智力,迅速学一些别人学不会的方法来争取自己的前程是很棒的。令我郁闷的是这种特别的规则而引发的浮躁心态以及对思考问题的态度。就拿SV来说,有个孩子学经济的,问想拿MCMC估计SV,但理解不了那个超长的函数,我就说这个叫做后验概率密度,他说什么叫后验密度?我说是通过贝叶斯原理得到的,他说贝叶斯不是统计学的吗?我要估计的是SV模型……尽最快速度达到目标,却忽略了积累与思索,这便得不偿失了。如果有负责的老师在这方面加以引导和帮助,我想这便是孩子们在大学,甚至是近十年的一大幸事。可惜的是这可遇而不可求,因为老师同样面临和孩子们一样的境况,只不过从保研变成了评职称而已,并且他们只能靠自己。

    如果真的不想要培养研究人才,拜托请不要变相用这种方式折磨孩子们的思维了。至少在绝大多数工作中,老板不会让你没事就去研究MCMC……而在看过评论之后唯一欣慰的是,COS的编辑们很热心负责的,他们抱着对统计认真的态度,在帮助在这条路上迷乱的人,尽管这很有限,但也让很多人满载而归了。

    深夜吐槽 尽请见谅~那我也来做点有用的吧。SV不多说,R上面也有包可以做,MCMC方法cos上也有文章介绍。如果你确实初学,没接触R也不想学非要用WinBUGS,这里有个paper
    http://www.mysmu.edu/faculty/yujun/Research/YuEJ2000.pdf
    我觉得这个东西,结合上面文章的内容,绝大多SV问题就应该解决了。另外有其他问题可以到cos论坛金融区去请教版主,版主在这方面可是很强的。

  23. 请问一个问题,我在看各种winBUGS的例题时,一般标准差都会有这么一句sigma<-sqrt(1/tau),其中tau是方差,但是不应该是sigma<-sqrt(tau)吗?

  24. 请问谢老师,用openbugs时,在load inits这步,总显示unable to generate initial values for node要怎么解决?

  25. MODEL:
    Y=a*T1*P1*W*E*P2/(1+P2/b)+c*T2+d

    Input data:
    T1= [0.9902 0.9696 0.9851 0.9341 0.9787 0.9502 0.8941 0.9987 0.9987 0.9897 0.9946 0.7887]
    T2= [18.261 22.945 22.074 24.285 22.474 23.742 25.379 20.631 19.377 21.735 18.711 11.493]
    P1= [0.5961 0.5777 0.5989 0.6213 0.6258 0.6478 0.6416 0.6554 0.6514 0.6553 1.0000 1.0000]
    P2= [19265.5 18877 18052.5 15476 20851.5 21657 22012 6481 21738 20884 17597.5 18969]
    W= [0.8755 0.8492 0.8791 0.9119 0.9176 0.9511 0.9417 0.9625 0.9552 0.9615 0.9525 0.9180]
    E= [0.2927 0.2817 0.3708 0.4006 0.3684 0.4377 0.4311 0.4939 0.4847 0.4762 0.4602 0.4043]

    Yo= [-6.23 -7.24 2.21 6.1 15.18 16.38 13.63 2.53 10.51 4.08 0.73 3.51]

    Cost function:
    J(p) = 1/2*[(Yo-Y(p))’Co ^(-1) (Yo-Y(p))+(p-pb) ‘ Pb^(-1) (p-pb)]
    where Yo represents the data vector, Y(p) is the model output vector, Co and Pb the error covariance matrix on data and model parameters, respectively. p is the parameter vector. pb is a priori values of p.
    pb = [0.06, 1600, 0.05, 1];

    My question is how to estimate the parameter (a, b, c, d) using WinBUGS?

  26. 请问 做SV-t模型的参数后想导出波动率theta[i](i=1,……,n)的值,但是n比较大,有什么办法吗?谢谢~

  27. 请问 做SV-t模型的参数估计后想导出波动率theta[i](i=1,……,n)的值,但是n比较大,有什么办法吗?谢谢~

发表评论

电子邮件地址不会被公开。 必填项已用*标注