winbugs操作方法(已翻译).docx

例 (引自于W inBUGS软件帮助手册中的Volum I) : George et al (1993) 讨论了分层模型的贝叶斯分析(Bayesian analysis of hierarchical mod2els) 。其中,第一层采用了共轭先验分布。该例考虑了相关的10个发电站的水泵,假设发生故障的水泵个数服从Poisson分布,即xi Poisson (i ti ) , i = 1, 2, 10. 其中,i 表示水泵i的发生故障率, ti 表示水泵运行时间的长度(单位:千小时) ,数据见下表:故障率的共轭先验分布假设为i Gamma (,) , i = 1, 2, 10.George et al (1993) 对超参数和先验假设为Exponential (1. 0) Gamma (0. 1, 1. 0)他们给出了的后验分布,但的标准后验分布无法给出。因而,他们使用Gibbs samp ler模拟得到的后验密度。以上表达式可用贝叶斯图建模方法表示成如下有向关系图(图1) ,在WinBUGS中称作Doodle模型。(注意:用alpha表示, xi 用x i 表示,i用lambda i 表示,其余类似。启动W inBUGS14,会出现两个窗口,关闭其中一个(L icence Agreement窗口) ,你会看到如下Win2BUGS主窗口(图2) :窗口简单明了,W inBUGS主窗口和W indows常用窗口结构类似,关闭、最小化等基本操作相同。WinBUGS主窗口最上面一行为标题栏( heading line ) , 下一行是菜单栏(menuheading line ) 有File, Tools, Edit, Attributes,Info, Model, Inference, Op tions, Doodle, Map,Window, Help12个菜单,最下面一行是状态栏( status line) 。步骤1: Doodle模型的建立和检验 1 建立Doodle模型在W inBUGS的使用中,Doodle模型非常特别,它以节点( nodes) 、箭头( edges)和平板( p lates)等图形方式出现(对应图1:包含alpha的小椭圆形和包含t i 的小矩形都为节点;箭头分实线箭头和双线空心箭头; 右侧和下侧边线较粗的大矩形称为“平板”) ,可以用图形的方式来构建模型。模型的每个节点,都含有特定的属性,如名称、类型、分布或逻辑函数的定义等( name, type, distribution, orlogical function definition) 。以下我们按照图1的形状,以从上到下、从左至右的顺序,利用Doodle菜单构建Doodle模型。(1)打开Doodle菜单,选择New (新建)命令,可以打开一个名为“New Doodle”的对话框。键入200、150和20,作为编辑窗口的显示宽度、显示高度和节点的宽度。这三个选项控制编辑窗口和节点的大小,输入数值的大小应与模型相适合。点击OK,一个无标题名( untitled1)的DoodleBUGS编辑窗口就产生了。(注:可以按默认值打开,窗口大小可以按调节windows窗口的一般方法调节) 。(2)在此DoodleBUGS编辑窗口( untitled1)中的空白处单击(左键,下同)鼠标,会自动产生一个椭圆形节点( node) (注:在此窗口中,不能按习惯任意点击,因为一点就会出现一个节点。多余节点删除与一般删除不同,先用鼠标将需删除的节点点击选中, 然后按住Ctrl键不放, 再按下delete 键或backspace键删除) ,窗口左上方也会同时出现七个蓝色的标题,分别为name (名称) , type (类型) ,density (密度) , mean (均值) , precision (精度) ,lower bound (下界) , upper bound (上界) ,光标在第一个名称栏( name)内闪烁,节点可以用鼠标点中后随意拖动。在名称栏( name)内,输入“alpha”(表示)作为此节点的名称。然后, name右侧的是类型栏( type) ,单击( type)不放它会自动出现下拉菜单,菜单中有三个选项:随机、逻辑和常量( stochastic,logical, and constant) ,第一项是默认值。由于alpha是随机变量,故选择默认值;再往右是密度栏( density) ,单击此栏会出现16种分布的下拉菜单栏(默认值是正态分布dnorm) 。因为假设alpha服从Exponential指数分布,所以选择dexp (表示Exponential分布)密度。接下来,在scale框中输入数值1. 0 (Exponential ( 1. 0) ) 。对于下界和上界( lower bound and upper bound)栏,不用输入数值,因为不需要对此节点的数值范围加以限制。这样,就完成了对模型中参数的节点属性定义。在窗口中的其它位置重新再点击鼠标,生成一个新的节点,命名为beta (表示) 。这个新的节点的颜色相比会显得“高亮”( highlighted) ,称之为选中状态。(注:节点只有处于选中状态时,才能对节点的属性进行修改和编辑) 。参照以上方法对此节点按假设条件(Gamma (0. 1, 1. 0) )进行属性设置。类型栏( type)中默认stochastic,在密度( densi2ty)类型中按要求选择dgamma分布,在形状参数框( shape)中输入0. 1,在尺度参数框( scale)中输入1. 0。创建一个新的节点theta i , name后输入theta i , type中默认stochastic , density中按要求选择dgamma分布。注意,对于此节点, shape和scale中不用输入值,因为按假设它们分别为alpha和beta(i Gamma (,) ) ,通过下面箭头关系的操作,程序会自动生成shape和scale中的值。保持节点theta i 处于选中状态,按住Ctrl键不放,点击alpha节点(内部) ,这时会在两个节点之间产生一个实线箭头(表示随机节点到随机节点) ,箭头的方向从alpha 节点指向theta i 节点将父节点与子节点进行连接这时,在形状框( shape)中就会自动出现alpha。同样操作,建立从beta节点指向theta i 节点的实线箭头,这时,在尺度框( scale)中就会自动出现beta。为了以后的应用,先将怎样进行箭头移动和删除的方法进行说明,想移动箭头,不要去点击箭头,通过移动节点,箭头会随着自动调整。要想删除箭头,先选中子节点(箭头所指) ,然后按住Ctrl键不放,再次点击父节点(箭尾) ,箭头删除要注意次序,否则会加上一个新箭头。按照以上方式继续进行下去,创建节点t i 。对这个节点, 选择它的类型( type ) 为常数( constant) ,这时,你会发现节点t i 的形状会变成矩形。(注意:属性栏也随之发生变化) 。接着创建并定义一个新的节点,命名为lambda i ,选择它的类型( type)为逻辑logical,在value栏中,输入逻辑表达式theta i * t i 。按前述方法创建从父节点theta i 到子节点lambda i 的箭头,注意,这次显示的箭头与以前有区别,呈双线空心箭头,这种箭头表示随机变量到逻辑变量之间的联系,区别于两个随机变量之间联系的实线箭头。同理,也可建立从父节点t i 到子节点lambda i 的双线空心箭头。最后,根据条件xi poisson (i ti ) ,创建节点x i , density中选择密度为dpois (表示poisson 分布) ,创建从父节点lambda i 到它的箭头,父节点的名称lambda i会自动出现在mean栏中。(3)建立矩形“平板”(p lates)。按下Ctrl键不放,在Doodle窗口内的空白处点击鼠标,就可创建一个矩形“平板”。用鼠标点击“平板”右侧或下侧较粗的边线,就可选中激活它,并且可以对它进行移动。想对“平板”的大小进行调整,将鼠标点击“平板”的右下角顶点进行拖动(如同对一般窗口调整大小) 。想删除“平板”,选中(点击“平板”右侧或下侧较粗的边线) ,按住Ctrl键不放,再按下delete键或backspace键删除。此例中,将创建的“平板”的大小调整到如图1所示,包含下半部分的四个节点。当一个“平板”被选中后(点击“平板”右侧或下侧较粗的边线) ,窗口的左上方会出现三个蓝色的标题:指标、起始和结束( index, from, up to) 。首先在指标文本框( index)中输入指标名称,它表示一列整数序列数值,从from中的数值到up to中的数值,在此例中, index中输入i, from中输入1, up to中输入N。这样就完成了例子的全部Doodle模型的构建。这时比较明智的选择应该是把刚才辛苦的工作保存下来,在File菜单中,选择Save as命令,然后在对话框中,输入便于记忆的文件名(用英文) ,WinBUGS会自动在文件名之后加上后缀名. odc,选择合适的目录保存,便于以后使用。 2 对Doodle模型进行检验( check)构建完Doodle模型后,下一步是对此模型进行检验,检验WinBUGS对它是否识别,检验模型在程序语法上是否正确。在Modle菜单中选Specifi2cation命令,会出现一个标题为“Specification Tool”的对话窗口 (用完后不要关闭,以下步骤还需使用) 。在窗口中点击“check model”按钮,如果你建立的Doodle模型正确,在W inBUGS主窗口的左下角的状态栏( status line)上就会出现一条信息“model is syntactically correct” 。否则,你必须对你的Doodle模型的每一个节点、每一个箭头和每一块平板进行检查,耐心检查它们的名称和属性,改正错误的地方。然后,重新对模型进行检验,直到W inBUGS 宣布“model is syntacticallycorrect”。步骤2:数据输入、模型编译和初始值设定 1 输入数据( data loaded)模型检验合格后,需要对数据进行定义和输入。例中的数据可以按以下方式进行定义list ( t = c ( 94. 3, 15. 7, 62. 9, 126, 5. 24,31.4, 1. 05, 1. 05, 2. 1, 10. 5) ,x = c (5, 1, 5, 14, 3, 19, 1, 1, 4, 22) ,N = 10)这种数据的表示格式被称为S - PLUS格式。观测值t和x被定义成数组(如有缺失数据,用NA表示) 。在File菜单中选择New命令。出现的窗口可以对文本进行编辑。将数据在英文输入状态下按以上形式分两行输入。输入完毕,检查确定无误后,选中输入内容(将光标点中以上输入内容中任意位置即可) ,然后在“Specification Tool”对话窗口中点击“load data”按钮。信息“data loaded”将会在状态栏显示。 2 对模型进行编译( comp ile)继续在“Specification Tool”的对话窗口中通过“comp ile”按钮完成。在对模型进行编译之前,注意在“comp ile”按钮右边上的文本框中的正整数应大于或等于默认值1 (这个数字表示使用者希望同时模拟互相独立的链的个数) 。模型的编译需要运行Gibbs samp ing,模型需要对数据的完整性和一致性进行检验。当模型编译成功后,信息“model comp iled”会出现在状态栏中(如果你输入数据有误,需要从模型检验( check)步骤开始重做) 。 3 对模型进行初始值设定( load inits)在编译好的模型真正运行之前,还需要对模型进行初始值设定( load inits) 。以下的语句可以执行这项任务。list ( alpha = 1, beta = 1 )注意,这只是模拟的初始值,并不需要一定要与我们实际期望的参数值很接近。将这条命令继续输入到数据窗口中(可换行,在刚才输入的两行数据下输入) , 选中这个语句, 按“SpecificationTool”的对话窗口中的“load inits”按钮,状态栏中会出现信息“this initial value does not correspond to astochastic node”。这时,“Specification Tool”的对话窗口中的“gen inits”按钮将被激活,按下“gen inits”按钮,状态栏上会出现信息“initial values generated,model initialized”。步骤3:变量监控和模型迭代 1 对变量进行监控设置(Monitor)打开Inference菜单,选择Samp les命令。会产生一个名为“Samp le Monitor Tool”的对话窗口,在node框中输入alpha,然后按下“set”; 输入beta,按下“set”;输入theta 1 ,按下“set”,再输入theta 2 , W inBUGS可以同时控制多个变量,但必须分别输入变量名和依次进行“set”操作。进行以上操作后,WinBUGS会在运行时,对输入的控制变量的运行过程及结果进行存储。 2 模型迭代运行( iterations, or updates)当以上步骤完成之后,就可以对模型进行迭代运行了。在Model菜单中选择Update命令,会出现一个“Update Tool”的对话窗口(图8) ,输入需要迭代( iterations/updates)的次数,按“update”按钮开始执行。迭代的速度非常的快,对此模型进行10 000次迭代运算,在150MHz Pentium lap top 的PC上需要5秒,而在a dual 200MHz Pentium Pro的机器上只需2秒钟。步骤4:模拟结果输出和收敛性判别再次选中“Samp le Monitor Tool”的对话窗口,在node框中,可以按需要输入你想观测的变量(节点) 的名称(步骤3 中已输入的控制变量) ,也可以输入万能替代符“* ”来表示“输入的所有需要