WinBUGS在统计分析中的应用

WinBUGS在统计分析中的应用（第一部分） By 齐韬 @ 2008/12/08 关键词：MCMC, R, WinBUGS, 空间统计 分类：统计软件 作者信息：Computational Mathematician in Annpro Analytic Technologies, Inc. 版权声明：本文版权归原作者所有，未经许可不得转载。原文可能随时需要修改纰漏，全文复制转载会带来不必要的误导，若您想推荐给朋友阅读，敬请以负责的态度提供原文链接；点此查看如何在学术刊物中引用本文 常规引用方式 齐韬. WinBUGS在统计分析中的应用（第一部分）. 统计之都, 2008.12. URL: http://cos.name/2008/12/statistical-analysis-and-winbugs-part-1/. BibTeX引用 @ARTICLE{, AUTHOR = {齐韬}, TITLE = {WinBUGS在统计分析中的应用（第一部分）}, JOURNAL = {统计之都}, YEAR = {2008}, month = {12}, URL = {http://cos.name/2008/12/statistical-analysis-and-winbugs-part-1/}, } 开篇词 首先非常感谢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，感受一下基本的操作流程，然后再考虑高级的操作。 第一步，打开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 附第二步中的代码如下： #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在统计分析中的应用 第一部分完 1. 胡江堂 on 2008/12/09 at 18:03 WinBugs，头一回听说啊。老齐在annpro，有空分享一下工作心得吧。 回复 2. 刘思喆 on 2008/12/10 at 11:13 土了，还以为是 windows 下的 debug 工具。原来是这 ！！ 回复 3. 谢益辉 on 2008/12/15 at 20:14 呵呵，能写点公式还是写点吧，看代码不知道做的问题是什么…… 回复 4. 谢益辉: 统计之都《本周导读》第三辑 | 统计之都 on 2008/12/15 at 20:23 [...] WinBUGS是贝叶斯统计的有力工具（而不是Windows下面的Debug工具），齐韬加入COS主站之后发 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 了第一篇文章WinBUGS在统计分析中的应用（第一部分），为我们讲述了WinBUGS的一些基本操作。 [...] 回复 5. 齐韬 on 2008/12/16 at 11:45 我会在第二部分中将相关的理论部分加上，看来得要学习一下怎么嵌入LaTex了. 回复 6. 王化儒 on 2008/12/16 at 16:57 嗯，楼主辛苦了！跟着一点一点学，期待着空间数据分析。。。 回复 7. 胡江堂 on 2008/12/18 at 11:09 贴一首诗吧，刚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. 回复 8. xjx on 2008/12/27 at 18:18 不知道在linux下能不能用? 回复 9. hong on 2009/04/25 at 21:55 大家好 麻烦问一下这个迭代次数如何选择.谢谢回答 回复 10. DJ on 2009/05/23 at 17:50 很好的教程，谢谢了。MCMC万岁！ 回复 11. 左伊秩訾 on 2009/05/26 at 22:35 前几天上课时听说了这个软件，真是及时雨！谢谢了！ 回复 12. icwei on 2009/06/01 at 21:35 小弟现正在剑桥mrc-bsu做postdoc，具体的项目就是BUGS的开发以及在生物及医学方面应用。WinBUGS这个软件是我两个老板David Spiegehalter,Dave Lunn 和其它一些牛人共同开发的。我们现在正在从WinBUGS 转向openBUGS,目的是将它做成open source的软件以应用在更广的领域。 我现在正在开发BUGS中的WBDiff部分并将它应用在二型糖尿病的动态系统的数据分析中。有兴趣或有问题的同学可以和我联系： chen.wei@mrc-bsu.cam.ac.uk 还有我们这里每年会举办3-4次BUGS的培训，2天的课程，在英国或能到英国出差的同学有兴趣的话可以参加，主讲人是David speigehalter 和 Dave Lunn。 回复 · DJ on 2009/06/07 at 22:08 楼上的大牛人啊！ 回复 · maple on 2009/11/09 at 17:14 这位仁兄真牛 回复 · 海涛 on 2009/11/14 at 10:53 请问在 winbugs中我编译FRAILTY模型时提示 educational version cannot do this model 难道这个软件还有别的版本？？？我用的版本是1.4 我编译的模型是帮助文件Examples Volume 1中的最后一个文件Cox regression with random effects 回复 13. zyn on 2009/06/16 at 15:42 希望能和您联系,我这边看了gibbs sampling有不少问题，不知道您可否提供帮助？ 回复 14. phlissia on 2009/09/18 at 11:10 弱弱的问一句，执行到第四步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”， 这要怎么弄呢？ 第一次接触这个软件，还不会使，请高手指点。 WinBUGS在统计分析中的应用（第二部分） By 齐韬 @ 2008/12/18 关键词：Bayesian Analysis, SAS, WinBUGS 分类：数据分析, 统计软件, 贝叶斯统计 作者信息：Computational Mathematician in Annpro Analytic Technologies, Inc. 版权声明：本文版权归原作者所有，未经许可不得转载。原文可能随时需要修改纰漏，全文复制转载会带来不必要的误导，若您想推荐给朋友阅读，敬请以负责的态度提供原文链接；点此查看如何在学术刊物中引用本文 常规引用方式 齐韬. WinBUGS在统计分析中的应用（第二部分）. 统计之都, 2008.12. URL: http://cos.name/2008/12/statistical-analysis-and-winbugs-part-2/. BibTeX引用 @ARTICLE{, AUTHOR = {齐韬}, TITLE = {WinBUGS在统计分析中的应用（第二部分）}, JOURNAL = {统计之都}, YEAR = {2008}, month = {12}, URL = {http://cos.name/2008/12/statistical-analysis-and-winbugs-part-2/}, } 第一节 WinBUGS数据分析案例 在这一节中，我将拿一个经典的研究数据，利用WinBUGS给出简单的分析。首先介绍一下这个数据：Seeds seed O. aegyptiaco 75 seed O. aegyptiaco 73 Bean Cucumber Bean Cucumber r n r/n r n r/n r n r/n r n r/n 10 39 0.26 5 6 0.83 8 16 0.5 3 12 0.25 23 62 0.37 53 74 0.72 10 30 0.33 22 41 0.54 23 81 0.28 55 72 0.76 8 28 0.29 15 30 0.5 26 51 0.51 32 51 0.63 23 45 0.51 32 51 0.63 17 39 0.44 46 79 0.58 0 4 0 3 7 0.43 10 13 0.77 这个数据来自Crowder (1978)。之后Breslow and Clayton (1993) 作为例子，也分析过这个数据。数据反映的是某一品种的豆类种子和某一品种的黄瓜种子，分别放在21个培养皿（plates）中分别培养，在根提取物aegyptiaco 75和aegyptiaco 73的作用下的出芽率的差异。其中r是出芽的个数，n是种子的个数，而r/n是出芽率。我们用random effect logistic regression模型来进行分析（注意，在Bayesian分析中，通常是将covariates看做是服从某一个分布的随机变量，这和一般意义上的GLM, GLME, LME中对于covariates解释是不同的）： 其中是种子的类型，是根提取物的类型，是交互项, 是给定的独立的 “noninformative” 先验参数。在Bayesian分析中，通常我们会定义一个DAG图(即Directed Acyclic Graph有向无圈图) 。我们可以在WinBUGS中通过设计DAG图来定义模型。不过这一节中我们还是用WinBUGS中的BUGS语言来定义模型，如何在WinBUGS中通过设计DAG图来定义模型我将在下一节中详细介绍，但是必须要说明的是BUGS语言比DAG图灵活，不过直观性不如后者。 模型 model { for( i in 1 : N ) { r[i] ~ dbin(p[i],n[i]) b[i] ~ dnorm(0.0,tau) logit(p[i]) <- alpha0 + alpha1 * x1[i] + alpha2 * x2[i] + alpha12 * x1[i] * x2[i] + b[i] } alpha0 ~ dnorm(0.0,1.0E-6) alpha1 ~ dnorm(0.0,1.0E-6) alpha2 ~ dnorm(0.0,1.0E-6) alpha12 ~ dnorm(0.0,1.0E-6) tau ~ dgamma(0.001,0.001) sigma <- 1 / sqrt(tau) } WinBUGS doodle模型 数据 list(r = c(10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15, 32, 3), n = c(39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41, 30, 51, 7), x1 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), x2 = c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1), N = 21) 初始值 list(alpha0 = 0, alpha1 = 0, alpha2 = 0, alpha12 = 0, tau = 1) 经过10000次迭代，我们得到参数的估计如下： node mean sd MC error 2.50% median 97.50% start alpha0 -0.5546 0.1941 0.007696 -0.9353 -0.5577 -0.1597 1001 alpha1 0.08497 0.3127 0.01283 -0.5814 0.09742 0.6679 1001 alpha12 -0.8229 0.4321 0.01785 -1.697 -0.8218 0.01641 1001 alpha2 1.356 0.2743 0.01236 0.8257 1.347 1.909 1001 sigma 0.2731 0.1437 0.007956 0.04133 0.2654 0.5862 1001 第二节 结合SAS做比较分析 下面我们用同样的数据在SAS中进行分析，看看结果相差多少： 首先生成数据： data seeds; input plate r n seed extract; datalines; 1 10 39 0 0 2 23 62 0 0 3 23 81 0 0 4 26 51 0 0 5 17 39 0 0 6 5 6 0 1 7 53 74 0 1 8 55 72 0 1 9 32 51 0 1 10 46 79 0 1 11 10 13 0 1 12 8 16 1 0 13 10 30 1 0 14 8 28 1 0 15 23 45 1 0 16 0 4 1 0 17 3 12 1 1 18 22 41 1 1 19 15 30 1 1 20 32 51 1 1 21 3 7 1 1 ; run; data seeds; set seeds; interact=seed*extract; run; proc print; run; 然后调用GENMOD过程，拟合经典的logistic regression回归模型 proc genmod data=seeds; /* Basic logistic regression WHITHOUT random plate effect */ title ‘Classical logistic regression’; model r/n= seed extract interact / dist=B; run; 得到 Analysis Of Parameter Estimates Parameter DF Estimate Standard Error Wald 95% Confidence Limits Chi-Square Pr > ChiSq Intercept 1 -0.5582 0.126 -0.8052 -0.3112 19.62 <.0001 seed 1 0.1459 0.2232 -0.2915 0.5833 0.43 0.5132 extract 1 1.3182 0.1775 0.9704 1.666 55.17 <.0001 interact 1 -0.7781 0.3064 -1.3787 -0.1775 6.45 0.0111 Scale 0 1 0 1 1 调用NLMIXED过程，拟合经典的logistic regression回归模型 proc nlmixed data=seeds; /* Cassical logistic regression using NLMIXED */ title 'Classical logistic regression with NLMIXED'; parms b0=-0.55 b1=0.15 b2=1.32 b12=-0.78; logitp=b0+b1*seed+b2*extract+b12*interact; p=exp(logitp)/(1+exp(logitp)); model r ~ binomial(n,p) ; run; 得到： Parameter Estimates Parameter Estimate Standard Error DF t Value Pr > |t| Alpha Lower Upper Gradient b0 -0.5582 0.126 21 -4.43 0.0002 0.05 -0.8202 -0.2961 -0.00000229 b1 0.1459 0.2232 21 0.65 0.5203 0.05 -0.3182 0.61 -8.82E-07 b2 1.3182 0.1775 21 7.43 <.0001 0.05 0.9491 1.6872 -0.00000161 b12 -0.7781 0.3064 21 -2.54 0.0191 0.05 -1.4154 -0.1408 -6.61E-07 调用NLMIXED过程，拟合经典带随机截距的logistic regression回归模型 proc nlmixed data=seeds; /* Logistic regression + RANDOM plate INTERCEPT */ title 'Logistic regression with random plate intercept (NLMIXED)'; parms b0=-0.55 b1=0.15 b2=1.32 b12=-0.78; logitp=b0+b1*seed+b2*extract+b12*interact+alpha; p=exp(logitp)/(1+exp(logitp)); random alpha ~ normal(0,varu) subject=plate out=seedmixed; model r ~ binomial(n,p) ; run; 得到： Parameter Estimates Parameter Estimate Standard Error DF t Value Pr > |t| Alpha Lower Upper Gradient b0 -0.5484 0.1666 20 -3.29 0.0036 0.05 -0.896 -0.2009 -0.00087 b1 0.097 0.278 20 0.35 0.7308 0.05 -0.483 0.677 0.00022 b2 1.337 0.2369 20 5.64 <.0001 0.05 0.8428 1.8313 -0.00018 b12 -0.8104 0.3852 20 -2.1 0.0482 0.05 -1.6139 -0.007 0.00037 1.07 0.295 varu 0.05581 0.05 20 9 0.05 -0.0527 0.1643 0.001515 当然了winBUGS的强大之处并不在于此，而是在处理诸如GLME（有些文献称GLMM），空间数据模型等计算复杂的模型，之后还会继续讨论。 参考文献： [1] Crowder M J (1978) Beta-binomial Anova for proportions. Applied Statistics. 27, 34-37. [2] Breslow N E and Clayton D G (1993) Approximate inference in generalized linear mixed models. Journal of the American Statistical Association. 88, 9-25. 最后再送出一本 书 关于书的成语关于读书的排比句社区图书漂流公约怎么写关于读书的小报汉书pdf ，供大家研究参考 Disease Mapping with WINBUGS and MLWin（djvu格式） WinBUGS在统计分析中的应用 第二部分完 相关文章 · WinBUGS在统计分析中的应用（第三部分） (35) · 我的求学之路：经济学、软件 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 、SAS (27) · WinBUGS在统计分析中的应用（第四部分） (0) · 分类模型的性能评估——以SAS Logistic回归为例(3): Lift和Gain (6) · R与SAS之争：一个导读 (27) 13 Responses to “ WinBUGS在统计分析中的应用（第二部分） ” 1. priss111 on 2008/12/19 at 18:09 谢老大，牛就一个字！ 回复 2. 谢益辉 on 2008/12/19 at 18:22 跟我有啥关系呀…… 回复 3. 王化儒 on 2008/12/20 at 20:30 谢谢楼主的讲座和书！ to priss111：谢老大和楼主都很牛–这么说就对啦！ 回复 4. 郑冰: 统计之都《本周导读》第四辑 | 统计之都 on 2008/12/21 at 17:53 [...] WinBUGS在统计分析中的应用（第二部分） [...] 回复 5. 陈堰平 on 2008/12/24 at 11:25 期待后续文章 回复 6. xjx on 2008/12/27 at 17:43 强人那, 期待…… 回复 7. xxyishui on 2009/04/29 at 13:56 很感谢楼主的文章，让我受益很多！ 我是winbugs的新手，我在做时遇到了一个问题，向您请教一下。模型是这样的： y=e+J*P, 期中 e ~ dnorm(u0,tau0), J~dnorm(u1,tau1), P~dpois(lanta), y的数据是已知的,我想要求出参数u0,u1,tau0,tau1,and lanta. 假设先验分布: u0,u1~dnorm, tau0,tau1~dgamma, lanta~ dbeta, and y~ dnorm,这些分布是有文献支持的. 但是我发现这样做出来的结果不好：首先我先假设参数的值，产生y的数据集，然后用y来求出这些参数，对比一下，发现差距比较大。我想是否我代码编错了？ model { for(i in 1:N) { yy[i] ~ dnorm (y[i],tau) cq[i] ~ dnorm (u1,tau1) bq[i] ~ dnorm(u0,tau0) dq[i] ~ dpois(lanta) y[i] <- bq[i]+cq[i]*dq[i] } u1 ~dnorm (0.01, 0.0002) u0 ~dnorm (1, 0.0002) tau0~dgamma (4, 2) tau1 ~dgamma (4, 2) lanta~dbeta (2,1) tau<- 1/(1/tau0+1/tau1*pow(1-exp(-lanta),2)) } 回复 · icwei on 2009/06/22 at 21:05 问题不是出在你的model上，而是出在你的实验设计上。你先假设参数的值，产生y的数据集，然后用y来求出这些参数，理论上来说，只有y无限多，且产生这些y的参数遍历其所有空间时，才能返回来求出这些参数的真值。举一个简单的例子：在normal distribution P(y|theta, tau)中，我已知theta=0,tau=1,用这样的参数来产生y；假如我产生了5次data, 非常非常巧，我5次data的值都是0，因此你的数据集y={0, 0, 0, 0, 0}.那么你能用这样的数据集来估计出真正的theta和tau么？答案是否定的，你只能估计出正确的theta,而tau却非常的不一样。所以说，你产生的data要能正确地反映出你model，你才能正确的估计它的参数。这是一个parameter identifiability 的问题。希望这样的解答对你有用。 回复 · priss111 on 2009/06/22 at 23:05 icwei您好，也很想请教您一个问题： 您说的”问题不是出在你的model上，而是出在你的实验设计上” 这句话中的‘出在实验设计上’怎么理解，能不能再通俗的解释一下。 。 从您这句话中也受到启发。 因为我也遇到了一个这样的问题， 当然，不排除还有其他问题的可能性。 · icwei on 2009/06/22 at 23:43 我的意思是原作者有一个model,他为了检验这个model的正确性，便设计了一个实验，这个实验就是先假设所有参数已知，用这些参数产生data，然后再用这些data来估计这些参数，对比求得的参数和原参数的值来证明model是否正确。这样的实验没有大问题，但实施的时候要小心，你需要的检验你所产生的这些data是否能反映model的所有特性。 实际上，原作者提出的验证model的方法并不正确，即使他能够得到和原参数一样的值也并不能说明这个model对于实际的data就是正确的。他犯了一个基本的错误，既他用frequentist的方法来验证bayesian的正确性。在frenquentist的概念中，model是存在真值的，所以你可以提出一个hypothesis去检验估计值是否等于真值；而在bayesian中，我们是用credit interval 来表示所求的参数，既真值并不存在，而是与我们的prior knowledge相关的一个范围。 实际问题中，我们更多的是只得到一些数据，然后用这些数据去估计我们的参数。正确的判断参数是否能够反映系统正确性的方法是看posterior predictive distribution, 既P(y_pred|y). 如果你的posterior predictive distribution能够正确的模拟你的data，那说明这个model所估计的参数是真正能够反映你系统的特性的；反之不能。 · icwei on 2009/06/23 at 00:29 不好意思，我学统计学的时候用英文学的，所以可能用中文表达起来有些难懂。举个不完全恰当但简单的例子吧，像我上面说的，我有dataset={0， 0， 0， 0， 0}，他们是从一个有噪音的系统产生的，噪音符合normal distribution，但均值和方差均不为0。那么我用一个normal distribution的model去fit这些data,得到的结果为mean=0, variance=0.这个model可以完美的fit我当前的data.可是当我用这个model求posterior predictive distribution的时候，我所得到的是P(y_pred|y)=N(0, 0)， 其所相对应得predictive dataset为y_pred={0, 0, 0, 0, 0, …..}. 可是当我们让这个系统继续产生更多的data时，我们会发现我们model产生的值并不能很好的反映系统的特性（尽管估计的参数能够很好的fit我们原始的dataset)。 那么，我们的疑问是：我们的model错了么？？？但在这个特殊的例子里，我们可以说，model没有错，只是他估计的参数错了。那么参数不正确是由什么引起的呢？那就是我们的实验设计。这个实验设计为只得到5个数据，而这5个数据不足以反映系统的全部特性，因此才估计出错误的参数。 因此，我们说，参数是否能代表系统的特性，不但model要能反映真实系统的特性，数据也要能够反映真实系统的特性，这涉及到实验设计的问题。通常情况下，只要考虑posterior predictive distribution是否能反映系统产生的未来数据的特性就可以了。 8. priss111 on 2009/06/23 at 23:37 谢谢icwen这么认真的回复，从中又受到一些启发。 哦，我基本上明白： 对于bayesian中的参数估计，需要用后验分布来检验所估计的参数是否能真实的反应实际情形；而其他类型参数估计的模拟实验，则直接用所估计到的参数与模拟时设定的参数进行比较，越吻合说明该估计方法对所估计的参数即有效也稳健。 不知道这样理解对不对。 另外，还想请教您一个问题：不只您是不是学数学的背景？ 回复 · icwei on 2009/06/24 at 01:33 你理解的基本正确,但注意不要混淆posterior distribution(后验分布)和posterior predictive distribution(后验估计分布？不知道这个中文名是否正确)的概念。他们的区别是：假设我们的参数为a,数据为y,那么后验分布是指P(a|y),而后验估计分布是指P(y_pred|y)，它等于P(y_pred|y)=integral(P(y_pred|a)*P(a|y))da，既你对未来数据的预测是基于当前model的likelihood，以后验分布作为piror，进行加权而得到的。他是对你估计参数的所有值得一个综合性考量。 我不是纯数学背景出身，我大本学的是自控和信号处理，从博士开始专攻统计信号处理，博士毕业后开始转入统计学在生物医学方面的数学建模，现在主要兴趣在Bayesian statistics, mcmc， hierarchical graphical model, differential equation 以及 stochastic system。有兴趣的话可以和我联系，我们可以共同探讨。我的email是 chen.wei@mrc-bsu.cam.ac.uk · “统计之都”现面向广大会员征稿！在投稿之前，请仔细阅读本站的“关于页面”以及“用户如何参加‘统计之都’的各项工作”。我们现在亟需精通HTML的会员参加文章编辑工作！ · 本站接受自愿捐赠，详情参见“统计之都捐赠说明”。目前的捐赠将于3月28日用于网站升级，感谢各位的支持！ 最新评论 · zong0jie: 补充一点，系统是xp sp3。 类似测试在fedora 12... · zong0jie: 在R2.10里面貌似循环速度大大加快，只比apply慢了30... · song: 我目前是学数学与应用数学的本科生，打算往考经济统计方面研究生... · Lynn: 呵呵~一不小心看到这里，赞一个！！！羡慕john同学的复... · Vincent Tinglu Xu: 好像有问题，厂商不希望>=250g被拒绝，所以... · xuweiwei: 我也想参加啊！！！ · myli: 谢谢，弄明白了。那个工作空间，原来R for begin... · Vincent Tinglu Xu: 感兴趣江堂兄做软件的时候数学是怎么学习的。我原来也是经济学的... · Vincent Tinglu Xu: 认同这个答案。 · 邱怡轩: 1、用R软件时可能经常要读取一些外部文件，引用这些文件有两种... 随机文章 · Google Visualization API 与在线数据分析 本周一共发布了2篇日志。 使用回归分析，样本过少时不妨好先作图看看 开源的计量经济学软件gretl 二、维基 本周维基上面内容变化不大。 三、论坛 关于邀请大家成为统计之都高校联系人 · 目前会员总数：74444，主题数:9694，帖子数:67720... 2009.01.11" 统计之都《本周导读》第七辑 我想写下一段自白,这自白既是我个人的,也具有普遍意义,因为一个人经历过的事情所有的人都可以经历。 · /*跟武汉博文视点合作，召集些身边的朋友，2009应届生，计算机背景，在毕业之前，讲讲自己求学、实习、找工作等的经历与感悟，文章将由电子工业出版社结集出版，在今天秋季学期开学之前出来。我是主编，也是作者之一，刚好经历跟大伙有重叠：经济学、软件...... 2009.08.9" 我的求学之路：经济学、软件工程、SAS 在静静地看了几遍David Freedman等著的《统计学》中关于试验设计的部分后，总觉得应该写点东西发泄一下。该书自从买来就一直放在书架很久没动，也懒得动，因为翻翻前面觉得太简单。最近心情比较平静，翻了翻试验设计部分，同时思考现实中的类似例子，觉得挺有味道的。 · 试验设计包括很多内容，比如大家可能都熟知的正交表构造、区组...... 2009.05.8" 如何设计一个试验 · 首先非常感谢COS论坛提供了这样一个良好的平台，敝人心存感激之余，也打算把一些学习心得拿出来供大家分享，文中纰漏之处还请各位老师指正。下面我将以WinBUGS的统计应用为题，分几次来谈一谈WinBUGS这个软件。其中会涉及到空间数据的分析、GeoBUGS的使用、面向R及SPLUS的接口包R2WinBUGS的使用、GIS与统计分析等等衍生出的话题。如有问题，请大家留下评论，我会调整内容，择...... 2008.12.8" WinBUGS在统计分析中的应用（第一部分） 标签云 Bootstrap Confusion Matrix COS Logistic回归 P值 R R语言 SAS Sensitiveity Specificity SPSS WinBUGS 中国人民大学 中心极限定理 主站 会议 假设检验 分析数据 分类模型 博客 可视化 回归 学习经历 导读 应用统计科学研究中心 收集数据 散点图 数据挖掘 整理数据 期望 概率论 混淆矩阵 相关 空间统计 统计之都 统计功课 统计学 统计学院 维基 表述数据 计量经济学 论坛 论文摘要 贝叶斯 重抽样 WinBUGS在统计分析中的应用（第三部分） By 齐韬 @ 2009/02/11 关键词：Bayesian Analysis, GeoBUGS, WinBUGS, 空间统计, 贝叶斯统计 分类：数据分析, 统计图形, 统计软件 作者信息：Computational Mathematician in Annpro Analytic Technologies, Inc. 版权声明：本文版权归原作者所有，未经许可不得转载。原文可能随时需要修改纰漏，全文复制转载会带来不必要的误导，若您想推荐给朋友阅读，敬请以负责的态度提供原文链接；点此查看如何在学术刊物中引用本文 常规引用方式 齐韬. WinBUGS在统计分析中的应用（第三部分）. 统计之都, 2009.02. URL: http://cos.name/2009/02/