 上传我的文档
 下载
 收藏
该文档贡献者很忙，什么也没留下。
 下载此文档
8 MC与MCMC方法
下载积分：100
内容提示：8 MC与MCMC方法
文档格式：PDF|
浏览次数：245|
上传日期： 10:44:42|
文档星级：
全文阅读已结束，如果下载本文需要使用
 100 积分
下载此文档
该用户还上传了这些文档
8 MC与MCMC方法
官方公共微信【图文】MCMC方法介绍_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
MCMC方法介绍
大小：3.13MB
登录百度文库，专享文档复制特权，财富值每天免费拿！
你可能喜欢苹果/安卓/wp
积分 36, 距离下一级还需 9 积分
道具: 彩虹炫, 涂鸦板, 雷达卡, 热点灯, 金钱卡下一级可获得
道具: 显身卡
购买后可立即获得
权限: 隐身
道具: 金钱卡, 彩虹炫, 雷达卡, 热点灯, 涂鸦板
开心签到天数: 3 天连续签到: 1 天[LV.2]偶尔看看I
我要做用mcmc方法求正态分布变点，结果程序老是运行不了，有没有人能给看看，下面是程序：
y1&-rnorm(8,0.8,2)
y2&-rnorm(6,1.2,2)
y3&-rnorm(6,0.4,2)
y&-c(y1,y2,y3)
mhsampler=function(numit,dat=y)
mchain=matrix(NA,6,numit)
mchain[,1]=c(0.2,0.4,0.6,1,5,10)
like1&-function(a,b,c)
sum(((y-c)^2)[a:b])
for(i in 2:numit)
y1hat&-sum(y[1:ck1])/ck1
y2hat&-sum(y[(ck1+1):ck2])/(ck2-ck1)
y3hat&-sum(y[(ck2+1):n])/(n-ck2)
ctheta1=mchain[1,i-1]
ctheta2=mchain[2,i-1]
ctheta3=mchain[3,i-1]
csigma=mchain[4,i-1]
ck1=mchain[5,i-1]
ck2=mchain[6,i-1]
(ctheta1=rnorm(1,y1hat,csigma/ck1))
(ctheta2=rnorm(1,y2hat,csigma/(ck2-ck1)))
(ctheta3=rnorm(1,y3hat,csigma/(n-ck2)))
csigma=1/rgamma(1,n-1,(like1(1,ck1,y1hat)+like1(ck1+1,ck2,y2hat)+like1(ck2+1,n,y3hat))/2)
(pk1=sample(x=seq(2,ck2-1),1))
(MHratio1&-exp((like1(1,pk1,ctheta1)+like1(pk1+1,ck2,ctheta2))/(-2*csigma))/exp((like1(1,ck1,ctheta1)+like1(ck1+1,ck2,ctheta2))/(-2*csigma))
(alpha1=min(1,MHratio1))
(ck1&-ifelse((runif(1))&alpha1,pk1,ck1))
(pk2=sample(x=seq(ck1+1,n-1),1))
(MHratio2&-exp((like1(ck1+1,pk2,ctheta2)+like1(pk2+1,n,ctheta3))/(-2*csigma))/exp((like1(ck1+1,ck2,ctheta2)+like1(ck2+1,n,ctheta3))/(-2*csigma))
(alpha2=min(1,MHratio2))
(ck2&-ifelse((runif(1))&alpha2,pk2,ck2))
mchain[,i]=c(ctheta1,ctheta2,ctheta3,csigma,ck1,ck2)
return(mchain)
apply(mhsampler(1000,dat=y),1,mean)
支持楼主：、
购买后，论坛将把您花费的资金全部奖励给楼主，以表示您对TA发好贴的支持
载入中......
不要沉啊！！！！！！！自己顶一下
y1&-rnorm(8,0.8,2)
y2&-rnorm(6,1.2,2)
y3&-rnorm(6,0.4,2)
y&-c(y1,y2,y3)
n=20
ck1=5
ck2=10
mhsampler=function(numit,dat=y)
{
mchain=matrix(NA,6,numit)
mchain[,1]=c(0.2,0.4,0.6,1,5,10)
like1&-function(a,b,c)
{
sum(((y-c)^2)[a:b])
}
for(i in 2:numit)
{
y1hat&-sum(y[1:ck1])/ck1
y2hat&-sum(y[(ck1+1):ck2])/(ck2-ck1)
y3hat&-sum(y[(ck2+1):n])/(n-ck2)
ctheta1=mchain[1,i-1]
ctheta2=mchain[2,i-1]
ctheta3=mchain[3,i-1]
csigma=mchain[4,i-1]
ck1=mchain[5,i-1]
ck2=mchain[6,i-1]
(ctheta1=rnorm(1,y1hat,csigma/ck1))
(ctheta2=rnorm(1,y2hat,csigma/(ck2-ck1)))
(ctheta3=rnorm(1,y3hat,csigma/(n-ck2)))
csigma=1/rgamma(1,n-1,(like1(1,ck1,y1hat)+like1(ck1+1,ck2,y2hat)+like1(ck2+1,n,y3hat))/2)
(pk1=sample(x=seq(2,ck2-1),1))
(MHratio1&-exp((like1(1,pk1,ctheta1)+like1(pk1+1,ck2,ctheta2))/(-2*csigma))/exp((like1(1,ck1,ctheta1)+like1(ck1+1,ck2,ctheta2))/(-2*csigma)))
(alpha1=min(1,MHratio1))
(ck1&-ifelse((runif(1))&alpha1,pk1,ck1))
(pk2=sample(x=seq(ck1+1,n-1),1))
(MHratio2&-exp((like1(ck1+1,pk2,ctheta2)+like1(pk2+1,n,ctheta3))/(-2*csigma))/exp((like1(ck1+1,ck2,ctheta2)+like1(ck2+1,n,ctheta3))/(-2*csigma)))
(alpha2=min(1,MHratio2))
(ck2&-ifelse((runif(1))&alpha2,pk2,ck2))
mchain[,i]=c(ctheta1,ctheta2,ctheta3,csigma,ck1,ck2)
}
return(mchain)
}
apply(mhsampler(1000,dat=y),1,mean)复制代码
热心帮助其他会员
总评分:&论坛币 + 10&
热心指数 + 1&
MHratio1，MHratio2最后少加了一个括号，你能给我发一份推到过程吗？，谢谢
zhou1_20 发表于
MHratio1，MHratio2最后少加了一个括号，你能给我发一份推到过程吗？，谢谢谢谢！！我最近没逛论坛，在准备开题，我给你发送了
Metropolis-Hasting算法有R的包的，不必自己写～
foozhencheng 发表于
Metropolis-Hasting算法有R的包的，不必自己写～你好，能讲讲用什么包么
Zjn 发表于
你好，能讲讲用什么包么可以查一下MHadaptive
好东西，谢谢分享
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
&nbsp&nbsp|
如有投资本站或合作意向，请联系（010-）；
邮箱：service@pinggu.org
投诉或不良信息处理：（010-）
论坛法律顾问：王进律师{"debug":false,"apiRoot":"","paySDK":"/api/js","wechatConfigAPI":"/api/wechat/jssdkconfig","name":"production","instance":"column","tokens":{"X-XSRF-TOKEN":null,"X-UDID":null,"Authorization":"oauth c3cef7c66aa9e6a1e3160e20"}}
{"database":{"Post":{"":{"contributes":[{"sourceColumn":{"lastUpdated":,"description":"专栏作者在工作中时常需要用到数据分析、统计学/统计推断，但并不是统计学科班出身，对于统计学的系统性学习仅仅止于本科课程概率论与数理统计。在工作中，作者发现需要更为深入地了解统计学的现代方法。并且随着大数据、机器学习等热门关键词时常见诸各类媒体，作者认为有必要赶一赶这个时髦。遂申请开设了此专栏，以分享作者作为一个统计学初级用户的学习历程。\n专栏头像是摩纳哥的Monte-Carlo赌场\nBy Fruitpunchline - Own work, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=","permission":"COLUMN_PUBLIC","memberId":,"contributePermission":"COLUMN_PRIVATE","translatedCommentPermission":"all","canManage":true,"intro":"分享科学数据分析的知识、经验与心得","urlToken":"data-analysis","id":12405,"imagePath":"v2-e53fdbd4cec6a21ded9551d8cfa10c84.jpg","slug":"data-analysis","applyReason":"0","name":"数据的背后","title":"数据的背后","url":"/data-analysis","commentPermission":"COLUMN_ALL_CAN_COMMENT","canPost":true,"created":,"state":"COLUMN_NORMAL","followers":106,"avatar":{"id":"v2-e53fdbd4cec6a21ded9551d8cfa10c84","template":"/{id}_{size}.jpg"},"activateAuthorRequested":false,"following":false,"imageUrl":"/v2-e53fdbd4cec6a21ded9551d8cfa10c84_l.jpg","articlesCount":5},"state":"accepted","targetPost":{"titleImage":"/v2-5c7cbb4b9fa300ac1de0f1dc3568fa3c_r.png","lastUpdated":,"imagePath":"v2-5c7cbb4b9fa300ac1de0f1dc3568fa3c.png","permission":"ARTICLE_PUBLIC","topics":[,27412],"summary":"终于，又有空添加一篇关于统计参数推断的文章了。 关于统计推断的数学，作者暂时没有挖掘到足够的深度，所以就不在此涉及更多数学的本质了。我们来谈一谈线性logit回归模型的参数推断吧。 所谓回归，如果以不太深刻的视角去审视，大概和曲线拟合有很大程度…","copyPermission":"ARTICLE_COPYABLE","translatedCommentPermission":"all","likes":0,"origAuthorId":0,"publishedTime":"T23:52:45+08:00","sourceUrl":"","urlToken":,"id":2504490,"withContent":false,"slug":,"bigTitleImage":false,"title":"基于MCMC方法的Logit回归模型参数推断","url":"/p/","commentPermission":"ARTICLE_ALL_CAN_COMMENT","snapshotUrl":"","created":,"comments":0,"columnId":12405,"content":"","parentId":0,"state":"ARTICLE_PUBLISHED","imageUrl":"/v2-5c7cbb4b9fa300ac1de0f1dc3568fa3c_r.png","author":{"bio":"不以物喜，不以己悲","isFollowing":false,"hash":"9a150d1993f1daba5c82","uid":515900,"isOrg":false,"slug":"astrojhgu","isFollowed":false,"description":"","name":"astrojhgu","profileUrl":"/people/astrojhgu","avatar":{"id":"da8e974dc","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},"memberId":,"excerptTitle":"","voteType":"ARTICLE_VOTE_CLEAR"},"id":579936}],"title":"基于MCMC方法的Logit回归模型参数推断","author":"astrojhgu","content":"终于，又有空添加一篇关于统计参数推断的文章了。 关于统计推断的数学，作者暂时没有挖掘到足够的深度，所以就不在此涉及更多数学的本质了。我们来谈一谈线性logit回归模型的参数推断吧。 所谓回归，如果以不太深刻的视角去审视，大概和曲线拟合有很大程度上的重合。但是回归在某种程度上更强调对未知数据的推断。很多人接触到的第一个回归模型大概是线性回归，也就是基于一系列的观测数据点对拟合出一条曲线，并用这个模型去推断对应于某个的某个未被观察到的的数值。在很多函数型计算器上这都是很基本的功能，喜欢玩计算器的同学一定不会错过。我们今天要讨论的模型更复杂一些。回忆普通线性回归，自变量和因变量都是连续型的变量，而在logit回归中因变量是离散的，确切地说，通常是一个二值变量（如果有同学发现有多值变量的话告诉我，反正我见到的例题里面都是二值的，anyway，这只是一个名称的问题），不失一般性，可以将其值映射到集合{0,1} 上。自变量通常是连续的，也可以有离散的。我们用数学的语言把问题表达出来：线性logit回归就是估计如下模型：也就是线性关系并非直接体现在自变量和因变量之间的联系，而是体现在自变量和因变量分布的参数的关系上。在这里，因和自变量之间是线性的。在这里，的表达式是,我们把它plot出来是这个样子的：而是logit函数的反函数，其表达式为，本文题图就是该函数的图形。这个函数在这里的作用就是把映射到这个区间，而这个区间正好是概率p的取值范围。以上作者罗里八嗦 一大堆，只是把logit回归的模型给大致交代了一下，那么如何用MCMC方法对其参数进行推断呢？这就需要构建一个参数的后验分布,其中。上面连乘号里面的因子是每一个的分布函数，实际上就是的，是参数的先验分布。写出这样一个后验分布之后，我们就能用某种采样方法对其进行采样，然后得到一个的样本，基于此可以抽取出参数的统计特征(比如期望、各阶矩量)。我们拿一个具体的例子看一下吧，这个例子来自于，其中数据从下载。数据格式是一个 CSV文件，它的头几行是这个样子的：##
gpa rank\n## 1
0 380 3.61
1 660 3.67
1 800 4.00
1 640 3.19
0 520 2.93
1 760 3.00
2 \n看样子是对一个考生录取情况的统计，一共有四列数据，分别是录取与否(录取的话取值1，否则取值0)、GRE分数、GPA，和考生来自的那个单位的声誉水平(从1到4，声誉依次递减)。因为我比较懒，所以我们跳过本应该一开头做的数据清洗等工作，直接跳到建模这一步，按照上面的公式，写出后验分布的模型：先验分布直接取均匀分布，这样就不出现在公式里面了。然后我们用里面提到的emcee方法对这样一个后验分布进行采样50000次，前1000次作为预烧不输出，得到一个样本，在进行进一步的分析之前我们先看看结果有没有收敛，以样本序号为横轴，参数为纵轴把数据plot出来： 长得很像白噪声，看上去已经收敛得非常好了。当然严格地说应当做定量的收敛分析，比如说计算自相关函数什么的，不过我比较懒，这里就忽略了。 长得很像白噪声，看上去已经收敛得非常好了。当然严格地说应当做定量的收敛分析，比如说计算自相关函数什么的，不过我比较懒，这里就忽略了。然后我们把pair plot给画出来，是这个样子的：对角线上的图是每一个变量的直方图，其他的图是四个变量两两之间的联合分布等高线图，从里到外分别代表的界限。好了，从图上我们能看出啥呢？大体上，也就是说GRE高的人、GPA高的人、学校好的人更容易被录取。而且注意到的概率远远小于千分之三，说明rank越小(学校越好)的人越容易被录取这个结论的置信度极高，果然出身很重要啊。而和都有不算非常小的概率可以小于0，说明尽管GRE和GPA都很重要，但和学校rank比起来都不是最重要的。和之间有轻微的反相关关系，也许可以解读为GPA高的人考GRE都不怎么上心？不好说，数据本身是一回事情，数据的解读又是另外一回事情，不能听风就是雨，将来报道出了偏差等于作者也有责任的。——果然还是分析错了，见下面勘误是截距，没有很明显的意义。好了，用上面这个例子解释了一下基于MCMC方法的logit回归参数推断问题。需要注意的是，上面做的这些其实忽略了很多步骤，比如说，至少应该把原始数据plot出来，直观感受一下用线性的模型是否合适，有没有outlier什么的，真的进行数据分析的时候这些步骤都是不应当忽略的，否则很容易落入对数据过度诠释的陷阱之中，有意无意地成为数据也会骗人的受害者。就到这里吧。ps，有人关心pair-plot是怎么画出来的，简单说一下：工具是python的matplotlib，据我所知seaborn也能画个大概，但是不能满足这里的要求，无法定量画出n-sigma的界线。算法的话，直方图不说了，直接调用histogram函数就行了。两两变量联合分布的contour是这么画的 ，先用histogram2d得到二维直方图的数值（一个二维数组），然后对数值进行排序（从大到小），然后从大到小累积得到累积分布，最后数出来99.7%,95%和68%的quantile，用matplotlib的contour函数在二维直方图的图上用quantile作为level画出直方图。h,xedges,yedges=np.histogram2d(x,y,bins=[...])\n\nlevels=contour_levels(h,[.997,.95,.68])\n\nax.contour(xc,yc,h,levels,linewidths=1,colors='black',linestyles=['dotted','dashed','solid'])\n其中contour_level的定义是这样的：def contour_levels(h,levels=[.997,.95,.68]):\n
d=h.flatten().tolist()\n
s=h.sum()\n
d.sort()\n
d=d[::-1]\n
cum_d=sp.cumsum(d)\n
result=[]\n
for i in levels:\n
n=next(((n,x) for (n,x) in enumerate(cum_d) if x&=i*s),None)[0]\n
result.append(d[n])\n
result.sort()\n
return result\n
\n勘误:上面说a和 b存在负相关，并试图解读为GPA高的人考GRE都不怎么上心，这是错误的，a和b只是两个系数，它们的负相关并不能导出GPA和GRE成绩这两个观测量之间的负相关。","updated":"T15:52:45.000Z","canComment":false,"commentPermission":"anyone","commentCount":5,"collapsedCount":0,"likeCount":7,"state":"published","isLiked":false,"slug":"","lastestTipjarors":[],"isTitleImageFullScreen":false,"rating":"none","titleImage":"/v2-5c7cbb4b9fa300ac1de0f1dc3568fa3c_r.png","links":{"comments":"/api/posts//comments"},"reviewers":[],"topics":[{"url":"/topic/","id":"","name":"统计推断"},{"url":"/topic/","id":"","name":"数理统计学"},{"url":"/topic/","id":"","name":"贝叶斯统计"}],"adminClosedComment":false,"titleImageSize":{"width":1280,"height":853},"href":"/api/posts/","excerptTitle":"","column":{"slug":"data-analysis","name":"数据的背后"},"tipjarState":"activated","tipjarTagLine":"真诚赞赏，手留余香","sourceUrl":"","pageCommentsCount":5,"tipjarorCount":0,"annotationAction":[],"hasPublishingDraft":false,"snapshotUrl":"","publishedTime":"T23:52:45+08:00","url":"/p/","lastestLikers":[{"bio":null,"isFollowing":false,"hash":"cb224e738c11bc89c543ed","uid":32,"isOrg":false,"slug":"hu-bin-36-78","isFollowed":false,"description":"","name":"胡斌","profileUrl":"/people/hu-bin-36-78","avatar":{"id":"da8e974dc","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},{"bio":"PRIS小蛀虫","isFollowing":false,"hash":"08e4c0b77e3e265a11e1","uid":84,"isOrg":false,"slug":"wai7niu8","isFollowed":false,"description":"","name":"邮递员小王","profileUrl":"/people/wai7niu8","avatar":{"id":"b7d61e6ec3c73b4a8ad375","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},{"bio":"警校生","isFollowing":false,"hash":"2f6d8f67dfdbdb096a2189fda9f492cd","uid":360800,"isOrg":false,"slug":"xiu-mu-44-33","isFollowed":false,"description":"","name":"朽木","profileUrl":"/people/xiu-mu-44-33","avatar":{"id":"v2-253e7b82e9c91eabdf0b2eb","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},{"bio":"经济学爱好者/滑雪爱好者/羽毛球爱好者","isFollowing":false,"hash":"302e6da83a29c394a55cb68d5772fce7","uid":951100,"isOrg":false,"slug":"li-hu-shi-kou","isFollowed":false,"description":"","name":"严嵩","profileUrl":"/people/li-hu-shi-kou","avatar":{"id":"f41de21d3fd6eacaf047e","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},{"bio":"All models are wrong, some are deadly","isFollowing":false,"hash":"780d99499aded1ef8632","uid":84,"isOrg":false,"slug":"ding-yi-88-84","isFollowed":false,"description":"做度量的","name":"丁一","profileUrl":"/people/ding-yi-88-84","avatar":{"id":"480cf022b","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false}],"summary":"终于，又有空添加一篇关于统计参数推断的文章了。 关于统计推断的数学，作者暂时没有挖掘到足够的深度，所以就不在此涉及更多数学的本质了。我们来谈一谈线性logit回归模型的参数推断吧。 所谓回归，如果以不太深刻的视角去审视，大概和曲线拟合有很大程度…","reviewingCommentsCount":0,"meta":{"previous":{"isTitleImageFullScreen":false,"rating":"none","titleImage":"/v2-d1e4afb838bfa0cbac698bba_r.jpg","links":{"comments":"/api/posts//comments"},"topics":[{"url":"/topic/","id":"","name":"概率论"},{"url":"/topic/","id":"","name":"概率论与数理统计"},{"url":"/topic/","id":"","name":"统计推断"}],"adminClosedComment":false,"href":"/api/posts/","excerptTitle":"","author":{"bio":"不以物喜，不以己悲","isFollowing":false,"hash":"9a150d1993f1daba5c82","uid":515900,"isOrg":false,"slug":"astrojhgu","isFollowed":false,"description":"","name":"astrojhgu","profileUrl":"/people/astrojhgu","avatar":{"id":"da8e974dc","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},"column":{"slug":"data-analysis","name":"数据的背后"},"content":"在和文章中，我们主要说了基于贝叶斯图和吉布斯采样方法的统计推断。然而吉布斯采样方法在下却几乎无法工作，其中最具有代表性的就是参数空间中存在多个相互独立的“概率孤岛”，孤岛之间跳转的概率非常小。这使得吉布斯采样器并非一种理想的不需调参的通用采样器。幸好除了吉布斯采样器之外，还存在其他选择。根据我的调研，有一种称为的算法是其中比较好用的一种。该方法是一种基于热力学系综思想的采样器。在每一步的采样中，该算法都根据上一轮采样得到的采样得到一组新的分布系综。相对于基于的吉布斯采样方法，emcee具有更高的计算效率，而相对于算法的吉布斯采样方法，emcee则具有不需要提供建议分布、不需要对参数进行精调的优点。我认为它在很多情况下都要优于吉布斯采样算法。当然emcee算法仍然不能很好地处理“概率孤岛” 的情况，只需要基于emcee方法做一些改进，这一问题就能得到很好的解决。这一改进就是引入所谓的“”。和emcee不同，parallel tempering是另一个范畴上层次更高的采样方法，它的底层可以是emcee，也可以是其它的采样方法，例如吉布斯采样方法。它引入了一个被称为“温度”的参数，这个参数对于系统的影响是这样的：它的数值越大，被采样的分布就越接近均匀分布()，从而概率孤岛之间的“海峡”就被填上了。这一方法同时对一系列温度参数调控下的概率分布进行采样，并且在某些时候根据一定判据，交换相邻的温度参数对应的马尔科夫链之间的参数的参数。最低的温度为1，对应于原始的分布。这一方法有两个需要调整的参数，就是温度列表和交换的频率。我们现在用基于emcee的parallel tempering 方法对一个双高斯混合模型进行采样。,,其中。这是一个明显的性质的分布，预计的数据点将主要集中在和 这两个区间内，然而两区间并不重叠，即, ,,这四个描述系统状态在那两个孤岛之间转移的概率非常接近于零。如果强行用吉布斯采样器对这个概率分布进行采样，在有限的时间内是无法达到平衡的。最后的结果就表现再数据点只出现在一个孤岛内，具体哪一个视初始条件而定。改为采用基于ensemble sampling 的parallel tempering方法对它进行采样，得到的不同温度的马尔科夫链数据的累积分布如下：不同颜色的线条代表不同的温度下的马尔科夫链的累积分布，黑色的线条就是我们想要采样的分布。可以发现温度越高的数据分布越趋向于均匀分布。这是累积分布函数，也就是概率密度的积分。直方图在此就不画了，从累积分布函数很容易想象出概率密度分布函数。从结果来看，显然是符合预期的。 题图：By Diacritica - Own work, CC BY-SA 3.0, ","state":"published","sourceUrl":"","pageCommentsCount":0,"canComment":false,"snapshotUrl":"","slug":,"publishedTime":"T23:14:09+08:00","url":"/p/","title":"吉布斯之外的其它选择","summary":"在和文章中，我们主要说了基于贝叶斯图和吉布斯采样方法的统计推断。然而吉布斯采样方法在下却几乎无法工作，其中最具有代表性的就是参数空间中存在多个相互独立的“概率孤岛”，孤岛之间跳转的概率非常小。这使得吉布斯采样器并非一种…","reviewingCommentsCount":0,"meta":{"previous":null,"next":null},"commentPermission":"anyone","commentsCount":0,"likesCount":1},"next":{"isTitleImageFullScreen":false,"rating":"none","titleImage":"/v2-ec4a09e45a799_r.jpg","links":{"comments":"/api/posts//comments"},"topics":[{"url":"/topic/","id":"","name":"贝叶斯网络"},{"url":"/topic/","id":"","name":"概率论与数理统计"},{"url":"/topic/","id":"","name":"机器学习"}],"adminClosedComment":false,"href":"/api/posts/","excerptTitle":"","author":{"bio":"不以物喜，不以己悲","isFollowing":false,"hash":"9a150d1993f1daba5c82","uid":515900,"isOrg":false,"slug":"astrojhgu","isFollowed":false,"description":"","name":"astrojhgu","profileUrl":"/people/astrojhgu","avatar":{"id":"da8e974dc","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false},"column":{"slug":"data-analysis","name":"数据的背后"},"content":"纸上得来终觉浅，绝知此事须躬行
—— 陆游《冬夜读书示子聿》引言 在此专栏的第一和第二篇文章中，作者大致阐述了基于贝叶斯网络和吉布斯采样算法从观测数据推断参数统计特征的概念和原理。然而，脱离具体实践操作的纯原理性论述对于理解这一统计工具的使用是不够。在本篇文章中，作者决定用一个具体的、实际的例子，结合一个统计推断工具——JAGS来阐述。JAGS软件 JAGS，即Just Another Gibbs Sampler，从名字看它是一个吉布斯采样工具。然而作者认为这个名字并不能很好地反映它的实质。实际上它应当被称为Just Another Bayesian Network Modeling Tool。因为从根本上说，它是一款贝叶斯网络建模工具，Gibbs采样算法只是对它所建立模型所代表的那个后验分布进行采样的一个算法。回忆本专栏的第三篇文章，我们完全可以用贝叶斯网络方法完成建模，但是不用Gibbs方法进行采样，而用ensemble method进行采样，最终可以得到同样的结果。以上是对这个软件名字的一些吐槽。JAGS软件是一款命令行界面的软件，运行在linux下（有没有windows我没有调研过）。JAGS读入一个建模文件 .bug（这个bug其实来自于JAGS软件的前身——Bayesian inference Using Gibbs Sampling）、一个可选的数据文件，并根据用户的命令行交互操作产生一系列采样结果。问题描述在本节中，作者将给出一个实际的统计推断问题的描述，在下一节中，给出JAGS的建模代码。这个问题的例子来自与作者多次引用的这本书的第8.2节。这个问题是这样的：某人开展了一项核物理实验，研究某种反应发生的概率和注入粒子能量之间的关系（如果这么表述还是太抽象，那么可以理解为用不同能量的粒子去轰击靶子，会得到两种结果，一种是打偏，另一种是打中，现在想知道打中的概率和该粒子的能量之间的关系）。实验结束后研究人员得到了一个数据表格，表格有三列，若干行，每一行代表粒子的能量E等于某个的情况下的数据。三列分别是1. 粒子能量E2. 打中的粒子数nrec3. 注入的粒子数 ninj 数据表格如下：从表格中可以看出，nrec始终小于或等于ninj，而且nrec/ninj的比值随能量的上升而上升，且在高能端出现了饱和。这是对数据的第一直观印象。下一步我们将本次实验中测得的反应比例nrec/ninj （在上一版本中，我用了反应概率这个词，这是不对的，概率是隐含的，不可测的，比例是可测的） 随E的关系plot出来看看，以便决定应该采用哪种模型来表示这一关系。对数据的直观感受 我们用jupyter notebook来画这个图，画图的命令是 结果是可以发现趋势和上面的直观分析是一致的。按照模型拟合的思路，我们现在只需要构建一个具有解析表达式的模型Eff=f(E)，并做最小二乘拟合或者假设一个泊松误差，进行chi^2拟合就得到结果了。但是根据我们之前的阐述可知，这样做得不到参数误差的信息。所以接下来我们要在贝叶斯统计的框架下对模型参数进行估计。参数估计的第一步就是建模，我们要思考，什么样的模型可以表示Eff~E的关系呢？这当然是一个以核物理为背景的实验数据分析问题， 这本书中的第8.2节已经给出了答案：用Phi函数来表示。Phi函数就是高斯概率密度函数的累积分布。它是一个特殊函数erf的变形，具体说 。把它plot出来是这个样子的：观察x&0的部分，是不是和上面的数据点的plot很像呢？所以我们可以假定这里的Eff就是能量E情况下的反应概率，或者说是效率，其中A、B、、是四个待估计参数。 这个公式代表了对Phi函数的平移和缩放。需要注意的是，这一方案（即用Phi函数近似反应概率和能量的关系）并不是唯一合理的方案，如果能基于物理原理得出一个公式，那会更好。接下来的的问题是如何在贝叶斯的框架下把这个模型给表达出来并和测量数据E、ninj、nrec给联系起来。基于物理图景建立模型 相对于最近热门的神经网络方法，贝叶斯网络的优点就是其可解释性，也就是我们能很好地基于内在机制进行建模。 我们先从物理的角度分析一下nrec、ninj和Eff的关系，然后给出JAGS模型定义语句。对于一个给定的能量E，我们都有一个内在的无法被观测到的隐含变量Eff，也就是一个概率，这个概率决定了我们发射ninj个粒子，能有多少个粒子发生反应，即nrec。nrec的本质是一个随机变量，我们能观察到的只是它的一个realization。那么我们要问nrec服从什么分布？显然，在这一问题中是最恰当的。关于什么是二项分布以及为什么我们在此选用这一分布： 二项分布是指进行N次独立的试验，每次试验的结果只有两个取值0和1，结果为1的概率是p，为0的概率是(1-p)，那么N次试验后结果为1的次数的分布情况，显然伯努利分布的模型和这一问题是非常匹配的。所以我们在jags的模型定义文件中写下：for (i in 1:length(nrec)) {\n
eff[i] &- A + (B-A)*phi((E[i]-mu)/sigma)
nrec[i] ~ dbin(eff[i],ninj[i])\n}\n 这四行代码很容易理解：第1行是一个循环，代表下标 i “跑遍” [1, 2, ..., length(nrec)] 。第2行是eff的表达式，就是上面公式的内容。这里用到了JAGS的内置函数phi。第3行是代表nrec分布的表达式，用到了JAGS的内置分布函数dbin即二项分布的概率函数，它有两个参数，第一个是每次实验得到结果1的概率p（也就是这里发生反应的概率），第二个参数是实验次数（也就是这里的注入粒子数ninj）。注意到在JAGS中用符号“ &- ” 代表生成确定性节点（近似理解为赋值），用符号 “ ～ ”代表某个变量遵从某个概率分布，也就是生成随机节点（stochastic node, 关于什么是节点、什么是随机节点、什么是确定性节点，参看本专栏的和篇文章）。需要注意的是，随机节点又可分为observed和unobserved，也就是具有观测值的和不具有观测值的。什么样的节点是observed的呢？见下面关于数据文件的描述。注意到，nrec是一个具有观测值的量，所以第3行代表的是在nrec具有当前观测值条件下，参数A、B、、取某一组具体值情况下的概率，换句话说，就是这四个参数的似然函数。有了似然函数，我们还需要有先验分布，根据经验，我们将这四个参数的先验分布选取为如下形式：A~dunif(0,1)\nB~dunif(0,1)\nmu~dunif(0,100)\nsigma~dunif(0,100)\n 再一次，我们用到了符号“ ~ ” 来生成随机节点，这是四个均匀分布，用到了JAGS的内置分布dunif。这个分布有两个参数，也就是均匀分布的上下限。至此，我们已经把模型文件给写好了。把上面的结果综合一下，model.bug文件的全部内容是：model {\nfor (i in 1:length(nrec)) {\n
nrec[i] ~ dbin(eff[i],ninj[i])\n
eff[i] &- A + (B-A)*phi((E[i]-mu)/sigma)\n}\nA~dunif(0,1)\nB~dunif(0,1)\nmu~dunif(0,100)\nsigma~dunif(0,100)\n}\n 然后我们还要告诉JAGS观测量的具体取值，所以我们就写一个data.R文件E &-c(2, 6, 10, 14, 18, 22, 26, 30, 34, 38, 42, 46, 50, 54, 58, 62, 66, 70, 74, 78, 82, 86, 90, 94, 98)\n\nnrec &-c(23, 71, 115, 159, 200, 221, 291, 244, 277, 221, 210, 182, 136, 119, 79, 81, 61, 44, 41, 32, 32, 31, 22, 18, 11)\n\nninj &-c(96, 239, 295, 327, 345, 316, 349, 281, 298, 235, 217, 185, 140, 121, 79, 81, 61, 45, 41, 32, 32, 31, 22, 18, 11)\n 熟悉R的同学可能看出来了，这就是R语言的向量声明语句。所以上文中提到的observed和unobserved stochastic node的区分就有了答案——凡是在data文件中指定了取值的节点就是observed的，否则就是unobserved的。JAGS只对unobserved stochastic node进行采样，observed stochastic node由于取值确定，所以不采样。而deterministic node的取值是从父节点计算出来的，所以也是有确定取值的，无须采样。用JAGS运行模型
然后我们就要把结果载入到JAGS中进行采样，也就是运行模型。为了使我们对JAGS的操作持久化，我们写一个.cmd文件：model in \"model.bug\"\ndata in \"data.R\"\ncompile, nchains(1)\ninitialize\nupdate 1000\nmonitor A\nmonitor B\nmonitor mu\nmonitor sigma\nupdate 30000\ncoda *\n 我们从上到下依次解释：model in ... 表示模型定义文件data in ... 表示数据定义文件compile, nchans(1)代表让JAGS将我们定义在model.bug中模型字符串编译为JAGS的内置数据结构，并检查语法。 nchains(1)代表我们只运行一条马尔科夫链initialize显然代表初始化update 1000代表让JAGS先采样1000次。因为MCMC方法存在一个收敛于平稳分布的问题，最初采样的那部分结果尚未收敛，不应当被采用，所以我们先跑1000次，而不记录数据。跑完1000次之后，我们接下来跑的结果都要被记录下来了。所以下面的四行 monitor ...就代表希望JAGS记录哪些参数的采样结果。update 30000代表采样30000次，得到一个容量为30000的参数样本。coda * 代表将结果存为文件。 然后在terminal中执行jags &run.cmd\n就得到了结果——两个文件，分别为CODAchain1.txt和CODAindex.txt，前者是一个文本形式的数据文件，一个数据占一行，CODAindex.txt的内容如下：A 1 30000B mu sigma
含义很明确，代表CODAChain1.txt中哪些行对应哪个参数。结果解读我们利用程序将结果给plot出来：从图中我们可以读出没一个参数的分布直放图、两两参数的互相关分布。为了得到定量结果，我们统计每一个参数的均值和标准差，结果如下：所以最终我们写出Eff~E的表达式：。参考文献： 作者： & Brian Weaver 题图取自：","state":"published","sourceUrl":"","pageCommentsCount":0,"canComment":false,"snapshotUrl":"","slug":,"publishedTime":"T15:54:54+08:00","url":"/p/","title":"一种进行贝叶斯统计推断的工具——JAGS","summary":"纸上得来终觉浅，绝知此事须躬行 —— 陆游《冬夜读书示子聿》引言 在此专栏的第一和第二篇文章中，作者大致阐述了基于贝叶斯网络和吉布斯采样算法从观测数据推断参数统计特征的概念和原理。然而，脱离具体实践操作的纯原理性论述对于理解这一统计工具的使…","reviewingCommentsCount":0,"meta":{"previous":null,"next":null},"commentPermission":"anyone","commentsCount":0,"likesCount":1}},"annotationDetail":null,"commentsCount":5,"likesCount":7,"FULLINFO":true}},"User":{"astrojhgu":{"isFollowed":false,"name":"astrojhgu","headline":"","avatarUrl":"/da8e974dc_s.jpg","isFollowing":false,"type":"people","slug":"astrojhgu","bio":"不以物喜，不以己悲","hash":"9a150d1993f1daba5c82","uid":515900,"isOrg":false,"description":"","profileUrl":"/people/astrojhgu","avatar":{"id":"da8e974dc","template":"/{id}_{size}.jpg"},"isOrgWhiteList":false,"badge":{"identity":null,"bestAnswerer":null}}},"Comment":{},"favlists":{}},"me":{},"global":{},"columns":{"next":{},"data-analysis":{"following":false,"canManage":false,"href":"/api/columns/data-analysis","name":"数据的背后","creator":{"slug":"astrojhgu"},"url":"/data-analysis","slug":"data-analysis","avatar":{"id":"v2-e53fdbd4cec6a21ded9551d8cfa10c84","template":"/{id}_{size}.jpg"}}},"columnPosts":{},"columnSettings":{"colomnAuthor":[],"uploadAvatarDetails":"","contributeRequests":[],"contributeRequestsTotalCount":0,"inviteAuthor":""},"postComments":{},"postReviewComments":{"comments":[],"newComments":[],"hasMore":true},"favlistsByUser":{},"favlistRelations":{},"promotions":{},"switches":{"couldAddVideo":false},"draft":{"titleImage":"","titleImageSize":{},"isTitleImageFullScreen":false,"canTitleImageFullScreen":false,"title":"","titleImageUploading":false,"error":"","content":"","draftLoading":false,"globalLoading":false,"pendingVideo":{"resource":null,"error":null}},"drafts":{"draftsList":[],"next":{}},"config":{"userNotBindPhoneTipString":{}},"recommendPosts":{"articleRecommendations":[],"columnRecommendations":[]},"env":{"isAppView":false,"appViewConfig":{"content_padding_top":128,"content_padding_bottom":56,"content_padding_left":16,"content_padding_right":16,"title_font_size":22,"body_font_size":16,"is_dark_theme":false,"can_auto_load_image":true,"app_info":"OS=iOS"},"isApp":false},"sys":{}}}