MCMC抽样---Metropolis算法

MCMC抽样—Metropolis算法

马尔科夫链蒙特卡洛抽样方法可追溯到1953年N.Metropolis等人在研究原子和分子的随机性运动问题时所引入的随机模拟方法。该方法被命名为Metropolis模拟算法，这个算法已被列为影响科学和工程技术发展的最伟大的十大算法之首。

Metropolis算法是MCMC的核心。MCMC的基本思想是构造一个遍历的马尔科夫链，使得其不变分布成为人们所需要的抽样分布。做到这点似乎相当复杂，但实际上由于人们可以非常灵活地选择简单的转移概率，所以构造该算法并不困难。

Metropolis算法的动机是接受-拒绝抽样方法。在接受-拒绝抽样方法中，我们有布标概率密度函数f，建议概率密度函数g和一个接受准则h(x,y)。同样地我们假设：f是全空间Ω\OmegaΩ上的目标概率密度函数，人们需要在Ω\OmegaΩ上产生样本马尔科夫链x1,x2,⋯,xt,⋯{x_1,x_2,\cdots,x_t,\cdots}x1,x2,⋯,xt,⋯，使得它的稳定分布恰好是这个目标概率密度函数f的分布。类似于接受-拒绝方法，在当前状态xtx_txt下产生下一个状态分两个步骤进行：
（1）从建议概率密度函数g(.∣xt)g(.|x_t)g(.∣xt)产生一个随机数y作为建议的下一个状态（注意，这是依赖当前状态的条件概率）；
（2）然后按均匀分布生成一个随机数r，如果r<h(x,y)r\lt h(x,y)r<h(x,y)则接受该建议的随机数，即xt+1=yx_{t+1}=yxt+1=y，否则放弃y而采用原状态xt+1=xtx_{t+1}=x_txt+1=xt。重复这个过程产生随机序列。很显然，这个过程所产生的随机序列做到了下一随机数仅依赖于当前的随机数，而和以前产生的随机数无关。

显然，在算法中建议的概率密度函数g主要目的是用来产生状态转移，即为每个状态构造出一个领域，并使邻域中的某个邻居倍选中。要使得它产生的序列是马尔科夫链，这要求建议概率密度函数g必须在整个Ω\OmegaΩ上都有定义，这通常是容易做到的。接受-拒绝抽样方法是最简单的情况：对所有的x,g(.∣x)=g(.)x,g(.|x)=g(.)x,g(.∣x)=g(.)，单MCMC推广了它，使它变成了一个条件概率密度函数。当然必须要求在x的某个邻域NxN_xNx里有g(y∣x)>0,∀y∈Nxg(y|x)>0,\forall y \in N_xg(y∣x)>0,∀y∈Nx，这样才能够使状态产生遍历性的转移。一种特殊的Metropolis算法需要g满足对称性，即满足条件g(y∣x)=g(x∣y)g(y|x)=g(x|y)g(y∣x)=g(x∣y)，这种对称性在大多数情况下可以很自然地做到，当然它也可以稍微减弱成下面的条件g(y∣x)>0⟺g(x∣y)>0g(y|x)\gt 0 \iff g(x|y) \gt 0g(y∣x)>0⟺g(x∣y)>0，即要求状态转移是可逆的。

算法接下来要做的是，将选中的状态和当前状态进行比较，以一定的概率接受其中之一为随机序列链的下一状态。这里与接受-拒绝方法的另一不同之处在于，接受概率不但取决于下一状态y，而且还与当前状态x有关。具体的方法是所谓的Metropolis-Hastings算法，该算法的接受概率定义为
h(x,y)=min1,f(y)g(x∣y)f(x)g(y∣x)h(x, y)=min{1, \frac{f(y)g(x|y)}{f(x)g(y|x)}} h(x,y)=min1,f(x)g(y∣x)f(y)g(x∣y)

党对称性条件满足时，上式简化成
h(x,y)=min(1,f(y)f(x))h(x, y)=min(1,\frac{f(y)}{f(x)}) h(x,y)=min(1,f(x)f(y))

这就是说，如果建议的下一状态y具有比当前状态x的概率较大，即f(y)≥f(x)f(y)\ge f(x)f(y)≥f(x)则肯定接受它；否则，接受概率是f(y)f(x)\frac{f(y)}{f(x)}f(x)f(y)。这里需要指出：MCMC算法只需知道f的相对值，即只要给出一个正比于f的函数即可，这种方便也是MCMC的优势之一，因为有些应用问题难以将f归一化。虽然这个Metropolis-Hastings算法看起来简单，但它却非常有用，连同所有改进的算法一起，它们已在许多学科领域里有着重要的应用。

下面给出Metropolis-Hastings算法的伪代码：

如果建议概率密度函数g满足对称性条件，则概率函数h简化为h(x,y)=min(1,f(y)f(x))h(x, y)=min(1,\frac{f(y)}{f(x)})h(x,y)=min(1,f(x)f(y))，此时算法被称为对称Metropolis-Hastings算法

MCMC算法举例

模拟掷双骰子的游戏，它的状态空间为Ω=2,3,⋯,12\Omega ={2,3,\cdots,12}Ω=2,3,⋯,12下面的表给出了未归一的目标概率密度函数f，为了比较采用不同的建议概率密度函数g对于MCMC的影响，我们用两种不同的方法来构造建议概率密度函数：极大邻域法和极小邻域法

x	2	3	4	5	6	7	8	9	10	11	12
f(x)	1	2	3	4	5	6	5	4	3	2	1

极小邻域法：这个方法为每个状态构造的邻域是由与它最相邻的状态组成，即邻域只有一个或者两个元素。具体说，对于每个状态x，定义建议概率密度函数为：
g(y∣x)={1/2,y=maxx−1,2,1/2,y=maxx+1,12,0,,othersg(y|x)=\left\{ \begin{array}{rcl} 1/2 &, & y=max{x-1,2}, & \\ 1/2 &, & y=max{x+1,12}, & \\ 0, &, & others \end{array} \right. g(y∣x)=⎩⎨⎧1/21/20,,,,y=maxx−1,2,y=maxx+1,12,others

注意到：g(3|2)=g(2|2)=1/2,g(11|12)=g(12|12)=1/2，这样的定义保证了对称性条件被满足，可以看到这样的建议矩阵是对称的，这样的状态转移方式明显是能够便利的，因为只要有足够多的迭代次数，它会在各状态之间不断游走，只有这样，MCMC采能够生成正确的数据，提供一个满足目标要求的不变分布。

close all
clear
f = [0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1];
d = zeros(1, 20000);
x = 5;for i = 1:20000u = rand;if x == 2if u < 0.5y = 3;else y = 2;endelseif x == 12if u < 0.5y = 11;elsey = 12;endelseif u < 0.5y = x - 1;elsey = x + 1;endendh = min(1, f(y) / f(x));u = rand;if u < hx = y;endd(i) = x;
enda = 1:1:12;
hist(d, a);

我们看到，从mcmc模拟结果几乎精确地再现了目标概率分布。

极大邻域法：对所有的x定义其邻域为Nx=ΩN_x = \OmegaNx=Ω，也就是说，对于任意状态x，任何其他状态都可能被建议。如果以等概率取所有可能状态，这样的建议矩阵如下：
G=[1/111/11⋯1/111/111/11⋯1/11⋮⋮⋮1/111/11⋯1/11]G = \left[ \begin{array}{l} 1/11 & 1/11 & \cdots &1/11\\ 1/11 & 1/11 & \cdots &1/11\\ \vdots& \vdots & &\vdots\\ 1/11& 1/11& \cdots& 1/11 \end{array} \right] G=⎣⎢⎢⎢⎡1/111/11⋮1/111/111/11⋮1/11⋯⋯⋯1/111/11⋮1/11⎦⎥⎥⎥⎤

这种邻域方式的建议概率密度函数g(.|x)是常量函数，它当然也满足对称性条件，这一，这种邻域选择方式将产生独立的抽样序列。

close all
clear
f = [0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1];
d = zeros(1, 20000);
x = 5;for i = 1:20000u = rand;if u < 1/11y = 2;elseif u < 2 / 11y = 3;elseif u < 3 / 11y = 4;elseif u < 4 / 11y = 5;elseif u < 5 / 11y = 6;elseif u < 6 / 11y = 7;elseif u < 7 / 11y = 8;elseif u < 8 / 11y = 9;elseif u < 9 / 11y = 10;elseif u < 10 / 11y = 11;elsey = 12;endh = min(1, f(y) / f(x));u = rand;if u < hx = y;endd(i) = x;
enda = 1:1:12;
hist(d, a);

我们看到，模拟结果同前面是一样的，这体现了MCMC的灵活性。

连续马尔科夫链的模特卡罗模拟
考虑密度函数为f=0.5x2e−xf=0.5x^2e^{-x}f=0.5x2e−x的伽马分布，我们要涉及的马尔科夫链是连续的状态空间Ω=[0,+∞)\Omega=[0,+\infty)Ω=[0,+∞)，这里要注意，目标概率密度函数有定义域的限制，也就是说，建议概率密度函数不能产生y<0y\lt 0y<0的状态；但同时要使得建议概率密度函数是对称的，即满足g(y∣x)=g(x∣y)g(y|x)=g(x|y)g(y∣x)=g(x∣y)，我们定义这样的建议概率密度函数：对所有的x，取半径为1，中心在x的均匀分布密度，但在任何状态不能小于0.如果y<0y\lt 0y<0，则令y=xy=xy=x，这样类似于上面离散例子的做法避免了不对称的情况

f = @(x)(0.5 * x * x * exp(-x));
d = zeros(1, 40000);
x = 2;
for i = 1:40000y = unifrnd(x-1, x+1);if y < 0y = x;endh = min(1, f(y) / f(x));u = rand;if u < hx = y;endd(i) = x;
enda = 0:0.2:20;
hist(d(20001:40000), a);

MCMC模拟练习

我们计算在举例1中的掷骰子问题在不同初始分布p0p_0p0下的概率分布ptp_tpt,例子中已经完成了对p0=(1,0,0)p_0=(1,0,0)p0=(1,0,0)情况的模拟，对同样的问题，做p0=(0,0,1)p_0=(0,0,1)p0=(0,0,1)情况下的模拟；对两种不同选择的情况下计算pt,t=1,2,⋯,10p_t,t=1,2,\cdots,10pt,t=1,2,⋯,10，再利用这些计算，计算更一般情况p0=(a,b,1−a−b)p_0=(a,b,1-a-b)p0=(a,b,1−a−b)其中a,b≥0,a+b<1a,b\ge 0, a+b\lt 1a,b≥0,a+b<1(事实上，矩阵乘法是线性的)。试随机模拟对于任何初始起点p0p_0p0，分布ptp_tpt收敛于不变分布。

p = [0.2, 0.3, 0.5];P = [0.3, 0.2, 0.5; 0.4, 0.2, 0.4; 0.4, 0.3, 0.3];n = 25000;p_final = p;for i = 1:np_final = p_final * P;endp_final

将上述p配置成[1, 0, 0], [0, 1, 0], [0, 0, 1], [0.2, 0.3, 0.5]等满足题意的初始分布，经过计算后得到的结果均为：
p_final =

0.3636    0.2397    0.3967

故该马尔科夫链的最终均收敛于平稳状态，与初始分布无关；

利用Metropolis算法从掷两个均匀流面骰子之和的样本空间Ω=2,3,⋯,12\Omega={2,3,\cdots,12}Ω=2,3,⋯,12取样，首先用Ω=2,3,⋯,12\Omega={2,3,\cdots,12}Ω=2,3,⋯,12上的均匀分布作为系统邻域，其次运用转移矩阵
G=[1/21/20⋯0001/201/2⋯00001/20⋯000⋮⋮⋮⋱⋮⋮⋮000⋯01/20000⋯1/201/2000⋯01/21/2]G=\left[ \begin{array}{l} 1/2& 1/2& 0& \cdots& 0& 0& 0\\ 1/2& 0& 1/2& \cdots& 0& 0& 0\\ 0& 1/2& 0& \cdots& 0& 0& 0\\ \vdots& \vdots& \vdots& \ddots& \vdots& \vdots& \vdots \\ 0& 0& 0& \cdots& 0& 1/2& 0\\ 0& 0& 0& \cdots& 1/2& 0& 1/2\\ 0& 0& 0& \cdots& 0& 1/2& 1/2 \end{array} \right] G=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡1/21/20⋮0001/201/2⋮00001/20⋮000⋯⋯⋯⋱⋯⋯⋯000⋮01/20000⋮1/201/2000⋮01/21/2⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤

给出随机游走作为系统邻域。记录下计算时间，并画出直方图和相关图，发表评论你的意见。

仿真代码参考上面的掷骰子的最大邻域方法与最小邻域方法，经过计算，
最小邻域方法计算速度更快，最大邻域方法耗时较慢

继续考察上例的取样问题，探讨MCMC方案的收敛速度。首先采用由上体转移矩阵所给出的随即又走的系统邻域分布。设样本序列是X1,X2,⋯,XNX_1,X_2,\cdots,X_NX1,X2,⋯,XN，计算样本均值X‾=1N∑i=1NXi\overline{X}=\frac{1}{N}\sum_{i=1}^{N}X_iX=N1∑i=1NXi，运行N=100的100个十里，画出样本均值的直方图，并计算样本方差。再次重复运行N=200的100个实例，以及N=1000的100个实例。这给出了样本均值估计随迭代次数增加怎样收敛的启发。现在除了最大邻域系统外，所有状态都选择等可能，比较两个不同邻域系统在其收敛性分布直方图和方差图。

从理论层面分析X‾\overline{X}X样本均值的理想状态是这个分布的期望值，通过计算很容易得到其均值为7，随着N的增大，样本均值越来约接近其理论均值7，在实例维度，样本均值近似于一个正态分布，仿真代码如下：

close all
%clearN = 10000;
M = 10000;
f = [0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1];
d1 = zeros(M, N);
x = 5;
tic
for j = 1:Mfor i = 1:Nu = rand;if u < 1/11y = 2;elseif u < 2 / 11y = 3;elseif u < 3 / 11y = 4;elseif u < 4 / 11y = 5;elseif u < 5 / 11y = 6;elseif u < 6 / 11y = 7;elseif u < 7 / 11y = 8;elseif u < 8 / 11y = 9;elseif u < 9 / 11y = 10;elseif u < 10 / 11y = 11;elsey = 12;endh = min(1, f(y) / f(x));u = rand;if u < hx = y;endd1(j, i) = x;end
end
tocd_sum = sum(d1, 2) / N;
fprintf('sequence_length = %d, std = %f\n', N, std(d_sum));
[dNum, dCen] = hist(d_sum, 20);
bar(dCen, dNum);

>> ex3_2
Elapsed time is 0.109177 seconds.
sequence_length = 100, std = 0.280242
>> ex3_2
Elapsed time is 0.726260 seconds.
sequence_length = 1000, std = 0.088843
>> ex3_2
Elapsed time is 7.015829 seconds.
sequence_length = 10000, std = 0.028191

从输出可知，随着序列长度的增加，X‾\overline{X}X的方差不断减小

模拟马尔科夫链转移矩阵
P=[c012e−300120120012e−1c2120012e−1c3]P=\left[ \begin{array}{l} &c_0 &\frac{1}{2}e^{-3} &0 &0\\ &\frac{1}{2} &0 &\frac{1}{2} &0\\ &0 &\frac{1}{2}e^{-1} &c_2 &\frac{1}{2}\\ &0 &0 &\frac{1}{2}e^{-1} &c_3 \end{array} \right] P=⎣⎢⎢⎡c0210021e−3021e−10021c221e−10021c3⎦⎥⎥⎤

其中cic_ici是互补的概率，试模拟几百次迭代以获得不变分布并画出直方图，另外手工计算其不变分布，并比较计算一下结果：
ωi=e−fiZ\omega_i=\frac{e^{-f_i}}{Z} ωi=Ze−fi

其中Z=∑i=03e−fi,f0=−1,f1=2,f3=0Z=\sum_{i=0}^{3}e^{-f_i},f_0=-1,f_1=2,f_3=0Z=∑i=03e−fi,f0=−1,f1=2,f3=0

c0 = 1 - 0.5 * exp(-3);
c2 = 1 - 0.5 * exp(-1) - 0.5;
c3 = 1 - 0.5 * exp(-1);P = [c0, 0.5*exp(-3), 0, 0; 0.5, 0, 0.5, 0; 0, 0.5*exp(-1), c2, 0.5; 0, 0, 0.5*exp(-1), c3];p = [1, 0, 0, 0];N = 10000;p_final = zeros(1, 3);
p_final = p * P;for i = 2:Np_final = p_final * P;
endp_final

从1000次开始，分布进入不变状态，

>> ex5p_final =0.6439    0.0321    0.0871    0.2369>> ex5p_final =0.6439    0.0321    0.0871    0.2369>>

通过手工计算
ω0=0.0321,ω1=0.6439,ω2=0.2369,ω3=0.0871\omega_0 = 0.0321,\omega_1=0.6439,\omega_2=0.2369,\omega_3=0.0871 ω0=0.0321,ω1=0.6439,ω2=0.2369,ω3=0.0871

和通过模拟迭代得到的不变分布一致。

MCMC抽样---Metropolis算法相关推荐

如何实现马尔可夫链蒙特卡罗MCMC模型、Metropolis算法？
原文链接:http://tecdat.cn/?p=2687 什么是MCMC,什么时候使用它? MCMC只是一个从分布抽样的算法. 这只是众多算法之一.这个术语代表"马尔可夫链蒙特卡洛&quo ...
MCMC、蒙特卡洛近似和Metropolis算法简介
MCMC 是Markov Chain Monte Carlo 的简称,但在传统模拟中有一个很重要的假设是样本是独立的(independent samples),这一点在贝叶斯统计尤其是高纬度的模型中很 ...
蒙特卡洛算法_MCMC、蒙特卡洛近似和Metropolis算法简介
MCMC 是Markov Chain Monte Carlo 的简称,但在传统模拟中有一个很重要的假设是样本是独立的(independent samples),这一点在贝叶斯统计尤其是高纬度的模型中很 ...
python概率密度函数参数估计_Python与项目反应理论：基于EM和MCMC的参数估计算法...
项目反应理论的开端早在上世纪初,智力测验的发明者比奈(也可能是西蒙)便发现了一条神奇的曲线,这条曲线的x轴是智力水平,y轴是试题正确率,而这是项目反应理论(以下简称IRT)的最初雏形.上世界五六十年 ...
已知函数的分布，如何使用metropolis 算法去得到目标样本函数
Metropolis算法是一种随机搜索方法,可以用来从已知函数的分布中抽取样本.它利用了Markov链的思想,主要通过迭代多次抽样,不断改变抽样点的位置,从而逐步优化出最优样本. 首先,我们需要定义一 ...
随机游动Metropolis算法拟合标准拉普拉斯分布
1. 题目分析使用随机游动的Metropolis抽样方法产生标准拉普拉斯分布的随机数,其中使用正态分布产生增量,比较由提议分布不同方差生成的链的差异,并比较每个链的接收概率. 2. 代码展示我们首 ...
马尔可夫链蒙特卡洛（一）:Metropolis算法
0. 引言马尔可夫链蒙特卡洛(Markov Chain Monte Carlo,MCMC)是一类从复杂概率分布中进行采样的方法.Metropolis算法是严格意义上第一个MCMC算法,它于1953年 ...
Metropolis算法求解积分问题
计算机图形学技术见计算机图形学技术 Metropolis算法有两种计算积分时的产生随机数的方法: 已知目标概率分布 p(x) ,使用蒙特卡洛方法对该目标函数进行抽样. 已知转移概率分布,给定初始概 ...
悼念！蒙特卡洛Metropolis算法贡献者之一Arianna Rosenbluth逝世
点击上方"3D视觉工坊",选择"星标" 干货第一时间送达作者丨陈彩娴编辑丨AI科技评论 AI科技评论消息,洛杉矶当地时间12月28日,Metropolis算 ...

MCMC抽样---Metropolis算法

MCMC抽样---Metropolis算法相关推荐

最新文章

热门文章