开通VIP，畅享免费电子书等14项超值服
首页
好书
留言交流
下载APP
联系客服
2021.04.26
马尔可夫链蒙特卡洛方法(MarkovChainMonteCarlo)，简称MCMC。其产生于20世纪50年代早期，是在贝叶斯理论框架下，通过计算机进行模拟的蒙特卡洛方法(MonteCarlo)。该方法将马尔可夫(Markov)过程引入到MonteCarlo模拟中，实现抽样分布随模拟的进行而改变的动态模拟，弥补了传统的蒙特卡罗积分只能静态模拟的缺陷。
Metropolis等人在1953年首次提出了基于马尔可夫链的蒙特卡罗方法，即Metropolis算法，并在最早的计算机上编程实现。Metropolis算法是首个普适的采样方法，并启发了一系列MCMC方法，所以人们把它视为随机模拟技术腾飞的起点。Metropolis算法这篇论文[1]被收录在《统计学中的重大突破》中，《ComputinginScience&Engineering》尝试列出了对20世纪科学与工程的发展和实践影响最大的十种算法：
MetropolisAlgorithmforMonteCarlo被列为十大算法之首。用于蒙特卡洛的Metropolis算法定义了一个收敛的马尔可夫链，其极限就是所需的概率分布。Metropolis算法及其推广算法已被称为蒙特卡洛马尔可夫链技术(MCMC)，因为这些算法模拟了一个马尔可夫链，从极限分布中获取抽样。
搞自然科学研究的人，MCMC应该是基本装备。所谓平生不识MCMC,便称英雄也枉然。
MCMC由两个MC组成，即马尔可夫链(MarkovChain，也简称MC)和蒙特卡洛法（MonteCarloSimulation，简称MC）。要弄懂MCMC的原理我们首先得搞清楚马尔科夫链的原理和蒙特卡罗法。
马尔可夫链的命名来自俄国数学家安德雷.马尔可夫以纪念其首次提出马尔可夫链和对其收敛性质所做的研究。
马尔可夫链是一组具有马尔可夫性质的离散随机变量的集合。所谓马尔可夫性质是指某一时刻状态转移的概率只依赖于它的前一个状态。
具体地，马氏链的数学定义:
对概率空间内的离散随机变量集合
若随机变量的取值都在可数集S内：
且随机变量的条件概率满足如下关系
则随机变量集合X被称为马尔可夫链，可数集S被称为状态空间，马尔可夫链在状态空间内的取值称为状态。
马尔可夫过程通常分为3类：
我们先来看马氏链的一个具体的例子。社会学家经常把人按其经济状况分成3类：下层(lower-class)、中层(middle-class)、上层(upper-class)，我们用1、2、3分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是0.65，属于中层收入的概率是0.28，属于上层收入的概率是0.07。事实上，从父代到子代，收入阶层的变化的转移概率如图2-1和图2-2
图2-1
图2-2
设状态空间S={1,2,3}(1、2、3分别代表人的收入阶层：下层、中层、上层),使用矩阵的表示方式，转移概率矩阵记为:
假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量:
那么他们的子女的分布比例将是
他们的孙子代的分布比例将是
第n代子孙的收入分布比例将是
假设初始概率分布为π0=[0.21,0.68,0.11]，则我们可以计算前n代人的分布状况如下
我们发现从第7代人开始，这个分布就稳定不变了，这个是偶然的吗？我们换一个初始概率分布π0=[0.75,0.15,0.1]试试看，继续计算前n代人的分布状况如下:
我们发现，到第9代人的时候,分布又收敛了。最为奇特的是，两次给定不同的初始概率分布，最终都收敛到概率分布π=[0.286,0.489,0.225]，也就是说收敛的行为和初始概率分布π0无关。这说明这个收敛行为主要是由概率转移矩阵P决定的。我们计算一下P的n阶矩阵为
我们发现，当n足够大的时候，这个P的n阶矩阵的每一行都是稳定地收敛到π=[0.286,0.489,0.225]这个概率分布。这个收敛现象并非是我们这个马氏链独有的，而是绝大多数马氏链的共同行为。这个性质同样不光是离散状态，连续状态时也成立。
马氏链定理：如果一个非周期马氏链具有转移概率矩阵P,状态空间S,i,j∈S且它的任何两个状态是连通的，那么
这就是马氏链的收敛定理。
由该定理可以得出如下结论：
<1>
<2>
<3>π是方程πP=π的唯一非负解。其中,
π称为马氏链的平稳分布。
请注意两种表示方式的不同：
这个马氏链的收敛定理非常重要，所有的MCMC方法都是以这个定理作为理论基础的。定理的证明相对复杂，一般的随机过程课本中也不给证明，所以我们就不用纠结它的证明了，直接用这个定理的结论就好了。我们对这个定理的内容做一些解释说明：
我们仅证明定理的第二个结论
证明：
上式两边取极限就得到
证明完毕。
给定一个马尔可夫链，若在其状态空间S存在概率分布
则π是该马尔可夫链的平稳分布。
细致平衡条件（DetailedBalanceCondition）：给定一个马尔科夫链，概率分布π和概率转移矩阵P，i，j∈S,(S是状态空间)，如果下面等式成立:
则此马尔科夫链具有一个平稳分布（StationaryDistribution)。
需要注意，这个细致平衡条件仅是一个充分条件，而不是必要条件，也就是说存在具有平稳分布的马尔科夫链不一定满足此细致平衡条件。
证明如下：
由条件推出
得
所以此马尔科夫链具有一个平稳分布。
从初始概率分布π0(x)出发，我们在马氏链上做状态转移，记Xi的概率分布为πi(x)，则有
由马氏链收敛的定理,概率分布πi(x)将收敛到平稳分布π(x)。假设到第n步的时候马氏链收敛，则有
所以
也就是说无论我们的起始状态是何种概率分布，在状态转移矩阵P的作用下，经过足够大的n步转移之后，最后都收敛到一个平稳分布，且该分布是唯一不变的。即平稳分布仅由状态转移矩阵P决定。
设马尔可夫链的状态空间S={ω1，ω2，...,ωi，...ωn}（可以理解为样本空间在随机过程中的表达），如果我们从一个具体的初始状态x0∈S（x0可能是ω2，也可能是ωi）开始,沿着马氏链按照概率转移矩阵做跳转，那么我们可以得到一个转移状态序列，马氏链终于露出了它的庐山真面：
尽管马尔可夫链最终会收敛到所需的分布，但是初始样本可能会遵循非常不同的分布，尤其是在起点位于低密度区域的情况下。结果，通常需要一个老化期，所以我们会舍弃前面的老化样本，得到我们需要的平稳样本集
这就是平稳分布π的状态样本序列。最后平稳分布的样本集可能达到几万、几十万甚至更大的规模。统计各个状态出现的频次，归一化后画出统计直方图，就可以模拟目标分布。
1>输入马尔可夫链状态转移矩阵P,状态空间S，设定状态转移收敛次数为n1,迭代次数为n2;
2>从任意简单概率分布采样得到初始状态值x0∈S;
3>fort=0ton1+n2-1:从条件概率分布p(x|x(t))中采样，得到样本x=x(t+1).去掉前面的老化样本，最终的样本集
即为我们需要的平稳分布对应的样本集。
这里大家会问：
简单概率分布可选择正态分布或均值分布，正态分布根据Box–Muller变换法或接受-拒绝采样法采样，均值分布可直接计算机采样；条件概率分布可由输入转移矩阵P得到概率数据，然后根据离散随机变量逆变换法采样。这些采样方法后面有详细介绍。
两种采样方法来理解：
第1种方法，从简单概率分布中采样10000个的初始状态样本，用最简单粗暴的方法对每一个初始样本都执行一次n(n>n1)步转移，状态之间的转移概率从该马尔科夫链的状态转移矩阵获取，最后依如下概率落入状态空间中的一个状态。
当初始状态样本足够多时，由大数定律，最终落入到平稳马氏链的各个状态的数量所占的比例逼近各个状态的平稳分布概率。但是这样一来计算量相当大，每次只利用了采样过程中的最后一个状态，且要运行10000次进入平稳状态的转移过程，数据的利用率很低，浪费也很大。
第2种方法，只需要对一个初始状态进行状态转移，就可以实现整个采样过程。我们把它进入稳态的那一刻记作t0，当t0时刻的状态向下一个时刻t1转移时，它依照稳态分布概率
（由n-步转移概率的性质Chapman–Kolmogorov等式可知，n-步转移矩阵是其之前所有转移矩阵的连续矩阵乘法，证明见附录1）
图2-3示意了进入平稳期的采样全过程（假设总共3个状态）：
图2-3
马尔可夫链的4个性质：不可约性、常返性、周期性和遍历性。
马尔可夫链在处理随机动力学时对问题建模的强大作用。
马尔可夫链预测理论在教育学、经济学、金融投资、生物学、农作物栽培、地质灾变，...特别是水资源科学中都得到了极为广泛地应用。
例如，马尔可夫链在排队理论、统计数据、生物种群进化建模、隐马尔可夫模型信息理论、语音识别、音字转换、词性标注、蒙特卡洛采样等都有应用。
蒙特卡洛(MonteCarlo，MC)方法,又称随机抽样法，统计试验法或随机模拟法。
蒙特卡罗法是由冯·诺伊曼等人在20世纪40年代为研制核武器的需要而首先提出并以摩纳哥的赌城MonteCarlo来命名的一种计算方法。
蒙特卡洛法是数学建模十大算法之一，是思想和技巧的艺术品。它的魔力在于，对于那些规模极大的问题，求解难度随着问题的维数（自变量个数）的增加呈指数级别增长、出现所谓的“维数的灾难”的问题，传统方法无能为力，而蒙特卡洛方法却可以独辟蹊径，基于随机仿真的过程给出近似的结果。
蒙特卡洛法是一类通过随机抽样过程来求取最优解的方法的总称，如果建立蒙特卡洛模型的过程没有出错，那么抽样次数越多，得到的答案越精确。蒙特卡洛法可以归纳为如下三个步骤：
对于本身就具有随机性质的问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。
实现从已知概率分布抽样。
构造了概率模型以后，按照这个概率分布抽取随机变量。这一般可以直接由软件工具包（见附录2，3，4）调用，或抽取均匀分布的随机数构造。这成为实现蒙特卡洛方法模拟实验的基本手段，也是蒙特卡洛方法被称为随机抽样的原因。
一般说来，构造了概率模型并能从中抽样后（即实现模拟实验后），我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。
蒙特卡洛方法应用非常广泛，其领域包括统计物理、天体物理、物理化学、数学、计算生物、统计学、经济学、金融证券、社会科学等等。
蒙特卡洛方法在解决物理和数学问题时时常被使用。当十分困难或不可能得到解析表达式，或不可能应用确定性算法时，该方法显得尤其重要和最为有用。
蒙特卡洛方法主要用于两个不同类别的问题：
无理数、不规则图形的面积、多重积分、求逆矩阵、解线性代数方程组、解积分方程、解某些偏微分方程边值问题和计算代数方程组、计算微分算子的特征值等等。
仿真、概率分布的生成（MCMC算法）。
解题思路：
构造一个单位正方形和一个单位圆的1/4，往整个区域内随机投入点，根据点到原点的距离判断点是落在1/4的圆内还是在圆外，从而根据落在两个不同区域的点的数目，求出两个区域的比值。如此一来，就可以求出1/4单位圆的面积，从而求出圆周率π。
解题步骤：
第1步，建模，构造概率过程，选择区域内随机投入点作为随机变量；
第2步，选择二维均匀分布Y=sqrt(h^2+v^2),h、v~U[0,1],作为已知概率分布抽样；
第3步，根据落入正方形样本数和1/4圆样本数，求出结果。
演示代码：
例子1
解题思路如下：
假定随机变量具有密度函数
根据数学期望公式得
构造统计数
根据大数定律
最终求出了定积分。
例子2求定积分
将原式变形为
解：采用蒙特卡洛方法，MATLAB模拟，
N=500;
x=unifrnd(0,2,N,1);
mean(2*exp(3*x.^2))
第一条语句设定了停止条件，共做N次MonteCarlo模拟；
第二条语句实现了在积分区间[0,2]均匀产生N个随机数；
第三条语句求样本均值，实现蒙特卡洛计算方法的面积逼近。
图4-1是一个中子穿过用于中子屏蔽的铅墙示意图。铅墙的高度远大于左右厚度。设中子是垂直由左端进入铅墙，在铅墙中运行一个单位距离然后与一个铅原子碰撞。碰撞后，任意改变方向，并继续运行一个单位后与另一个铅原子碰撞。这样下去，如果中子在铅墙里消耗掉所有的能量或者从左端逸出就被视为中子被铅墙挡住，如果中子穿过铅墙由右端逸出就视为中子逸出。如果铅墙厚度为5个单位，中子运行7个单位后能量耗尽，求中子逸出的几率。
这个问题并不复杂，但不容易找到一个解析表达式。而用模拟仿真的方法求解却可以有满意的结果。解题步骤如下：
第1步，建模，构造概率过程，选择中子碰撞后的角度作为随机变量；
第2步，选择角度的均匀分布作为已知概率分布抽样；
第3步，根据假设的总中子样本数和逸出的中子样本数，求出结果。
解：具体做法如下：
我们关心的是一次碰撞后，中子在水平方向行进了多少，所以行进方向是正负θ的结果是一样的，我们就只考虑θ是正的情形。由于中子运行的方向θ是随机的，我们用计算机抽取在0到π间均衡分布的随机数，模拟1000000个中子在铅墙里行进的情形，看看这些中子与铅原子碰撞7次后，有多少超过了铅墙的右端。
实现伪代码
n=1000000;m=0;flag=1;fori=1:nx=1;forj=1:7angle=pi*rand;x=x+cos(angle);ifx<0k=0;flag=0;endendifx>5&flag==1k=1;elsek=0;endm=m+k;flag=1;end运行结果：m/n=1.3%
我们运行程序得出逸出铅墙的中子的可能性约为1.3%，有了这个数字，我们可以报告安全部门，如果数字不能达到安全要求，我们则要加厚铅墙。
对于给定的概率分布p(x)，我们希望能有便捷的方式生成它对应的样本。由于马氏链能收敛到平稳分布，于是一个很的漂亮想法是：如果我们能构造一个转移矩阵为P的马氏链，使得该马氏链的平稳分布恰好是p(x),那么我们从任何一个初始状态x0出发沿着马氏链转移,得到一个转移状态序列
如果马氏链在第n步已经收敛了，于是我们就得到了平稳分布p(x)的样本集
这个绝妙的想法在1953年被Metropolis等人想到了，为了研究粒子系统的平稳性质，Metropolis考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法，即Metropolis算法，并在最早的计算机上编程实现。
从贝叶斯(Bayesian)的观点来看，模型中的观测变量和参数都是随机变量，因此，样本x与参数θ的联合分布可以表示为：
根据贝叶斯理论，可以通过样本x的信息对θ的分布的进行更新，后验概率为：
其中π(θ)先验概率；f(x|θ)似然函数（即样本的概率分布）；分母是标准化常量。后验概率=(似然函数*先验概率)/标准化常量=标准似然度*先验概率。
贝叶斯学习中经常需要进行三种积分运算：
1>后验概率计算中需要归范化计算：
2>如果有隐变量z∈Z，后验概率的计算需要边缘化计算：
3>如果有一个函数g(θ)，可以计算该函数的关于后验概率分布的数学期望：
此积分值为样本x的函数，因此可以对g(θ)进行推断。
该期望估计
当维数很高时，直接进行这三类积分是非常困难的，MCMC方法就是为此目的而诞生，为这些计算提供了一个通用的有效解决方案。
对很多贝叶斯推断问题来说，有时候后验分布过于复杂，使得积分没有显示结果，数值方法也很难应用；有时候需要计算多重积分（比如后验分布是多元分布时）。这些都会带来计算上的很大困难。这也是在很长的时期内，贝叶斯统计得不到快速发展的一个原因。1990年代MCMC计算方法引入到贝叶斯统计学之后，一举解决了这个计算的难题。MCMC方法的最新发展使计算大型分层模型成为可能，这些模型需要对数百到数千个未知参数进行积分。可以说，近年来贝叶斯统计的蓬勃发展，特别是在各个学科的广泛应用和MCMC方法的使用有着极其密切的关系。
MCMC算法是一个普适的采样方法，它的出现带来了随机模拟技术的腾飞。在随机变量生成技术中，MCMC是一种相当高级的方法，可以从一个非常困难的概率分布中获得样本。更出乎意料的是，可以用MCMC从一个未经标准化的分布中获得样本，这来自于定义马尔可夫链的特定方式，马尔可夫链对这些归一化因子并不敏感。
对高维（多个自变量）积分的数值近似计算是MCMC方法的一种重要的应用，很多统计问题，比如计算概率、矩（期望、方差等）都可以归结为定积分的计算。
MCMC的精髓在于构造合适的马氏链。
MCMC是一种简单有效的计算方法，在很多领域到广泛的应用，如统计学、贝叶斯(Bayes)问题、计算机问题等。
MCMC背后的基本思想是构造一个遍历的马尔可夫链，使得其平稳的、收敛的极限分布成为人们所需要的目标概率分布。
设马尔可夫链的状态空间S={ω1，ω2，...ωn}（可以理解为样本空间在随机过程中的表达），从某个初始状态x0出发，用构造的马尔可夫链随机游走，产生状态样本序列:
舍弃前面的老化(不稳定)的样本集合
应用遍历定理，确定正整数m和n,得到平稳分布的状态样本集合：
该样本集可能达到几万、几十万、几千万甚至更大的规模。统计该状态样本集中对应的各个状态ω1，ω2...出现的频次，归一化后画成直方图。如图3-1所示，绿色部分就是状态概率分布（即抽样概率分布），红色部分是目标概率分布。t趋近∞时，抽样概率分布可以无限逼近目标概率分布。如图4-1所示。
图4-1
平时我们接触比较多的场景是，给定一堆样本数据，求出这堆样本的概率分布p(x)。而采样刚好是个逆命题：给定一个概率分布p(x)，如何生成满足条件的样本？
通过样本可以计算各类统计量（如均值、方差、矩等），而此类统计量就是欲求问题的近似解。例如求积分，根据大数定律，相当于求某概率分布采样样本的均值。
蒙特卡罗模拟有一个危险的缺陷：如果必须输入一个模式中的随机数并不像设想的那样是随机数，而却构成一些微妙的非随机模式，那么整个的模拟和预测结果都可能是错的。
对于给定的概率分布来说，采样必须是随机的，这个非常重要。随机采样对于人来说，是不能也，非不为也。依靠人来采样很难得到真正的随机样本，人为地通过概率分布函数计算得到一批样本注定是徒劳的。
真随机数生成器（TrueRandomNumberGenerator,TRNG）是一种通过物理过程而不是计算机程序来生成随机数字的设备。通常是基于一些能生成随机的“噪声”信号的微观现象，如热力学噪声、光电效应和量子现象等，其物理过程在理论上是完全不可预测的，并且这种不可预测性要经过实验检验。真随机数生成器通常由换能器、放大器和模拟数字转换器组成。其中换能器用来将物理过程中的某些效果转换为电信号，放大器及其电路用来将随机扰动的振幅放大到宏观级别，而模拟数字转换器则用来将输出变成数字，通常是二进制的零和一。通过重复采样这些随机的信号，一系列的随机数得以生成。
硬件随机数发生器通常每秒仅产生有限数量的随机位。为了增加可用的输出数据速率，通常将它们用于生成“种子”，以用于更快的密码安全伪随机数生成器，然后生成器以更高的数据速率生成伪随机输出序列。
伪随机数发生器是依靠计算机按照一定算法模拟产生的，其结果是确定的，是可见的。通过这种方式产生的“随机数”并不随机，是伪随机数。贝尔实验室的里德博士告诫人们记住伟大的诺依曼的忠告:“任何人如果相信计算机能够产生出真正的随机的数序组都是疯子。”
伪随机数存在两大问题
真随机数的生成对技术的要求比较高，在实际应用中往往使用伪随机数就够了。在使用伪随机数时，要尽量克服其存在的两大问题。对于第一个问题，不能从本质上加以改变，但只要递推公式选的比较好，随机数间的相互独立性是可以近似满足的；对于第二个问题，因为用蒙特卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所用随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。
我们通过采样获取足够多的样本，然后统计样本空间中的样本出现的频次，来模拟和逼近目标概率分布的。采样样本越多，越逼近目标概率分布。如图5-1所示。
几万或几十万的样本集是小目标；几百万或几千万的样本集很常见；在AI和大数据时代，目标概率分布有时会涉及到几百个、几千个参数，几百亿的样本集都出现了。
图5-1
这与算法和采样策略有关。
这与算法有关，要从理论上证明。
例如：线性同余伪随机数生成器
X是伪随机序列；
m,m>0表示模量；
a,0 c,0≤c X0，0≤X0 如果c=0,称为乘法同余发生器，如果c!=0,称为混合同余发生器； Yn则是[0,1)内服从均匀分布的随机变量。 均匀分布随机数可以直接由计算机模拟产生。 Box–Muller变换是一种利用均匀分布快速产生符合标准正态分布伪随机数对的一种方法。基本思想是先得到服从均匀分布的随机数，再将服从均匀分布的随机数转变为服从正态分布。 假设变变量U1和变量U2是（0，1]均匀分布的随机数， 且U1和U2彼此独立，令： 则Z0和Z1就是服从N（0，1）的标准正态分布随机数，并且Z0和Z1相互独立。 假设变量u和变量v是在[-1，1]上的均匀分布随机量， u和v相互独立，令： 故而，随机数z0和z1计算后得出如下结果： z0和z1是服从分布N（0，1）的随机数，并且z0和z1相互独立。 逆变换法就是利用累积分布函数求逆，将所求概率分布的抽样转化为随机均匀分布的抽样。 设所求连续分布采样的连续随机变量X的概率密度函数为f(x)。 第1步，求概率累积分布函数F(x) F(x)是单调递增的连续函数，令 则U也是随机变量，且U∈[0,1]。 第2步，求反函数 第3步，随机产生均匀分布U[0,1]，求出对应的X，就可得到所求概率分布的采样。 设所求离散分布采样的离散随机变量的概率密度为 第1步，生成均匀分布U[0,1]随机数。 第2步，求出对应的X 也就是如下分段式： X就是所求概率分布的采样。 任意离散分布都可以画在轮盘上，采样只需要随机地旋转轮盘即可。如图5-2所示。 该轮盘模拟了累积分布函数，随机转动轮盘相当于随机均匀采样，最后轮盘的停止位置就是采样值。 图5-2 优点：随机数可以直接由均匀分布导出。 缺点：很多复杂累积分布函数求逆困难。 实际随机分布抽样中会遇到如下问题： 接受-拒绝采样法就派上了用场。基本思想是选择一个容易抽样的随机分布函数q(x)和常数k,使得k*q(x)≥p(x)（这里p(x)可以进行非归一化，乘以常数Z，以下省略），然后按照一定的方法拒绝某些抽样，达到接近p(x)分布的目的。 如图5-3所示： 图5-3 假设我们选取已知的高斯分布q(x)=Gassian(1.4,1.2)，我们按照一定的方法拒绝某些样本，达到接近p(x)分布的目的。 具体步骤如下： 第1步，选择已知的分布函数q(x)，确定常量k，使得p(x)总在kq(x)的下方。 第2步，在x轴上，从q(x)分布抽样取得a。但是a不一定留下，会有一定的几率被拒绝。 第3步，在y轴上，从均匀分布(0,1)中抽样得到u。如果u.kq(a)>p(a)，就是落到了灰色的区域中，拒绝；否则接受这次抽样。 第4步，重复2、3步。 模拟结果如图5-4所示： 图5-4 在高维（多个自变量）分布的情况下，接受-拒绝采样会出现如下一些问题： 1>合适的提议分布q(x)比较难以找到，并且当kq(x)与p(x)相差太多时，会导致拒绝率很高，无用计算增加。 2>很难确定一个合理的常数k值,并且k值往往会很大； 3>所需的样本量随空间维数增加而指数增长。 重要性采样是统计学中估计某一分布性质时使用的一种方法，该方法从与原分布不同的新分布中采样，而对原分布的性质进行估计。重要性采样常用于蒙特卡洛积分。 重要性采样是蒙特卡洛方法中的一个重要策略。该方法不改变统计量，只改变概率分布。重要性采样算法就是在有限的采样次数内，尽量让采样点覆盖对积分贡献很大的点或者以一种受控的方式改变仿真，以便增加稀少事件的数目。它是一种降方差的仿真方法。 其中，q(x)为提议分布(proposaldistribution)，p(x)/q(x)可以看做importanceweight。 优点：如果我们无法从分布p(x)中抽样,或者从中抽样的成本很高,那么我们还可以求E[f(x)]吗？答案是可以的,我们可以从一个简单的分布q(x)中进行抽样得到x,求f(x)p(x)/q(x)的期望，进而可以求得函数积分的近似解。 缺点：在高维空间里找到一个这样合适的q(x)非常难。 拒绝采样法和重要性采样法面临两个问题： 1>如果提议分布q(x)与目标分布p(x)差距过大，模拟采样效果就不太好； 2>构造一个与目标分布p(x)相似的提议分布q(x)很困难。 我们知道，对于任意的目标分布，我们可以设计一种接受-拒绝的概率，从而使得所有接受的点都是该分布的样本。拒绝采样法和重要性采样法采用了提议函数q(x)和接受概率实现了这种思想，然而问题在于，由于其提议分布q(x)是固定的，这对q(x)的要求很高，试想下，如果q与真实p的差距过大，那么几乎每个样本都会被拒绝掉，效率很低。那如何设计更好的q呢？一个巧妙的办法是设计一个自适应的提议分布q(x’|x)，它以上一次接受的样本作为条件，然后经过转移propose一个新的样本，如果我们能够约束这个转移不会改变目标分布，那么我们就能快速的找到目标分布p(x)的样本。如图5-5所示，这也正是MCMC的思想。 图5-5 另一方面，大多数拒绝采样法和重要性采样法都遭受“维数的诅咒”，其中拒绝的概率随维数的增加呈指数增长。MCMC方法都没有这个“维数诅咒”的问题，当要采样的分布维数很大时，它通常是唯一可用的解决方案。结果，MCMC方法成为了从分层贝叶斯模型和当今许多学科中使用的其它高维统计模型中生成样本的首选方法。MCMC采样是一种普适采样方法，它可以解决其它采样方法无法搞定的很多问题。 提议分布q(x,x’)=q(x’|x)通常要满足以下三个条件： 1>对于固定的x,q(.)是一个概率密度函数; 2>对于x,x’∈S,q(x,x’)的值要能计算; 3>对于固定的x,能够方便地从q(x,x’)中产生随机数。 从理论上讲，提议函数q(.)的选取是任意的，但在实际计算中，提议函数的选取对于算法效率的影响是相当大的。一般认为提议函数与目标分布越接近，模拟的效果越好。 q(.)函数的选择一般有两种： 1>对称:如高斯分布，均匀分布，这种情况下称为随机游走; 2>非对称：如log-normal，采样值倾向于选择高概率的点，因为我们不希望采样点来回震荡。 MCMC采样法的关键是构建马氏链，马氏链的收敛性质主要由状态转移矩阵P决定,所以基于马氏链做采样的关键问题是如何构造状态转移矩阵P,使得马氏链具有一个平稳分布，且该平稳分布π(x)恰好是我们要的分布p(x)。如何能做到这一点呢？ 答案的关键在细致平衡条件 其中i，j∈S(S是状态空间)，P是状态转移矩阵。 不幸的是，对于平稳分布π(x)，即我们要的分布p(x)，要找到满足细致平衡条件的状态矩阵P很困难。没关系，我们随机找一个已知的马尔科夫链状态转移矩阵Q，其中Q(i,j)=Q(j|i)是从状态i转移到状态j的提议分布条件概率，i,j∈状态空间S，它可能一开始并不满足细致平衡条件，即： 我们可以对上式做一个改造，使细致平衡条件成立。方法是引入一个α(i,j),使上式可以取等号，即： 要使该等式成立，只要满足下面两式即可 这样，我们就得到了平稳分布π(x)对应的马尔科夫链状态转移矩阵P 我们有一般称α(i,j)为接受概率 其取值在(0,1)之间，可以理解为一个概率值。 也就是说，基于平稳分布的状态转移矩阵P可以通过任意一个马尔科夫链状态转移矩阵Q以一定的接受概率获得。在接受-拒绝采样法中，那里是以一个常见分布通过一定的接受-拒绝概率得到目标概率分布，这里是以一个常见的马尔科夫链状态转移矩阵Q通过一定的接受-拒绝概率得到目标转移矩阵P，两者解决问题的思路是类似的。 因为α(i,j)<1和α(j,i)<1,我们可以将两者同时放大，且不破环等式成立。 1>当α(i,j)≤α(j,i)时，则α(i,j)按照α(j,i)放大到1（概率不可能大于1）的比例等比例放大: 2>当α(i,j)>α(j,i)时，则α(j,i)按照α(i,j)放大到1的比例等比例放大: 所以将两者同时放大有 此时α(i,j)就是Metropolis-Hasting算法的接受概率。其中π(x)=p(x),Q(i,j)=Q(j|i)是提议分布条件概率，α(i,j)在取到1的情况下能实现链的满概率跳转，否则以放大1/(π(i)Q(i,j))倍的概率进行跳转。 当提议分布对称时，即Q(j,i)=Q(i,j)，接受概率 此时α(i,j)就是Metropolis算法的接受概率。其中π(x)=p(x)。 由任选的提议分布Q(i,j)和接受概率α(j,i),可以构造一个满足细致平衡条件的马氏链，该马氏链具有一个平稳分布。构造平稳分布的马氏链的过程就是我们实现MCMC采样的过程。 MCMC算法包括Metropolis算法、Metropolis-Hasting采样法、吉布斯采样法等。 已知目标概率分布p(x)，如何对其进行采用？ 我们分别采用接受-拒绝算法和Metropolis算法进行采样，并对两者做了对比。 接受-拒绝算法的具体步骤如下： 1>首先选择一个容易抽样的分布作为提议分布，记为q(x)，接着确定一个常数C，其中C>1，使得在x的定义域上都有p(x)≤Cq(x)。 2>生成服从提议概率密度函数q(x)的随机数y。 3>接着生成一个服从均匀分布U(0,1)的随机数u。 4>计算接受概率 如果u<α(y)，则接受y;否则，拒绝。 5>不断重复2~4就可以生成一组服从目标分布p(x)的随机数序列。 Metropolis算法就是改进的接受-拒绝算法，具体步骤如下： 1>首先构造提议分布q(.|x(t)),这里提议分布是对称的，即q(y|x)=q(x|y)，例如，我们可以选均值为u=x(t)的正态分布为提议分布。设置t=0和初始状态x0。 2>从提议概率分布q(.|x(t))中产生一个随机数y作为下一个状态。 如果u≤α，则接受候选随机数y,赋值x(t+1)=y，y作为提议分布的新状态，意味着此操作将鼓励跳转到目标分布的高密度区。 如果u>α，则拒绝候选随机数，赋值x(t+1)=x(t)，表明f(x)在此y状态为低密区,以概率的形式被拒绝。 5>增量设置t=t+1。 6>重复步骤2~步骤5，我们就可以得到一组仅依赖于上一个随机数并且与前面随机数无关的随机序列 Metropolis算法与接受-拒绝算法不同之处在于： 1>Metropolis算法的提议分布是自适应的条件概率分布，产生的随机数是由上一次的状态值转移得到；接受-拒绝算法的提议分布是固定的，没有状态转移的概念。 2>两者的接受概率定义不同。 Metropolis算法中唯一的限制是提议分布必须是对称的，满足条件的这类分布有正态分布、柯西分布、t分布和均匀分布等。为了能够使用非对称分布，或者目标分布的支持集是不对称的，如x∈(0,+∞),我们将考虑使用Metropolis-Hastings算法。 例子已知：(贝叶斯推断:一个简单的投资模型)假设有5种股票被跟踪记录了250的交易日每天的表现,在每一个交易日,收益最大的股票被标记出来。用Xi表示股票i在250个交易日中胜出的天数,则记录到得的频数(x1,x2,...,x5)为随机变量(X1,X2,...,X5)的观测。基于历史数据,假设这5种股票在任何给定的一个交易日能胜出的先验机会比率为1:(1β):(12β):2β:β,这里β∈(0,0.5)是一个未知的参数。在有了当前这250个交易日的数据后,使用Bayes方法对此比例进行更新。 求：根据观测的数据求β目标分布随机抽样的样本集，进而估计β？ 解：根据前述，(X1,X2,...,X5)在给定β的条件下服从多项分布,概率向量为 因此β服从后验分布为（根据联合密度的多项分布公式，把x1,x2,...,x5看作常数） 我们不能直接从此后验分布中产生随机数。一种估计β的方法是产生一个链,使其平稳分布为此后验分布,然后从产生的链中抽样来估计β. 我们这里使用随机游动Metropolis算法,可以初始化β=0.2,提议分布为均匀分布U[0,1],收敛次数burn=1000,迭代次数m=5000,接受的概率为 因为后验分布即为所求的目标分布f(x),可得 观测数据(x1,x2,x3,x4,x5)是5只股票服从概率向量p在250天中胜出的频次。 Metropolis算法的应用受到了提议概率分布的对称形式的限制，因此，Hastings把提议概率分布从对称分布形式推广到一般情况，形成了Metropolis-Hastings算法，简称M-H算法，也称之为通用梅特罗波利斯算法。 1>首先构造提议概率分布q(.|x(t))。设置t=0和初始状态x0。 2>从提议分布q(.|x(t))中产生一个随机数y作为下一个状态。 如果u≤α，则接受候选随机数y,赋值x(t+1)=y; 如果u>α，则拒绝候选随机数，赋值x(t+1)=x(t)。 6>重复步骤2~步骤5，我们就可以得到一组仅依赖于上一个随机数并且与前面随机数无关的随机序列。 第1次采样，选择初始状态x0 图5-6a 第2次采样，接受候选样本x1 第3次采样，接受候选样本x2 图5-6c 第4次采样，拒绝候选样本x3 图5-6d 第5次采样，接受候选样本x4 图5-6e 马尔可夫链蒙特卡罗法的目标分布通常是多元联合概率分布 如果条件概率分布 中所有k个变量全部出现，其中 那么称这种条件概率分布为满条件分布（fullconditionaldistribution)。 满条件分布有以下性质： 1>对任意的x,x’和任意的IK,有 2>对任意的x,x’和任意的IK,有 M-H算法中，可以利用上述性质，简化计算，提高计算效率。具体地，通过满条件分布概率的比 代替计算联合概率的比 而前者更容易计算。 在Metropolis-Hastings算法中，通常需要对多元变量分布进行抽样，有时对多元变量分布的抽样是困难的。 可以对多元变量的每一变量的条件分布依次分别进行抽样，从而实现对整个多元变量的一次抽样，这就是单分量Metropolis-Hastings(single-componentMetropolis-Hastings）算法。 假设马尔可夫链的状态由k维随机变量表示 马尔可夫链在时刻i的状态 为了生成容量为的样本集合 单分量Metropolis-Hastings算法的基本做法：假设有一个k维的随机向量，现想要构造一条k维向量n个样本的马尔可夫链序列，那么随机初始化一个k维向量，然后固定这个向量其中的k-1个元素，抽取剩下的那个元素，生成转移状态的随机数，这样循环k次，就把整个向量更新了一遍，也就是生成了一个新的k维向量样本，把这个整体重复n次就得到了一条马尔科夫链。 单分量M-H算法由下面的k步迭代实现M-H算法的一次迭代 马尔可夫链的转移概率为 如图5-7所示，单分量M-H算法的迭代过程。目标是对含有两个变量的随机变量x进行抽样 如果变量x1或x2更新，那么在水平或垂直方向产生一个移动，连续水平和垂直移动产生 一个新的样本点。 图5-7 例子逻辑回归应用 考虑54位老年人的智力测试成绩(WecheslerAdultIntelligenceScale,WAIS,0-20分)。研究的兴趣在于发现老年痴呆症。 我们采用如下简单的逻辑回归模型: 则似然函数为 考虑β0,β1的先验分布π为独立的正态分布: 从而我们需要从此分布中产生随机数。 我们考虑如下三种抽样方法. 提议分布取为 则算法如下： cor(beta[(burnin+1):m,1],beta[(burnin+1):m,2])=-0.954. 由似然函数,容易计算得到Fisher信息阵为 算法如下： 在计算exp(η)时,由于在R中,exp(710)=Inf,因此我们对η>700,取近似log(1+exp(η))≈η. 对链的收敛诊断可以看出,链混合的效率很高。 在M-H算法中,按分量进行逐个更新,其优势在于应用方便,不需要考虑调节参数.算法如下: y<-wais[,2]x<-wais[,1]m<-10000beta0<-c(0,0)#initialvaluemu.beta<-c(0,0)#priors.beta<-c(100,100)#priorprop.s<-c(1.75,0.2)#sdofproposalnormalbeta<-matrix(nrow=m,ncol=2)acc.prob<-c(0,0)current.beta<-beta0for(tin1:m){for(jin1:2){prop.beta<-current.betaprop.beta[j]<-rnorm(1,current.beta[j],prop.s[j])cur.eta<-current.beta[1]+current.beta[2]*xprop.eta<-prop.beta[1]+prop.beta[2]*xif(sum(prop.eta>700)>0){print(t);stop;}if(sum(cur.eta>700)>0){print(t);stop;}loga<-(sum(y*prop.eta-log(1+exp(prop.eta)))-sum(y*cur.eta-log(1+exp(cur.eta)))+sum(dnorm(prop.beta,mu.beta,s.beta,log=TRUE))-sum(dnorm(current.beta,mu.beta,s.beta,log=TRUE)))u<-runif(1)u<-log(u)if(u 吉布斯采样算法是StuartGeman和DonaldGeman两兄弟在1984年提出来的，以物理学家JosiahWillardGibbs的名字命名。这个算法在现代贝叶斯分析中占据重要位置。 吉布斯采样算法是单分量M-H算法的特殊情况，但是更容易实现，因而被广泛使用。 吉布斯采样特别适合于采样贝叶斯网络的后验分布，因为贝叶斯网络通常被指定为条件分布的集合。 1>吉布斯采样算法的提议分布是当前变量xj，j=1,2,…，k的满条件概率分布 2>吉布斯采样算法的接受概率α=1 一般的采样方法，生成一维随机变量并不困难，对于高维分布 生成各分量非独立的随机向量就非常困难。吉布斯采样法很擅长解决这类问题，其广泛用于多元变量联合分布的抽样和估计。 吉布斯采样法把一次一个n维样本变为“n次一维”样本，但前提是我们事先知道完全条件概率分布 数学上称为半共轭，为方便常简写为 不同于M-H算法，该算法不引入新的提议分布q(.)，仅通过条件分布得到以给定分布p(x)为平稳分布的马氏链的转移概率。 把条件分布当作提议分布q(.)，则有 接受概率 以接受概率为1的方式实现链的满概率跳转。 图5-8 假设有一个二维概率分布p(x,y),考察x坐标相同的两个点A(x1,y1),B(x1,y2)。如图5-8所示。 由条件概率公式推得 上面两式右边相等，所以 故 从上面等式，可以看出，在x=x1这条平行于y轴的直线上，如果使用条件分布p(y|x1)作为任意两点之间的转移概率，那么任意两点之间的转移满足细致平稳条件。同样地，如果在y=y1这条直线上任意取两个点A(x1,y1),C(x2,y1),也有如下等式： 于是，我们可以构造平面上任意两点之间的转移概率矩阵Q 根据上式，我们很容易验证对平面上任意两点X、Y(如图5-8的两点A、B),满足细致平稳条件 因此，这个二维空间上的马氏链将收敛到平稳分布p(x,y)。 马氏链的状态转移只是轮换地沿着x轴和y轴做转移，于是得到样本(x0,y0),(x0,y1),(x1,y1),(x1,y2),(x2,y2),...马氏链收敛后，最终得到的样本就是p(x,y)的样本。坐标轴轮换采样是不强制要求的，也可以在x轴和y轴之间随机地选一个坐标轴，然后按条件概率做转移，马氏链也是一样收敛的。 我们很容易把二维推广到高维的情形，对于二维细致平稳条件 如果x1变为多维X1，可以看出推导过程不变，所以多维细致平稳条件同样是成立的。 此时转移矩阵Q由条件分布p(y|X1)定义。上式只是说明了一个坐标轴的情况，和二维情形类似，很容易验证对所有的坐标轴都有类似的结论。 所以，n维空间对于概率分布p(x1,x2,...,xn)可以定义如下转移矩阵 1>如果当前状态为(x1,x2,...,xn)，马氏链状态转移的过程中，只能沿着坐标轴做转移，转移概率定义为 2>其它无法沿着坐标轴进行的跳转，转移概率都设置为0。 于是，我们可以把GibbsSampling算法从二维的概率分布p(x,y)推广到n维的概率分布p(x1,x2,...,xn)。 例子已知：贝叶斯后验分布有如下形式 求：利用Gibbs采样法获取满足该后验分布的样本序列 解： 由上面两式，求得条件概率分布 这里A(y),B(x)是关于y,x的函数，当指定y,x时，它们是常数。 求得的条件概率分布随机变量x,y都符合正态分布，即 按已知正态分布f(x|y,data)和f(y|x,data)在两个坐标轴上不停轮换采样就可以得到后验分布的样本集。 图5-9Bibbs采样点轨迹图 从图5-9的实验结果可以看出采样点大都集中在高密区。 当采样空间的维数比较高时，接受-拒绝采样和重要性采样都不适用，原因是当维数升高时很难找到适合的提议分布，采样效率比较差。没有比较就没有伤害，高维Gibbs采样算法的优势就体现出来了，其优点为： MCMC方法依赖于模拟的收敛性，下面介绍三种常用的收敛性诊断方法。 这种方法的思路是选取多个不同的初值，同时产生T条马尔科夫链，记第j条链的标准差估计为Sj，链内方差的均值为W，链间方差为B，则： R的值趋近1，则表明MCMC模拟收敛，且比较稳定。 GewekeTest由一系列的Z检验组成，其基本思路是：先对所有样本(假设有M个)做一次Z检验，然后去掉最前面的c个样本，用剩余的M-c个样本重复上述检验，再继续去掉最前面的c个样本，用剩余的M-2c个样本重复上述检验，依次类推，重复K次这样的检验，直到M-Kc。 该证明内容我们用定理来表述 定理：对于n≥0，i，j∈S,(S是状态空间)，则有 该定理可通过数学归纳法和Chapman-Kolmogorov等式来证明。 大家都知道，对于矩阵乘法，下面式子显然成立。 其中n≥0,m≥0. Chapman-Kolmogorov等式 其中n≥0,m≥0，i，j∈S，S是状态空间。 下面只需证明Chapman-Kolmogorov等式，然后就可以用数学归纳法证明我们的定理了。 证明方法一： 证明方法二： 根据该等式用数学归纳法证明定理略。 MCNP(MonteCarloNParticleTransportCode)是由美国洛斯阿拉莫斯国家实验室(LosAlamosNationalLaboratory)开发的基于蒙特卡罗方法的用于计算三维复杂几何结构中的中子、光子、电子或者耦合中子/光子/电子输运问题的通用软件包，也具有计算核临界系统（包括次临界和超临界系统）本征值问题的能力。 目前，MCNP以其灵活、通用的特点以及强大的功能被广泛应用于辐射防护与射线测定、辐射屏蔽设计优化、反应堆设计、（次）临界装置实验、医学以及检测器设计与分析等学科领域，并得到一致认可。 1906年AndreyAndreyevichMarkov引入了马尔可夫链的概念 Markov,A.A.(1906).Rasprostraneniezakonabol’shihchiselnavelichiny,zavisyaschiedrugotdruga.IzvestiyaFiziko-matematicheskogoobschestvapriKazanskomuniversitete,15(135-156),18 1953年梅特罗波利斯将蒙特卡洛方法引入马尔可夫链 1957年马尔可夫决策过程被提出 Bellman,R.(1957).AMarkoviandecisionprocess.JournalofMathematicsandMechanics,679-684. 1966年LeonardE.Baum等学者提出了隐马尔可夫模型（HiddenMarkovModel，HMM） Baum,L.E.;Petrie,T.(1966).StatisticalInferenceforProbabilisticFunctionsofFiniteStateMarkovChains.TheAnnalsofMathematicalStatistics.37(6):1554–1563. 1975年Baker将HMM用于语音识别 Baker,J.(1975).TheDRAGONsystem—Anoverview.IEEETransactionsonAcoustics,Speech,andSignalProcessing.23:24–29. 1980年《Markovrandomfieldsandtheirapplications》出版，详细描述了马尔可夫随机场的理论和应用 Kindermann,R.,&Snell,J.L.(1980).Markovrandomfieldsandtheirapplications(Vol.1).AmericanMathematicalSociety. 1984年Stuart和DonaldGeman兄弟描述了Gibbs抽样算法 Geman,S.;Geman,D.(1984).StochasticRelaxation,GibbsDistributions,andtheBayesianRestorationofImages.IEEETransactionsonPatternAnalysisandMachineIntelligence.6(6):721–741. 1988年JudeaPearl提出马尔可夫毯的概念 Pearl,J.(2014).*Probabilisticreasoninginintelligentsystems:networksofplausibleinference*.MorganKaufmann. 1989年JamesD.Hamilton1989年机应用了制转换模型 Hamilton,J.(1989).Anewapproachtotheeconomicanalysisofnonstationarytimeseriesandthebusinesscycle.Econometrica.57(2):357–84. 20世纪80年代HMM开始用于分析生物序列（DNA） Bishop,M.andThompson,E.(1986).MaximumLikelihoodAlignmentofDNASequences.JournalofMolecularBiology.190(2):159–165. 1991年Lovejoy研究了部分可观测马尔可夫决策过程（POMDP） Lovejoy,W.S.(1991).AsurveyofalgorithmicmethodsforpartiallyobservedMarkovdecisionprocesses.AnnalsofOperationsResearch.28(1):47–65. 1995年D.P.Bertsekas和J.N.Tsitsiklis讨论了一类用于不确定条件下的控制和顺序决策的动态规划方法 BertsekasD.P.andTsitsiklis,J.N.(1995).Neuro-dynamicprogramming:anoverview.Proceedingsof199534thIEEEConferenceonDecisionandControl.1:560-564. 1998年S.Brin和L.Page提出PageRank算法 Page,L.(1998).Thepagerankcitationranking:bringingordertotheweb.StanfordDigitalLibrariesWorkingPaper,9(1),1-14. 2017年YoshuaBengio等研究者提出了GibbsNet Lamb,A.etal.(2017).GibbsNet:IterativeAdversarialInferenceforDeepGraphicalModels.arXiv:1712.04120. 2017年JiamingSong,ShengjiaZhao和StefanoErmon研究了生成对抗的训练方法来对马尔可夫链（Markovchain）的转移算子（transitionoperator）进行学习 Song,J.;Zhao,S.;Ermon,S.(2017).GENERATIVEADVERSARIALLEARNINGOFMARKOVCHAINS.ICLR. 1970年Hastings完成了MCMC方法史上里程碑的论文-MonteCarloSamplingMethodsUsingMarkovChainsandtheirApplications<<马尔可夫链蒙特卡洛抽样方法及其应用>>.解决了核物理研究中碰到的粒子分布高纬计算问题，攻克了常规蒙特卡洛方法遇到的维度瓶颈问题，将Metropolis算法一般化，并推广为一种统计模拟工具，形成了Metropolis-Hastings算法。 1986年AdrianSmith做了关于分层模型的系列学术演讲。 1987年Tanner和Wong在论文中采用基于多个条件分布进行模拟的方法，这种思路等价于从联合分布进行模拟，基本具备了Gibbs采样的雏形。 1989年AdrianSmith第一次详细阐释了Gibbs抽样的本质。 1990年Gelfand和Smith发表论文Sampling-basedapproachestocalculationmarginaldensities.<<基于抽样的边际分布计算方法>>，将Gibbs抽样的本质阐述得更为深刻和完整，成为主流统计学界大规模使用MCMC方法的真正起点。 1995年Green提出得逆跳的马尔可夫链蒙特卡洛模拟(ReversibleJumpMarkovChainMonteCarlo,RJMCMC)，可被视为MCMC方法第二次革命的开端，所构建的马氏链不仅可以在一个模型的参数空间内进行转移，还可以在不同模型(维度可以不同)之间跳跃，从而为贝叶斯模型选择提供了强大工具。