winbugs操作方法(已翻译)

1、例(引自于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 5

2、51.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 inB

3、UGS14,会出现两个窗口 ，关闭其中一个(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)。步

4、骤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的形状，以

5、从上到下、从左至右的顺序，利用Doodle菜单构建Doodle模型。(1)打开Doodle菜单选择New (新建)命令，可以打开一个名为“ New Doodle "的对话框。键入200、150和 20,作为编辑窗口的显示宽度、显示高度和节点的宽度。这三个选项控制编辑窗口和节点的大小，输入数值的大小应与模型相适合。点击0K, 个无标题名（untitledl）的DoodleBUGS编辑窗口就产生了。（注:可以按默认值打开 ,窗口大小可以按调节 windows 窗口的一般方法调节 ） 。在此DoodleBUGS编辑窗口（ untitledl）中的空白处单击（左键，下同）鼠标，会自动产生一个

6、椭 圆形节点 （ node） （注 :在此窗口中 ,不能按习惯任意点击 ,因为一点就会出现一个节点。 多余节点 删除与一般删除不同，先用鼠标将需删除的节点点击选中，然后按住Ctrl键不放，再按下delete键或backspace键删除），窗口左上方也会同时出现七个蓝色的标题，分别为name （名称），type（类型 ） ,density （ 密度 ） , mean （均值） , p recision （ 精度 ） ,lower bound （ 下界） , upper bound （ 上界 ） , 光标在第一个名称栏（name）内闪烁，节点可以用鼠标点中后随意拖动。在名称栏（name）内输入&q

7、uot;alpha”（表示）作为此节点的名称。然后，name右侧的是类型栏（type）,单击（type）不放它会 自动出现下拉菜单 ,菜单中有三个选项 :随机、 逻辑和常量 （ stochastic,logical, and constant） , 第一 项是默认值。由于alpha是随机变量，故选择默认值;再往右是密度栏（density）,单击此栏会出现 16种分布的下拉菜单栏（默认值是正态分布dnorm）。因为假设alpha服从Exponential指数分布， 所以选择dexp （表示Exponential分布）密度。接下来，在scale框中输入数值 1. 0 （ a Exponential

8、 （ 1. 0）。对于下界和上界（lower bound and upper bound）栏，不用输入数 值,因为不需要对此节点的数值范围加以限制。这样，就完成了对模型中参数a的节点属性定义。在窗口中的其它位置重新再点击鼠标，生成一个新的节点，命名为beta （表示3 ）。这个新的节点的颜色相比会显得“高亮” （ highlighted） , 称之为选中状态。（注:节点只有处于选中状态时 ,才能对节点的属性进行修改和编辑 ） 。参照以上方法对此节点按假设条件 （3 Gamma （0. 1, 1.0）进行属性设置。类型栏（type）中默认stochastic,在密度（densi2ty）类型中按要

9、求选择 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节点（内部），这时会在两个节点之间产生

10、一 个实线箭头 （表示随机节点到随机节点） ,箭头的方向从 alpha 节点指向 theta i 节点将父节点与子节点进行连接这时 ，在形状框（shape）中就会自动出现alpha。同样操作，建立从beta节点指向theta i 节点的实线箭头，这时，在尺度框（scale）中就会自动出现beta。为了以 后的应用 ,先将怎样进行箭头移动和删除的方法进行说明 ,想移动箭头 ,不要去点击箭头 ,通过 移动节点，箭头会随着自动调整。要想删除箭头，先选中子节点（箭头所指）撚后按住Ctrl键不放，再次点击父节点 （箭尾 ） ，箭头删除要注意次序 ，否则会加上一个新箭头。按照以上方式继续进行下去，创建节点

11、 t i 。对这个节点 ， 选择它的类型 （ type ） 为常数（ con2stant） ，这时 ，你会发现节点 t i 的形状会变成矩形。（注意 :属性栏也随之发生变化 ） 。接着创建并定义一个新的节点，命名为lambda i ,选择它的类型（type）为逻辑logical,在value栏中，输入逻辑表达式 theta i *t i 。按前述方法创建从父节点 theta i 到子节点 lambda i 的箭 头，注意，这次显示的箭头与以前有区别，呈双线空心箭头 ，这种箭头表示随机变量到逻辑变量之间的联系 ，区别于两个随机变量之间联系的实线箭头。同理，也可建立从父节点 t i 到子节点lam

12、bda i 的双线空心箭头。最后,根据条件Xipoisson （0 i ti）,创建节点x i , density中选择密度为dpois （表示poisson分布），创建从父节点lambda i 到它的箭头，父节点的名称 lambda i会自动出现在 mear栏中。建立矩形“平板”（p lates）。按下Ctrl键不放，在Doodle窗口内的空白处点击鼠标 ，就可创建一个矩形“平板”。用鼠标点击“平板”右侧或下侧较粗的边线，就可选中激活它 ，并且可以对它进行移动。 想对“平板” 的大小进行调整 ，将鼠标点击 “平板” 的右下角顶点进行拖动 （如同对一般窗口调整大小 ） 。想删除“平板” ,选中

13、（点击“平板”右侧或下侧较粗的边线） ,按住Ctrl键不放，再按下delete键或backspace键删除。此例中，将创建的"平板”的大小调整到如 图1所示 ,包含下半部分的四个节点。当一个“平板”被选中后（点击“平板”右侧或下侧较粗的边线 ） ， 窗口的左上方会出现三个蓝色的标题:指标、起始和结束 （ index， from， up to） 。首先在指标文本框（index）中输入指标名称，它表示一列整数序列数值，从from中的数值到up to 中的数值，在此例中，index中输入i, from中输入1, up to中输入N。这样就完成了例子的全部 Doodle模型的构建。这时比较明

14、智的选择应该是把刚才辛苦的工作保存下来，在File菜单中，选择Save as命令撚后在对话框中，输入便于记忆的文件名（用英文）,WinBUGS会自动在文件 名之后加上后缀名.ode,选择合适的目录保存，便于以后使用。2 对Doodle模型进行检验（cheek）构建完Doodle模型后,下一步是对此模型进行检验，检验Win BUGS对它是否识别,检验模型在程序语法上是否正确。在Modle菜单中选Speeifi2eati on命令 ，会出现一个标题为“ Specification Tool ”的对话窗口 （用完后不要关闭 ，以下步骤还需使 用）。在窗口中点击“ cheek model ”按钮,如果

15、你建立的Doodle模型正确,在W in BUGS主窗口 的左下角的状态栏 （status line）上就会出现一条信息" model is syntactically correct ”。否则， 你必须对你的Doodle模型的每一个节点、每一个箭头和每一块平板进行检查，耐心检查它们的名称和属性，改正错误的地方。然后，重新对模型进行检验？?,直到W in BUGS宣布“ model is syntacticallycorrect ”。步骤 2:数据输入、模型编译和初始值设定 1 输入数据 （ data loaded） 模型检验合格后 ，需要对数据进行定义和输入。 例中的数据可以按 以

16、下方式进行定义 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命令。出现的窗口可以对文本进行编辑。将数据在英文输入状态下按 以上形式分两行输入。 输入完毕 ,检查确定无误后 ,选中输入内容 （将光标点中以上输入内容中 任意位置即可），然后在“Specific

17、ation Tool对话窗口中点击“load data"按钮。信息"data loaded” 将会在状态栏显示。 2 对模型进行编译 （ comp ile）继续在“ Specification Tool ”的对话窗口中通过“ comp ile ”按钮完成。在对模型进行编译之 前,注意在“ comp ile ”按钮右边上的文本框中的正整数应大于或等于默认值 1 （这个数字表示 使用者希望同时模拟互相独立的链的个数）。模型的编译需要运行 Gibbs samp ing,模型需要对数据的完整性和一致性进行检验。当模型编译成功后,信息“ model comp iled ”会出现在状态

18、栏中（如果你输入数据有误，需要从模型检验（check）步骤开始重做）。 3 对模型进行初始值设定 （ load inits） 在编译好的模型真正运行之前,还需要对模型进行初始值设定 （ load inits） 。以下的语句可以执行这项任务。list （ alpha = 1, beta = 1 ） 注意 ,这只是模拟的初始值 ,并不需要一定要与我们实际期望的参数值 很接近。 将这条命令继续输入到数据窗口中（可换行 ,在刚才输入的两行数据下输入） , 选中这个语句，按“ Specificati on Tool的对话窗口中的“ load in its ”按钮，状态栏中会出现信息“ thisiniti

19、al 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，然后按下“s

20、et”；输入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 中已输入的控制变量 ) ,也可以输入万能替代符