转载自：http://site.douban.com/182577/widget/notes/10567181/note/280112466/

####贝叶斯集锦这个系列目的是想收集一些使用R的贝叶斯统计案例
####没什么系统性也没啥方法介绍，大家见谅。

MCMCpack包的一个例子：贝叶斯线性回归

从1960年代贝叶斯统计学派复兴到今天贝叶斯统计的广泛运用，MCMC方法起到了极重要的作用。作为一种计算手段，MCMC以模拟的方法解决了贝叶斯方法中后验分布的计算问题。

MCMCpack包提供了（主要对社会科学研究）进行贝叶斯推断和贝叶斯计算的工具（特别是MCMC）。 MCMCpack包的设计思想是针对特定的模型运用MCMC方法。 MCMCpack 当前包括18种统计模型: 线性回归模型 (linear regression with Gaussian errors, a singular value decomposition regression, and regression for a censored dependent variable), 离散选择模型 (logistic regression, multinomial logistic regression, ordinal probit regression, and probit regression), 测量模型(a one-dimensional IRT model, a k-dimensional IRT model, a k-dimensional ordinal factor model, a k-dimensional linear factor model, a k-dimensionl mixed factor model, and a k- dimensional robust IRT model), 计数数据模型(a Poisson regression model), 生态推断模型(a hierarchical ecological inference model and a dynamic ecological inference model), 变点问题的时间序列模型 (a binary change-pointmodel, a probit change-point model, an ordinal probit change-point model, and a Poisson change-point model). 其中特别是测量模型尤其适合采用MCMC处理。这个包同时还包括了下列分布Dirichlet, inverse Gamma, inverse Wishart, noncentral Hypergeometric, Wishart 的概率密度和伪随机数生成器。

MCMCpack包在函数的设计上尽可能的和相应问题的其它R函数在语法上一致。 以线性回归为例，大家很熟悉： lm(y ~ x1 + x2 + x3, data = mydata)， MCMCpack包的相应的函数是： MCMCregress(y ~ x1 + x2 + x3, data = mydata)。

下面讨论一个例子。

1.贝叶斯线性回归 
来自R包dataset的数据集swiss。这个数据集是1888年对瑞士的47个法语地区的社会经济指标的调查。

library(MCMCpack)

## Loading required package: coda
## Loading required package: lattice
## ## ## Markov Chain Monte Carlo Package (MCMCpack)

data(swiss)
swiss.posterior1 <- MCMCregress(Fertility ~ Agriculture + Examination + Education + 
    Catholic + Infant.Mortality, data = swiss)
summary(swiss.posterior1)

## 
## Iterations = 1001:11000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
## 
## Mean SD Naive SE Time-series SE
## (Intercept) 67.021 11.0813 0.110813 0.111734
## Agriculture -0.172 0.0731 0.000731 0.000731
## Examination -0.259 0.2606 0.002606 0.002606
## Education -0.872 0.1892 0.001892 0.001892
## Catholic 0.104 0.0360 0.000360 0.000360
## Infant.Mortality 1.074 0.3958 0.003958 0.003819
## sigma2 54.050 12.6860 0.126860 0.145950
## 
## 2. Quantiles for each variable:
## 
## 2.5% 25% 50% 75% 97.5%
## (Intercept) 45.5320 59.5653 67.060 74.3160 88.871
## Agriculture -0.3179 -0.2212 -0.171 -0.1236 -0.027
## Examination -0.7659 -0.4306 -0.258 -0.0862 0.249
## Education -1.2428 -0.9983 -0.871 -0.7454 -0.499
## Catholic 0.0315 0.0801 0.104 0.1276 0.175
## Infant.Mortality 0.2859 0.8167 1.073 1.3377 1.855
## sigma2 34.5771 45.0633 52.351 61.0374 83.851

MCMCregresss默认的参数值是无信息先验（采取这样的先验，参数估计的结果（在贝叶斯回归里就是我们所求的参数均值），和最小二乘法求得的回归系数是一致的）。 MCMC包的模型函数返回由coda包所定义的mcmc对象（见str(swiss.posterior1)）。MCMCpack要依赖coda来实现后验分布的摘要和模拟数据的收敛诊断。

摘要包括后验均值，标准差，均值标准误以及分位数。

作为对比，用最小二乘法的普通线性回归：

swiss.lm <- lm(Fertility ~ Agriculture + Examination + Education + Catholic + 
    Infant.Mortality, data = swiss)
summary(swiss.lm)

## 
## Call:
## lm(formula = Fertility ~ Agriculture + Examination + Education + 
## Catholic + Infant.Mortality, data = swiss)
## 
## Residuals:
## Min 1Q Median 3Q Max 
## -15.274 -5.262 0.503 4.120 15.321 
## 
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) 
## (Intercept) 66.9152 10.7060 6.25 1.9e-07 ***
## Agriculture -0.1721 0.0703 -2.45 0.0187 * 
## Examination -0.2580 0.2539 -1.02 0.3155 
## Education -0.8709 0.1830 -4.76 2.4e-05 ***
## Catholic 0.1041 0.0353 2.95 0.0052 ** 
## Infant.Mortality 1.0770 0.3817 2.82 0.0073 ** 
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.17 on 41 degrees of freedom
## Multiple R-squared: 0.707, Adjusted R-squared: 0.671 
## F-statistic: 19.8 on 5 and 41 DF, p-value: 5.59e-10

用coda包的plot命令，可以表现模拟trace图和后验的分布密度。

plot(swiss.posterior1[, 1:2], col = 4)

 

2.模型选择

MCMCpack包采用贝叶斯因子进行模型选择。 设数据y可以生成模型M1和M2，则贝叶斯因子BF12=P(y|M1)/P(y|M2)。

swiss.posterior1 <- MCMCregress(Fertility ~ Agriculture + Examination + Education + 
    Catholic + Infant.Mortality, data = swiss, marginal.likelihood = "Chib95", 
    b0 = 0, B0 = 0.1, c0 = 2, d0 = 0.11)
# 对 improper prior不能计算marginal.likelihood，所以了取了一个weakly
# #informative prior

swiss.posterior2 <- MCMCregress(Fertility ~ Agriculture + Education + Catholic + 
    Infant.Mortality, data = swiss, marginal.likelihood = "Chib95", b0 = 0, 
    B0 = 0.1, c0 = 2, d0 = 0.11)

swiss.posterior3 <- MCMCregress(Fertility ~ Agriculture + Examination + Catholic + 
    Infant.Mortality, data = swiss, marginal.likelihood = "Chib95", b0 = 0, 
    B0 = 0.1, c0 = 2, d0 = 0.11)
bf <- BayesFactor(swiss.posterior1, swiss.posterior2, swiss.posterior3)

第2个模型最好，第2个模型的摘要如下：

summary(swiss.posterior2)

## 
## Iterations = 1001:11000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
## 
## Mean SD Naive SE Time-series SE
## (Intercept) 3.7080 3.1297 0.031297 0.032899
## Agriculture 0.1064 0.0733 0.000733 0.000733
## Education -0.4672 0.1692 0.001692 0.001692
## Catholic 0.0784 0.0376 0.000376 0.000376
## Infant.Mortality 3.1195 0.2586 0.002586 0.002586
## sigma2 94.5259 21.5607 0.215607 0.247799
## 
## 2. Quantiles for each variable:
## 
## 2.5% 25% 50% 75% 97.5%
## (Intercept) -2.400 1.6036 3.711 5.800 9.803
## Agriculture -0.038 0.0564 0.107 0.155 0.249
## Education -0.800 -0.5814 -0.466 -0.352 -0.134
## Catholic 0.006 0.0529 0.078 0.104 0.154
## Infant.Mortality 2.619 2.9446 3.118 3.293 3.631
## sigma2 61.204 79.2795 91.626 106.441 143.902

贝叶斯集锦：MCMCpack包相关推荐

  1. 贝叶斯集锦：从MC、MC到MCMC

    转载自: #####一份草稿 贝叶斯计算基础 一.从MC.MC到MCMC 斯坦福统计学教授Persi Diaconis是一位传奇式的人物.Diaconis14岁就成了一名魔术师,为了看懂数学家Fell ...

  2. 贝叶斯集锦：R和JAGS的交互

    转载自:http://site.douban.com/182577/widget/notes/10567181/note/295466672/ Markov chain Monte Carlo (MC ...

  3. 贝叶斯集锦：贝叶斯统计基础

    转载自:http://site.douban.com/182577/widget/notes/10567181/note/294041203/ 1.从贝叶斯定理到贝叶斯统计推断 (1)贝叶斯统计简史 ...

  4. 贝叶斯优化python包_贝叶斯全局优化（LightGBM调参）

    这里结合Kaggle比赛的一个数据集,记录一下使用贝叶斯全局优化和高斯过程来寻找最佳参数的方法步骤. 1.安装贝叶斯全局优化库 从pip安装最新版本 pip install bayesian-opti ...

  5. 贝叶斯集锦：贝叶斯派和频率派的一个例子

    转载自:http://site.douban.com/182577/widget/notes/10567181/note/278503359/ 这个例子的主要目的在于探讨贝叶斯派和频率派适用的具体情境 ...

  6. 贝叶斯优化python包_《用贝叶斯优化进行超参数调优》

    TPE CMAES 网格搜索 随机搜索 贝叶斯优化 用贝叶斯优化进行超参数调优 @QI ZHANG · JUL 12, 2019 · 7 MIN READ 超参数调优一直是机器学习里比较intract ...

  7. 贝叶斯优化python包_贝叶斯优化

    万壑松风知客来,摇扇抚琴待留声 1. 文起 本篇文章记录通过 Python 调用第三方库,从而调用使用了贝叶斯优化原理的 Hyperopt 方法来进行超参数的优化选择.具体贝叶斯优化原理与相关介绍将在 ...

  8. 朴素贝叶斯关于naivebayes包核函数等

    关于naivebayes包: 默认 S3 方法: naive_bayes(x, y, prior = NULL, laplace = 0, usekernel = FALSE, usepoisson ...

  9. 贝叶斯优化python包_Bayesian2D-用贝叶斯优化方法求任意二维函数的最大值或最小值的软件包-Juhan Raidal...

    作者:Juhan Raidal ### 作者邮箱:juhantiitus.raidal@gmail.com ### 首页:https://github.com/JRaidal/Bayesian2D # ...

最新文章

  1. 神经网络入门——14多层感知机
  2. [Ubuntu 12.10] Openstack 多节点安装--前期准备网络拓扑
  3. 攻击链路识别——CAPEC（共享攻击模式的公共标准）、MAEC（恶意软件行为特征）和ATTCK（APT攻击链路上的子场景非常细）...
  4. mybatis中ResultSetHandler的设计与实现
  5. 浏览器获取设备信息_过滤获取日志和浏览器信息
  6. wxWidgets：使用文本模板
  7. python数据预测_python时间序列预测股票走势
  8. 用户变量和系统变量的区别_环境变量的用户变量与系统变量的区别
  9. win8 无法打开任务管理器
  10. CVTE 2017 秋季校招笔试题回忆（C++后台）
  11. 《Spring Data实战》——2.2　定义查询方法
  12. visual studio编Java,如何用Visual Studio编译Java源代码
  13. 以贝叶斯思维看待世界
  14. CSDN 上传资源已经存在
  15. golang cond
  16. 3D折纸效果怎么实现？
  17. 51单片机74ls273并行输出地址c语言程序,跑马灯/输入输出接口（片选地址74LS273）...
  18. 阿里首推“数据安全合作伙伴计划” 构建数据安全生态
  19. C++ 演讲比赛流程管理系统
  20. php导购系统,php源码：鱼福CMS淘客导购系统V1.2–免费开源

热门文章

  1. 2021音视频技术大会北京站开幕
  2. Python调用自己写的模块
  3. 计算机网络的三种通讯模式（单播，广播，组播）小结
  4. Termux第一篇之ssh使用
  5. Ubuntu14.04 安装pip
  6. 深入剖析Android音频(三)AudioPolicyService
  7. bat之启动与禁用网卡
  8. lisp 读取样条曲线座标点_MATLAB插值绘制曲线
  9. array转list_Java面试题Array和ArrayList有何区别？
  10. ubuntu linux配置bond 网卡绑定 多个bond配置多网关