R语言使用Metropolis-Hastings采样算法自适应贝叶斯估计与可视化

By 大数据部落, 数理统计, 机器学习, 计算机科学与技术, , , , , , , ,

如果您可以写出模型的似然函数，则 Metropolis-Hastings算法可以负责其余部分（即MCMC ）。

WeChat Tencent QQ email print

由Kaizong Ye，Liao Bao撰写

我写了r代码来简化对任意模型的后验分布的估计。具体如下：


我写了r代码来简化对任意模型的后验分布的估计。具体如下：

可下载资源

     完整程序、数据和文档（word） 

 
 
a<-pars[1]      #截距
 
b<-pars[2]      #斜率
 
sd_e<-pars[3]   #残差
 
if(sd_e<=0){return(NaN)}
 
 
 
log_likelihood<-sum( dnorm(data[,2],pred,sd_e, log=TRUE) )
×

1.1 基本概念

  • 从一个概率分布（目标分布[公式]）中得到随机样本序列

  • 这个序列可以用于：a) 近似估计分布[公式]; b) 计算积分(期望)

  • 用于高维分布取样

  • 缺陷: MCMC固有缺点, 样本自相关性

1.2 优势

  • 可以从任意的概率分布[公式]中取样，只要满足条件：函数[公式]成比例于[公式]的密度。

  • 更宽松的要求：[公式]仅需要与[公式]的密度成比例。

1.3 要点

  • 生成样本值序列；样本值产生得越多，这些值的分布就越近似于[公式]

  • 迭代产生样本值：下一个样本的分布仅仅取决于当前样本值（马尔科夫链特性）

  • 接受/拒绝概率：接受计算出的值为下一个样本值/拒绝并重复使用当前样本值；基于[公式]，接受概率通过比较[公式]（当前值）和[公式]（备选值）得出

1.4 Metropolis Algorithm (对称分布)

  • Input: [公式]，与目标分布[公式]成比例的函数


1. 初始化：

  • [公式]: 取任意点作为第一个样本值

  • [公式]: 任意的概率密度函数；给定上一个样本值[公式]，计算出下一个样本值[公式]； 一般我们设定[公式]是对称的，即[公式]；我们也会令[公式]服从以[公式]为中心的[公式]分布，即越靠近[公式]，概率越大

  • 将样本序列制成随机游走（random walk）；[公式]: 假定密度/跳跃分布(jumping distribution)


2. 迭代[公式]

  • 通过密度函数[公式]生成备选样本值[公式]

  • 计算接受备选样本值的比率: [公式]；由于[公式]又与[公式]成比例，因此[公式]。

  • 若[公式]，[公式]比[公式]更接近[公式]，令[公式]；否则，以[公式]的概率接受备选样本值[公式]；若备选样本值被拒绝，令[公式]。

1.5 缺陷

  • 样本自相关：相邻的样本会相互相关，虽然可以通过每[公式]步取样的方式来减少相关性，但这样的后果就是很难让样本近似于目标分布[公式]。

  1. 自相关性可通过增加步调长度（jumping width，与jumping function [公式]的方差有关）来控制，但同时也增加了拒绝备选样本的几率[公式]。

  2. 过大或过小的jumping size会导致slow-mixing Markov Chain, 即高度自相关的一组样本，以至于我们需要得到非常大的样本量[公式]才能得到目标分布[公式]。

  • 初始值的选择：尽管Markov Chain最后都会收敛到目标分布，然而初始值的选择直接影响到运算时间，尤其是把初始值选在在了“低密度”区域。因此，选择初值时，最好加入一个“burn-in period”（预烧期，预选期）。

1.6 优势

  • 抗“高维魔咒”（curse of dimensionality）：维度增加，对于rejection sampling方法来说，拒绝的概率就是呈指数增长。而MCMC则成为了解决这种问题的唯一方法

  • 多元分布中，为了避免多元初始值以及[公式]选择不当而导致的问题，Gibbs sampling是另外一个更适合解决多元分布问题的MCMC 方法。Gibbs sampling从多元分布的各个维度中分别选择初始值，然后这些变量分别同时取样。

1.7 衍生算法

Metropolis-Hastings算法的目的是根据目标分布[公式]生成一个状态集。算法运用了渐进于平稳分布[公式]的Markov过程使得[公式]。

Markov过程由它的状态转换概率所决定的。[公式]为从任意状态[公式]转化到另外一个任意状态[公式]的转换概率。当它符合以下两点是，它有唯一的平稳分布[公式]:

  • 平稳分布[公式]的存在性：充分但非必要条件为细致平衡（detailed balance）：每一个转换[公式]都是可逆的，即对任意连个状态[公式]和[公式]，在状态[公式]时[公式]的概率等于在状态[公式]时[公式]的概率，即[公式]

  • 平稳分布[公式]的唯一性：由Markov过程的遍历性可得，即每一个转换[公式]需满足：

    • 非周期性的（aperiodic）：系统不会再某一个固定时间区间内回到原状态

    • 正递归的（positive recurrent）：回到同一个状态的步数的期望是有限值

Metropolis-Hastings算法需要通过构建转换概率来构建一个Markov过程且需要满足以上两点来使得Markov过程的平稳分布[公式]就是目标分布[公式]。算法衍生从满足细致平衡条件开始：


[公式],

即 [公式] (1)

这样就把转换的细致平衡转化为两步：提议和接受-拒绝。

提议分布[公式]为给定状态[公式]时提出状态[公式]的条件概率，接受分布[公式]为接受所提出的状态[公式]的概率。转换概率就能写成它们的乘积：

[公式] (2)

由(1)和(2)就有：

[公式]

之后就要选一个满足上面条件的接受函数[公式]了。一般我们选Metropolis选项，即：

[公式]

满足了Metropolis算法的基本步骤，即大于1的时候总是接受，小于1的时候按比例接受和拒绝。

则Metropolis-Hastings算法可以概述为：

  1. 初始化：随机选取初始状态[公式]

  2. 根据[公式]随机选取状态[公式]

  3. 根据[公式]选择接受或拒绝状态[公式]；若接受，系统转换到状态[公式]，否则停留在状态[公式]

  4. 返回步骤2，直到生成了[公式]个状态

  5. 保存状态[公式]，返回步骤2


原则上说，被保存的状态是从分布[公式]中生成的，因为步骤4保证了这些状态是非相关的。取值[公式]需根据不同的特性选取，如提出的分布以及，它必须是Markov过程自相关时间的阶。

关于[公式]的选取需具体问题具体分析和调整。

先验：

 
epsilon<-pars[3]    #残差
 
 
 
prior_a<-dnorm(a,0,100,log=TRUE)     ##所有的非信息性先验
 
prior_b<-dnorm(b,0,100,log=TRUE)     ## 参数.
 
prior_epsilon<-dgamma(epsilon,1,1/100,log=TRUE)

1）现在让我们模拟一些数据以进行运行测试：

x<-runif(30,5,15)
 
y<-x+rnorm(30,0,5) ##斜率=1, 截距=0, epsilon=5

2）Metro Hastings 完成所有工作。

MH(li_func=li_reg,pars=c(0,1,1),

视频

贝叶斯推断线性回归与R语言预测工人工资数据

探索见解 ➜

去bilibili观看 ➜

探索更多视频 ➜

3）您可以使用plotMH（）查看所有模型参数的后验

plot(mcmc)

绘制所有参数之间的相关性。

4）输出后验置信区间。

BCI
 
#              0.025    0.975
 
# a       -5.3345970 6.841016
 
# b        0.4216079 1.690075
 
# epsilon  3.8863393 6.660037

接下来，我想提供一种直观的方法来可视化此算法运行的情况。

主要思想是从分布中抽取样本。积分很重要，贝叶斯定理本身：

P（θ| D）= P（D |θ）P（θ）/ P（D）

其中P（D）是观察数据的无条件概率。由于这不依赖于推断的模型（θ）参数，因此P（D）是归一化常数。

因此，我们有一个非归一化的概率密度函数，我们希望通过随机抽样来估计。对于复杂的模型而言，随机抽样本身的过程通常很困难，因此，我们使用马尔可夫链来探索分布。我们需要一个链，如果运行时间足够长，它将作为目标分布的随机样本整体。我们构建的马尔可夫链的这种特性称为 遍历性。Metropolis-Hastings算法是构建这种链的一种方法。

步骤

  1. 在参数空间k_X中选择一些起点
  2. 选择一个候选点k_Y〜N（k_X，σ）。这通常称为提议分布。
  3. 移至候选点的概率为：min（π（k_Y）/π（K_X），1）
  4. 重复。

最受欢迎的见解

1.使用R语言进行METROPLIS-IN-GIBBS采样和MCMC运行

2.R语言中的Stan概率编程MCMC采样的贝叶斯模型

3.R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样

4.R语言BUGS JAGS贝叶斯分析 马尔科夫链蒙特卡洛方法（MCMC）采样

5.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

6.R语言Gibbs抽样的贝叶斯简单线性回归仿真分析

7.R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数

8.R语言使用Metropolis- Hasting抽样算法进行逻辑回归

9.R语言中基于混合数据抽样(MIDAS)回归的HAR-RV模型预测GDP增长

以下代码通过简单的正态目标分布演示了此过程。

 
###     Metropolis-Hastings 可视化                #######
 
 
 
k_X = seed; ##将k_X设置为种子位置
 
 
 
 
for(i in 1:iter)
 
{
 
track<-c(track,k_X)    ## 链
 
k_Y = rnorm(1,k_X,prop_sd) ##候选点
 
 
 
## -- 绘制链的核密度估计
 
 
 
 
 
lines(density(track,adjust=1.5),col='red',lwd=2)
 
 
 
## -- 绘制链
 
 
plot(track,1:i,xlim=plot_range,main='',type='l',ylab='Trace')
 
 
 
 
## -- 绘制目标分布和提议分布 
 
 
 
curve(dnorm(x,k_X,prop_sd),col='black',add=TRUE)
 
abline(v=k_X,lwd=2)
 
 
 
 
## 接受概率为a_X_Y 
 
if (log(runif(1))<=a_X_Y)
 
 
 
points(k_Y,0,pch=19,col='green',cex=2)
 
 
 
## 调整提议
 
if(i>100)
 
prop_sd=sd(track[floor(i/2):i])

R语言使用Metropolis- Hasting抽样算法进行逻辑回归

阅读文章 ➜

该算法实现中的一个普遍问题是σ的选择。当σ接近目标分布的标准偏差时，将发生有效混合（链收敛到目标分布）。当我们不知道这个值时。我们可以允许σ根据到目前为止的链历史记录进行调整。在上面的示例中，将σ更新为链中某些先验点的标准偏差值。

输出为多页pdf，可以滚动浏览。

随时关注您喜欢的主题

注意：请注意，前100次左右的迭代是目标分布的较差表示。在实践中，我们将“预烧”该链的前n个迭代-通常是前100-1000个。

顶部显示了目标分布（蓝色虚线）和通过MCMC样本对目标进行的核平滑估计。第二面板显示了链的轨迹，底部显示了算法本身的步骤。

可下载资源

关于作者

Kaizong Ye是拓端研究室（TRL）的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文，如需帮助请联系我们！

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498
微信扫一扫
关注拓端助手

关注有关新文章的微信公众号

永远不要错过任何见解。当新文章发表时，我们会通过微信公众号向您推送。

技术干货

最新洞察

This will close in 0 seconds