MCMC、蒙特卡洛近似和Metropolis算法简介

简介: MCMC、蒙特卡洛近似和Metropolis算法简介

MCMC 是Markov Chain Monte Carlo 的简称,但在传统模拟中有一个很重要的假设是样本是独立的(independent samples),这一点在贝叶斯统计尤其是高纬度的模型中很难做到。所以MCMC的目的就是运用蒙特卡洛模拟出一个马可链(Markov chain)。

640.jpg

如今,概率建模风靡一时,但是当我第一次了解它时,总有一件事情困扰我。许多贝叶斯建模方法都需要计算积分,而我看到的任何工作示例似乎都使用高斯或伯努利分布,原因很简单如果您尝试使用比这更复杂的方法,它将成为分析的噩梦。将贝叶斯模型限制在“表现良好”的分布的小子集中,可能会极大地阻碍你对问题建模的能力,所以我们必须找到克服这一限制的方法。

蒙特卡洛近似

如果我不想分析计算某个讨厌的积分怎么办?可以使用蒙特卡洛近似。

我们知道,我们可以通过使用目标分布的样本值计算期望通过使用目标分布的样本值计算样本均值。为什么重要?那么,期望是什么呢?

640.jpg

连续随机变量的期望。同样的过程也适用于离散的情况,只要改变求和的积分。

这种估计积分的方法由中心极限定理提供了一些很好的保证。首先,这是期望的无偏估计,其次,我们可以计算估计的方差。

640.jpg

使用蒙特卡罗样本计算积分是非常好的,但是我们如何从目标分布中抽取样本呢?绘制高斯或均匀样本很容易,但np.random会让你失望。

画样本最简单的方法是使用逆CDF方法但这依赖于获得逆CDF函数它通常没有一个很好的解析形式只对一维随机变量有意义。

Metropolis算法是许多马尔可夫链蒙特卡洛(MCMC)采样方法的组成部分之一。当您可以访问的只是目标分布的pdf时,它使我们能够绘制样本。

MCMC方法需要注意的是我们不再采用独立样本所以我们不能保证估计的方差如何随着样本数量的增加而减少。如果样本是独立的,中心极限定理告诉我们估计值的方差将与样本数量(N)成反比地减少。对于MCMC,我们可以通过将样本数量从N调整为N_eff来对其贴现。N_eff(几乎)总是小于N,与链中样本的相关性有关。

Metropolis采样

Metropolis算法的步骤如下:

1.从目标分布域或先前分布的域中均匀采样起点。

2.在那时pdf。

3.根据一些状态转换函数,从当前位置走一步,为新样本提出建议。

4.计算新的pdf值。

5.计算新pdf /旧pdf的值。

6.如果比率大于1,请接受该步骤。

7.如果比率小于1:

     1.采样一个统一的随机变量。

     2.如果比率大于随机样本,请接受该步骤。

8.否则拒绝该建议,将当前职位作为示例添加并采取新的步骤。

请注意,在5–8中描述的过程等效于以概率为min(1,p(new)/ p(old))的伯努利概率接受样本,记住这个后面会用到……

Metropolis为何起作用?

对于任何MCMC方法,我们都希望确保一种称为详细平衡或可逆性的属性。

640.jpg

详细平衡

如果π满足,则π是马尔可夫链(1)的平稳分布。我们可以通过对等式两边求和来证明这一点。如果我们可以保证详细的平衡,那么我们也知道我们正在从马尔可夫链的固定分布中取样,我们将其作为目标分布。

详细平衡的思想是,由于每次转移的概率“质量”从状态i到状态j的转移与从状态j到状态i的转移是相同的,因此在链的每次转移之后,我们保持在平稳分布 。

现在,让我们展示一下Metropolis算法如何满足这一条件……

为了找到满足详细平衡的p(i,j),我们首先提出任意的“跳跃”概率q(i,j),然后仅通过接受概率为α(的“跳跃”来获得p(i,j)。i,j)。当“跳转”被拒绝时,状态保持j = i。这种“接受”的想法并不是Metropolis算法独有的,它存在于MCMC的大多数变体中。

640.jpg

跳跃概率推导

这取决于α是有效概率分布。因此,α的有效形式为:

640.jpg

Metropolis-Hastings 跳跃概率

如果跳跃概率是对称的,则可以简化为:

640.jpg

否则,它可以保留为完整形式,称为Metropolis-Hasting MCMC。

现在我们可以保证详细的平衡,我们可以让马尔可夫链式接管。如果马尔可夫链是遍历的(所有状态都是不可约的),那么在某个时候,该链将到达平稳分布,并且我们能够从目标分布中获取样本。

你可能还注意到,由于alpha是π(j)/π(i)的函数。这样的结果意味着目标分布无需标准化。当将Metropolis用于难以计算证据项的贝叶斯后验估计时,这特别有用。

实现的注意事项

Metropolis算法的通用版本称为“随机行走Metropolis”,其中建议的状态为当前状态,再加上均值为零且协方差矩阵为σ²I的多元高斯。σ应选择为足够使得足够多的样本被拒绝大。这是为了确保对目标分布进行良好的探索。

要注意的第二件事是老化的概念。在马尔可夫链到达平稳分布之前采集的样本应删除,因为它们在链收敛之前不能代表目标分布。确定要删除的样本数量有些困难,但通常,要删除的样本为10–20%(或有效的样本为10–100)。

Numpy代码实现

defmetropolis(pi, dims, n_samples, burn_in=0.1, var=1):
theta_=np.random.randn(dims)*varsamples= []
whilelen(samples) <n_samples:
theta=theta_+np.random.randn(dims)*varratio=pi(theta)/pi(theta_)
ifnp.random.rand(1) <ratio:
sample=thetatheta_=thetasamples.append(sample)
else:
sample=theta_samples.append(sample)   samples=np.array(samples)
returnsamples[int(samples*burn_in):,:]

我们可以看到这在两个高斯函数的和上的表现(注意这是一个非正态分布)。

fromscipy.statsimportmultivariate_normaldefmake_pdf(mean1, mean2, cov1, cov2):   pdf1=multivariate_normal(mean1, cov1)
pdf2=multivariate_normal(mean2, cov2)   defpdf(x):
returnpdf1.pdf(x) *pdf2.pdf(x)   returnpdfpdf=make_pdf(
    [3, 3],
    [-1, -1],
np.array([[1,0.1],[0.1,1]], dtype=float),
np.array([[1,0.1],[0.1,1]], dtype=float),
)
samples=metropolis(pdf, 2, 1_000, 0.1)

640.gif

上面的gif显示了算法是如何遍历分布的,偶尔会在分布的两种不同模式之间跳转。注意,这也突出了metropolis算法的一个弱点,它处理相对较差的多模型分布。这是由于要探索一种新模式,步骤必须足够大,以便从一种模式希望到另一种模式。这要么需要一个大的步长,要么需要一个模态紧密分布。修改如哈密顿量MCMC可以帮助解决这一问题,但一般来说,这是大多数MCMC方法的一个问题。

目录
相关文章
|
6月前
|
机器学习/深度学习 算法 TensorFlow
机器学习算法简介:从线性回归到深度学习
【5月更文挑战第30天】本文概述了6种基本机器学习算法:线性回归、逻辑回归、决策树、支持向量机、随机森林和深度学习。通过Python示例代码展示了如何使用Scikit-learn、statsmodels、TensorFlow库进行实现。这些算法在不同场景下各有优势,如线性回归处理连续值,逻辑回归用于二分类,决策树适用于规则提取,支持向量机最大化类别间隔,随机森林集成多个决策树提升性能,而深度学习利用神经网络解决复杂模式识别问题。理解并选择合适算法对提升模型效果至关重要。
239 4
|
2月前
|
算法 Java 数据安全/隐私保护
国密加密算法简介
国密指国家密码局认定的国产密码算法,主要包括SM1、SM2、SM3、SM4等,并持续完善。SM1是对称加密算法,加密强度与AES相当,需加密芯片支持;SM2是非对称加密,基于ECC算法,签名和密钥生成速度优于RSA;SM3为杂凑算法,安全性高于MD5;SM4为对称加密算法,用于无线局域网标准。本文提供使用Java和SpringBoot实现SM2和SM4加密的示例代码及依赖配置。更多国密算法标准可参考国家密码局官网。
110 1
|
27天前
|
存储 算法 安全
ArrayList简介及使用全方位手把手教学(带源码),用ArrayList实现洗牌算法,3个人轮流拿牌(带全部源码)
文章全面介绍了Java中ArrayList的使用方法,包括其构造方法、常见操作、遍历方式、扩容机制,并展示了如何使用ArrayList实现洗牌算法的实例。
12 0
|
2月前
|
机器学习/深度学习 数据采集 存储
一文读懂蒙特卡洛算法:从概率模拟到机器学习模型优化的全方位解析
蒙特卡洛方法起源于1945年科学家斯坦尼斯劳·乌拉姆对纸牌游戏中概率问题的思考,与约翰·冯·诺依曼共同奠定了该方法的理论基础。该方法通过模拟大量随机场景来近似复杂问题的解,因命名灵感源自蒙特卡洛赌场。如今,蒙特卡洛方法广泛应用于机器学习领域,尤其在超参数调优、贝叶斯滤波等方面表现出色。通过随机采样超参数空间,蒙特卡洛方法能够高效地找到优质组合,适用于处理高维度、非线性问题。本文通过实例展示了蒙特卡洛方法在估算圆周率π和优化机器学习模型中的应用,并对比了其与网格搜索方法的性能。
267 1
|
4月前
|
机器学习/深度学习 人工智能 自然语言处理
算法金 | 秒懂 AI - 深度学习五大模型:RNN、CNN、Transformer、BERT、GPT 简介
**RNN**,1986年提出,用于序列数据,如语言模型和语音识别,但原始模型有梯度消失问题。**LSTM**和**GRU**通过门控解决了此问题。 **CNN**,1989年引入,擅长图像处理,卷积层和池化层提取特征,经典应用包括图像分类和物体检测,如LeNet-5。 **Transformer**,2017年由Google推出,自注意力机制实现并行计算,优化了NLP效率,如机器翻译。 **BERT**,2018年Google的双向预训练模型,通过掩码语言模型改进上下文理解,适用于问答和文本分类。
144 9
|
4月前
|
算法
Raid5数据恢复—Raid5算法简介&raid5磁盘阵列数据恢复案例
Raid5算法也被称为“异或运算”。异或是一个数学运算符,它应用于逻辑运算。异或的数学符号为“⊕”,计算机符号为“xor”。异或的运算法则为:a⊕b = (¬a ∧ b) ∨ (a ∧¬b)。如果a、b两个值不相同,则异或结果为1。如果a、b两个值相同,异或结果为0。 异或也叫半加运算,其运算法则相当于不带进位的二进制加法。二进制下用1表示真,0表示假。异或的运算法则为:0⊕0=0,1⊕0=1,0⊕1=1,1⊕1=0(同为0,异为1),这些法则与加法是相同的,只是不带进位。 异或略称为XOR、EOR、EX-OR,程序中有三种演算子:XOR、xor、⊕。使用方法如下z = x ⊕ y z
Raid5数据恢复—Raid5算法简介&raid5磁盘阵列数据恢复案例
|
3月前
|
算法
【算法】贪心算法简介
【算法】贪心算法简介
|
3月前
|
算法
【算法】递归、搜索与回溯——简介
【算法】递归、搜索与回溯——简介
|
6月前
|
算法 关系型数据库 C语言
卡尔曼滤波简介+ 算法实现代码(转)
卡尔曼滤波简介+ 算法实现代码(转)
74 4
|
6月前
|
算法 项目管理
R语言实现蒙特卡洛模拟算法
R语言实现蒙特卡洛模拟算法