R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯估计

MCMC是从复杂概率模型中采样的通用技术。

由Kaizong Ye,Weilong Zhang撰写

1.蒙特卡洛

2.马尔可夫链

3.Metropolis-Hastings算法

如果需要计算有复杂后验pdf p(θ| y)的随机变量θ的函数f(θ)的平均值或期望值


您可能需要计算后验概率分布p(θ)的最大值。

 

解决期望值的一种方法是从p(θ)绘制N个随机样本,当N足够大时,我们可以通过以下公式逼近期望值或最大值

将相同的策略应用于通过从p(θ| y)采样并取样本集中的最大值来找到argmaxp(θ| y)。


解决方法

1.1直接模拟

1.2逆CDF

1.3拒绝/接受抽样

如果我们不知道精确/标准化的pdf或非常复杂,则MCMC会派上用场。


热门课程

R语言数据分析挖掘必知必会

面对扑面而来的数据浪潮,包含Google、Facebook等国际企业,都已采用R语言进行数据分析

探索课程

马尔可夫链

为了模拟马尔可夫链,我们必须制定一个 过渡核T(xi,xj)。过渡核是从状态xi迁移到状态xj的概率。

 马尔可夫链的收敛性意味着它具有平稳分布π。马尔可夫链的统计分布是平稳的,那么它意味着分布不会随着时间的推移而改变。

Metropolis算法

 对于一个Markov链是平稳。基本上表示

处于状态x并转换为状态x’的概率必须等于处于状态x’并转换为状态x的概率

 

或者

方法是将转换分为两个子步骤;候选和接受拒绝。

令q(x’| x)表示 候选密度,我们可以使用概率 α(x’| x)来调整q  。

候选分布 Q(X’| X)是给定的候选X的状态X’的条件概率,

和 接受分布 α(x’| x)的条件概率接受候选的状态X’-X’。我们设计了接受概率函数,以满足详细的平衡。

该 转移概率 可以写成:

 

即,当接受度大于1时,我们总是接受,而当接受度小于1时,我们将相应地拒绝。因此,Metropolis-Hastings算法包含以下内容:

  1. 初始化:随机选择一个初始状态x;
  2. 根据q(x’| x)随机选择一个新状态x’;

3.接受根据α(x’| x)的状态。如果不接受,则不会进行转移,因此无需更新任何内容。否则,转移为x’;

4.转移到2,直到生成T状态;

5.保存状态x,执行2。

原则上,我们从分布P(x)提取保存的状态,因为步骤4保证它们是不相关的。必须根据候选分布等不同因素来选择T的值。 重要的是,尚不清楚应该使用哪种分布q(x’| x);必须针对当前的特定问题进行调整。


属性

Metropolis-Hastings算法的一个有趣特性是它 仅取决于比率

是候选样本x’与先前样本xt之间的概率,

是两个方向(从xt到x’,反之亦然)的候选密度之比。如果候选密度对称,则等于1。


使用R语言进行Metroplis-in-Gibbs采样和MCMC运行分析

阅读文章


马尔可夫链从任意初始值x0开始,并且算法运行多次迭代,直到“初始状态”被“忘记”为止。这些被丢弃的样本称为预烧(burn-in)。其余的x可接受值集代表分布P(x)中的样本

 


Metropolis采样

一个简单的Metropolis-Hastings采样

让我们看看从 伽玛分布 模拟任意形状和比例参数,使用具有Metropolis-Hastings采样算法。

下面给出了Metropolis-Hastings采样器的函数。该链初始化为零,并在每个阶段都建议使用N(a / b,a /(b * b))个候选对象。

基于正态分布且均值和方差相同gamma的Metropolis-Hastings独立采样

  1. 从某种状态开始xt。代码中的x。
  2. 在代码中提出一个新的状态x’候选
  3. 计算“接受概率”
  4. 从[0,1] 得出一些均匀分布的随机数u;如果u <α接受该点,则设置xt + 1 = x’。否则,拒绝它并设置xt + 1 = xt。


随时关注您喜欢的主题


MH可视化


画图

设置参数。

修改图,仅包含预烧期后的链

初始值

第一个样本 vec 是我们链的初始/起始值。我们可以更改它,以查看收敛是否发生了变化。


选择方案

如果候选密度与目标分布P(x)的形状匹配,即q(x’| xt)≈P(x’)q(x’|),则该算法效果最佳。 xt)≈P(x’)。如果使用正态候选密度q,则在预烧期间必须调整方差参数σ2。

通常,这是通过计算接受率来完成的,接受率是在最后N个样本的窗口中接受的候选样本的比例。

如果σ2太大,则接受率将非常低,因为候选可能落在概率密度低得多的区域中,因此a1将非常小,且链将收敛得非常慢。

示例2:回归的贝叶斯估计

Metropolis-Hastings采样用于贝叶斯估计回归模型。


设定参数


DGP和图


正态分布拟然


为什么使用对数

我们为什么要做这个?强烈建议这样做,因为许多小概率相乘的概率会变得很小。在某个阶段,计算机程序会陷入数值四舍五入或下溢问题。

似然函数中概率的对数,这也是我求和所有数据点的概率(乘积的对数等于对数之和)的原因。

因此, 当您编写概率时,请始终使用对数


示例:绘制斜率a的似然曲线


先验分布

这三个参数的均匀分布和正态分布。


后验

先验和概率的乘积是MCMC将要处理的实际量。此函数称为后验函数。同样,这里我们使用和,因为我们使用对数。


Metropolis算法

该算法是 后验密度采样最常见的贝叶斯统计应用之一 。

上面定义的后验。

  1. 从随机参数值开始
  2. 根据某个候选函数的概率密度,选择一个接近旧值的新参数值
  3. 以概率p(new)/ p(old)跳到这个新点,其中p是目标函数,并且p> 1也意味着跳跃
  4. 请注意,我们有一个 对称的跳跃/候选分布 q(x’| x)。

标准差σ是固定的。

所以接受概率等于

实施

(e)输出接受的值,并解释。

算法的第一步可能会因初始值而有偏差,因此通常会被丢弃来进行进一步分析(预烧期)。令人感兴趣的输出是接受率:候选多久被算法接受拒绝一次?候选函数会影响接受率:通常,候选越接近,接受率就越大。但是,非常高的接受率通常是无益的:这意味着算法在同一点上“停留”,这导致对参数空间(混合)的处理不够理想。

我们还可以更改初始值,以查看其是否更改结果/是否收敛。


小结


总结:


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

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

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

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498