winbugs操作方法已翻译

例(引自于W in BUGS软件帮助手册中的 Volum I) : George et al (1993)讨论了分层模型的贝 叶斯分析(Bayesian analysis of hierarchical mod 2els)。其中，第一层采用了共轭先验分布。该例 考虑了相关的10个发电站的水泵，假设发生故障的水泵个数服从 Poisson分布，即为Poisson ( 0 i ti ) , i = 1,2, ?10.其中,0 i表示水泵i的发生故障率，ti表示水泵运行时间的长度(单位:千小时)，数据见下表：水泵123455?S910ri15. ： 35 24 3L 4 L 05105 2 11Q 551.5143 I? i1422故障率的共轭先验分布假设为0Gamma ( a , 3 ) , i = 1,2, ?10.George et al (1993)对超参数 a 和 3 先验假设为 a Exponential (1.0)3 Gamma (0. 1,1. 0)他们给出了 3的后验分布，但3的标准后验分布无法给出。因而 他们使用Gibbs samp ler模拟 得到a的后验密度。以上表达式可用贝叶斯图建模方法表示成如下有向关系图(图1)，在Win BUGS中称作Doodle模型。(注意:a用alpha表示,xi 用 x i 表示,入用 lambda i 表示,其余类 似。启动W inBUGS14,会出现两个窗口 ，关闭其中一个(L icence Agreement窗口 )，你会看到如下Win 2BUGS主窗口(图2)窗口简单明了，W in BUGS主窗口和W in dows常用窗口结构类似，关闭、 最小化等基本操作相同。Win BUGS主窗口最上面一行为标题栏(headi ng line ),下一行是菜单栏(menuheading line )有 File, Tools, Edit, Attributesnfo, Model, Inferenee, Op tions, Doodle,Map,Window, Help12个菜单，最下面一行是状态栏(status line)。步骤1: Doodle模型的建立和检验1 建立Doodle模型在W in BUGS的使用中，Doodle模型非常特别，它以节点(nodes)、箭头(edges)和平板(p lates)等图形方式出现(对应图1:包含alpha的小椭圆形和包含t i 的小矩形都为节点；箭头分 实线箭头和双线空心箭头；右侧和下侧边线较粗的大矩形称为“平板”)，可以用图形的方式来构建模型。模型的每个节点，都含有特定的属性，如名称、类型、分布或逻辑函数的定义等(name, type, distribution, orlogical function definition)。m I mdb etv以下我们按照图1的形状，以从上到下、从左至右的顺序，利用Doodle菜单构建Doodle模型。(1)打开Doodle菜单选择New (新建)命令，可以打开一个名为“ New Doodle 的对话框。键入200、150和 20,作为编辑窗口的显示宽度、显示高度和节点的宽度。这三个选项控制编辑窗口和节点的大小，输入数值的大小应与模型相适合。点击0K, 个无标题名（untitledl）的DoodleBUGS编辑窗口就产生了。（注:可以按默认值打开 ,窗口大小可以按调节 windows 窗口的一般方法调节 ） 。在此DoodleBUGS编辑窗口（ untitledl）中的空白处单击（左键，下同）鼠标，会自动产生一个椭 圆形节点 （ node） （注 :在此窗口中 ,不能按习惯任意点击 ,因为一点就会出现一个节点。 多余节点 删除与一般删除不同，先用鼠标将需删除的节点点击选中，然后按住Ctrl键不放，再按下delete键或backspace键删除），窗口左上方也会同时出现七个蓝色的标题，分别为name （名称），type（类型 ） ,density （ 密度 ） , mean （均值） , p recision （ 精度 ） ,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 （ a Exponential （ 1. 0）。对于下界和上界（lower bound and upper bound）栏，不用输入数 值,因为不需要对此节点的数值范围加以限制。这样，就完成了对模型中参数a的节点属性定义。在窗口中的其它位置重新再点击鼠标，生成一个新的节点，命名为beta （表示3 ）。这个新的节点的颜色相比会显得“高亮” （ highlighted） , 称之为选中状态。（注:节点只有处于选中状态时 ,才能对节点的属性进行修改和编辑 ） 。参照以上方法对此节点按假设条件 （3 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（ 0 i Gamma （ a , 3 ），通过下面箭头关系的操作，程序会自动生成shape和scale中的值。保持节点 theta i 处于选中状态，按住Ctrl键不放，点击alpha节点（内部），这时会在两个节点之间产生一 个实线箭头 （表示随机节点到随机节点） ,箭头的方向从 alpha 节点指向 theta i 节点将父节点与子节点进行连接这时 ，在形状框（shape）中就会自动出现alpha。同样操作，建立从beta节点指向theta i 节点的实线箭头，这时，在尺度框（scale）中就会自动出现beta。为了以 后的应用 ,先将怎样进行箭头移动和删除的方法进行说明 ,想移动箭头 ,不要去点击箭头 ,通过 移动节点，箭头会随着自动调整。要想删除箭头，先选中子节点（箭头所指）撚后按住Ctrl键不放，再次点击父节点 （箭尾 ） ，箭头删除要注意次序 ，否则会加上一个新箭头。按照以上方式继续进行下去，创建节点 t i 。对这个节点 ， 选择它的类型 （ type ） 为常数（ con2stant） ，这时 ，你会发现节点 t i 的形状会变成矩形。（注意 :属性栏也随之发生变化 ） 。接着创建并定义一个新的节点，命名为lambda i ,选择它的类型（type）为逻辑logical,在value栏中，输入逻辑表达式 theta i *t i 。按前述方法创建从父节点 theta i 到子节点 lambda i 的箭 头，注意，这次显示的箭头与以前有区别，呈双线空心箭头 ，这种箭头表示随机变量到逻辑变量之间的联系 ，区别于两个随机变量之间联系的实线箭头。同理，也可建立从父节点 t i 到子节点lambda i 的双线空心箭头。最后,根据条件Xipoisson （0 i ti）,创建节点x i , density中选择密度为dpois （表示poisson分布），创建从父节点lambda i 到它的箭头，父节点的名称 lambda i会自动出现在 mear栏中。建立矩形“平板”（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会自动在文件 名之后加上后缀名.ode,选择合适的目录保存，便于以后使用。2 对Doodle模型进行检验（cheek）构建完Doodle模型后,下一步是对此模型进行检验，检验Win BUGS对它是否识别,检验模型在程序语法上是否正确。在Modle菜单中选Speeifi2eati on命令 ，会出现一个标题为“ Specification Tool ”的对话窗口 （用完后不要关闭 ，以下步骤还需使 用）。在窗口中点击“ cheek model ”按钮,如果你建立的Doodle模型正确,在W in BUGS主窗口 的左下角的状态栏 （status line）上就会出现一条信息 model is syntactically correct ”。否则， 你必须对你的Doodle模型的每一个节点、每一个箭头和每一块平板进行检查，耐心检查它们的名称和属性，改正错误的地方。然后，重新对模型进行检验？?,直到W in BUGS宣布“ 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 ） 注意 ,这只是模拟的初始值 ,并不需要一定要与我们实际期望的参数值 很接近。 将这条命令继续输入到数据窗口中（可换行 ,在刚才输入的两行数据下输入） , 选中这个语句，按“ Specificati on Tool的对话窗口中的“ load in its ”按钮，状态栏中会出现信息“ thisinitial value does not correspond to astochastic node ”。这时 , “Specification Tool ”的对话窗口中 的“ gen inits ”按钮将被激活，按下“ gen in its ”按钮,状态栏上会出现信息 “ ini tial values gen erated, model initialized ”。步骤 3:变量监控和模型迭代1 对变量进行监控设置（Monitor）打开Inference菜单，选择Samp les命令。会产生一个名为“Samp le Monitor Tool的对话窗口，在node框中输入alpha，然后按下“set”；输入beta,按下“set”;输入theta 1 ,按下“ set”，再输入theta 2 , ?W in BUGS可以同时控制多个变量，但必须分别 输入变量名和依次进行set”操作。进行以上操作后，WinBUGS会在运行时，对输入的控制变 量的运行过程及结果进行存储。 2 模型迭代运行 ( iterations, or updates) 当以上步骤完成之后 ,就可以对模型进行迭代运行了。在Model菜单中选择Update命令，会出现一个Update Tool”的对话窗口(图8)，输入需要迭代(iterations/updates)的次数，按update”按钮开始执行。迭代的速度非常的快，对此模型进行10000次迭代运算，在 150MHz Pentium lap top的 PC上需要 5秒，而在 a dual 200MHz Pentium Pro 的机器上只需 2秒钟。步骤 4:模拟结果输出和收敛性判别再次选中“Samp le Monitor Tool的对话窗口，在node框中，可以按需要输入你想观测的变量(节点) 的名称 (步骤3 中已输入的控制变量 ) ,也可以输入万能替代符“*”来表示“输入的所有需要控制的变量”,然后 ,点击“ stats” ,会同时产生所有控制的变量的统计量模拟估计结果;“density ”按钮可以对连续变量画出光滑的核密度估计( kernel density estimate) , 对离散变量出直方图(histogram);按trace”按钮，被选择作为控制的变量的样本路径会由动态轨迹图(ady2namic trace p lot) 显示出来。注：单击“ coda”按钮可以打开两个输出窗口。一个窗口，我们将其称为“ ind”窗口，它包含了控制变量的名称和列数的范围 ( a range of row num 2bers) 。这些列数相应的输入行出现在 另一个窗口“out ”窗口中，这个窗口包含控制变量的模拟值。这两个窗口给出了模拟数据的ASCII表示方式，这种表示方式可以复制到其它应用软件中使用，例如Excel或S - PLUS。另外，软件会自动产生模型的源程序：(在Doodle窗口激活的条件下，在Doodle菜单中,点击WriteCode命令)(2)收敛性的判别 ：对于任何以 MCMC 为基础的完全概率模型分析,关于 MCMC 模拟的监测和收敛性的判断都是非常重要的一步。一个 MCMC 模拟称为收敛的 ,是指模拟结果来自于真的 Markovchain 的稳定或目标分布。能从理论上证明收敛性的情形非常少,大多数情况 ,我们只能满足于有效收敛 ( effective convergence) 。Log (日志)窗口会记录以上操作的信息丄og窗口可以通过选择Info菜单中的Open Log命令打开。