为了正常的体验网站,请在浏览器设置里面开启Javascript功能!
首页 > WinBUGS在统计分析中的应用

WinBUGS在统计分析中的应用

2019-02-24 10页 doc 40KB 27阅读

用户头像

is_348269

暂无简介

举报
WinBUGS在统计分析中的应用WinBUGS在统计分析中的应用(第一部分) 开篇词 首先非常感谢COS论坛提供了这样一个良好的平台,敝人心存感激之余,也打算把一些学习心得拿出来供大家分享,文中纰漏之处还请各位老师指正。下面我将以WinBUGS的统计应用为题,分几次来谈一谈WinBUGS这个软件。其中会涉及到空间数据的分析、GeoBUGS的使用、面向R及SPLUS的接口包R2WinBUGS的使用、GIS与统计分析等等衍生出的话题。如有问题,请大家留下评论,我会调整内容,择机给予回答。 第一节 什么是WinBUGS? WinBUGS对于研究Bayesian...
WinBUGS在统计分析中的应用
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目前是一款免费的软件,去下载就好了。不过要用高级功能(如GeoBUGS)的话,还是去注册一下好了,挺方便的。系统会立即把注册码发到你邮箱(真是好人啊)。不过只可以用一个月。这倒无妨,到时在注册一下就好了。 第四节 初试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 附第二步中的代码如下: #MODELmodel{ 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}#DATAlist(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 VALUESlist(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数据分析案例 在这一节中,我将拿一个经典的研究数据,利用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) 继续阅读
/
本文档为【WinBUGS在统计分析中的应用】,请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑, 图片更改请在作品中右键图片并更换,文字修改请直接点击文字进行修改,也可以新增和删除文档中的内容。
[版权声明] 本站所有资料为用户分享产生,若发现您的权利被侵害,请联系客服邮件isharekefu@iask.cn,我们尽快处理。 本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用。 网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。

历史搜索

    清空历史搜索