随机模拟方法简介
随机模拟方法又称为蒙特卡罗方法（Monte Carlo Method）。蒙特卡洛模拟方法的原理是当问题或对象本身具有概率特征时，可以用计算机模拟的方法产生抽样结果，根据抽样计算统计量或者参数的值；随着模拟次数的增多，可以通过对各次统计量或参数的估计值求平均的方法得到稳定结论。由于涉及到时间序列的反复生成，蒙特卡洛模拟法是以高容量和高速度的计算机为前提条件的，因此只是在近些年才得到广泛推广。
从上面的描述我们不难看出随机模拟方法主要是针对那些确定算法不好解或者解不出来的情况。因此是一种典型的寻求近似解的方法。
基本思路
- 针对实际问题建立一个简单易行的概率统计模型，使问题所求的解为该模型的概率分布或者数字特征，比如：某个事件的概率或者是某个随机变量的期望值。
- 对模型中的随机变量建立抽样方法，在计算机上进行模拟试验，得到足够的随机抽样，并对相关事件进行统计
- 对试验结果进行分析，给出所求解的估计及其精度（方差）的估计
蒙特卡罗方法在数学中的应用
通常蒙特卡洛方法通过构造匹配一定规则的随机数来解决数学上的各种问题。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，蒙特卡洛方法是一种有效的求出数值解的方法。一般蒙特卡洛方法在数学中最常见的应用就是蒙特卡洛积分。下面是蒙特卡罗方法的两个简单应用：
圆周率
蒙特卡洛方法可用于近似计算圆周率：让计算机每次随机生成两个0到1之间的数，看以这两个实数为横纵坐标的点是否在单位圆内。生成一系列随机点，统计单位圆内的点数与总点数，（圆面积和正方形面积之比为PI:4，PI为圆周率），当随机点获取越多时，其结果越接近于圆周率（然而准确度仍有争议：即使取10的9次方个随机点时，其结果也仅在前4位与圆周率吻合）。用蒙特卡洛方法近似计算圆周率的先天不足是：第一，计算机产生的随机数是受到存储格式的限制的，是离散的，并不能产生连续的任意实数；上述做法将平面分区成一个个网格，在空间也不是连续的，由此计算出来的面积当然与圆或多或少有差距。
蒙特卡罗积分
举个简单的例子，假设我们要求f(x)的积分∫baf(x)dx
但是f(x)的形式比较复杂不好求积分，使用蒙特卡罗积分方法转化为∫baf(x)q(x)q(x)dx
把q(x)看作是x在区间[a,b]内的概率分布，前面的分数部分看做一个函数，然后再q(x)下抽取n个样本，当n足够大时，可以利用均值来近似（大数定理），1n∑if(xi)p(xi)
因此只要q(x)比较容易采样就行。
随机模拟三个要素：随机数，逻辑模型，反复试验。其基本思路就是要把待解决的问题转化为一种可以通过某种采样方法可以解决的问题，因此随机模拟方法的核心就是如何对一个概率分布得到样本，即抽样（sampling）。
常见的采样方法
直接采样
通过对均匀分布采样，实现对任意分布采样。
一般而言均匀分布 Uniform(0,1)的样本是相对容易生成的。 通过线性同余发生器可以生成伪随机数，我们用确定性算法生成[0,1]之间的伪随机数序列后，这些序列的各种统计指标和均匀分布 Uniform(0,1) 的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。
而我们常见的概率分布，无论是连续的还是离散的帆布，都可以基于Uniform(0,1)的样本生成。
直接采样步骤：
- 从Uniform(0,1)随机产生一个样本z
- 另z=h(y)，其中h(y)为y的累积概率分布CDF
- 计算y=h?1(z)
- 结果y为对p(y)的采样
注：需要知道累积概率分布的解析表达式，且累积概率分布函数存在反函数
但是如果h(x)不能确定或者没有无法解析求逆则直接抽样不再合适。对于复杂的现实模型其实是不常用的。
接受-拒接采样
简称拒绝采样，基本思想：假设我们需要对一个分布f(x)进行采样，但是却很难直接进行采样，所以我们想通过另外一个容易采样的分布g(x)的样本，用某种机制去除掉一些样本，从而使得剩下的样本就是来自与所求分布f(x)的样本。
- 给定目标分布密度π(x)
- 建议密度（proposal density）q(x)和常数M，使得
对q(x)采样比较容易
q(x)的形状比较接近π(x)
对任意x，有π(x)≤Mq(x)（包络原则）
- 通过对q(x)采样实现对π(x)采样
采样过程：
产生样本X～q(x)，和U～Uniform[0,1]
若U≤π(X)/Mq(x)，则接受X，接受的样本服从分布π(x)
等价于：
Y=Mq(X)U，若Y≤π(X)，则接受X
在高维空间，接受-拒绝采样会出现两个问题，一是合适的proposal density q(x)比较难找，而是很难确定一个准确的M值。这两个问题会导致拒绝率很高，无用计算增加。
重要性采样
通过从已知采样的概率q(x)采样，近似积分
I=∫f(x)π(x)dx=∫f(x)π(x)q(x)q(x)dx
从q(x)中抽取N个样本，上述式子就约等于1N∑Nif(xi)π(xi)q(xi)。这相当于给每一个样本赋予了一个权重w(xi)=π(xi)q(xi)，q(x)大意味着概率大，那么N里面含有这个样本的xi就多，即这些样本的权重大，所以称之为重要性抽样。下面这个链接里关于重要性采样的理解写的很形象：
http://blog.sina.com.cn/s/blog_4e5740460100cw5b.html
采样过程：
- 选择一个容易抽样的分布q(x)，产生N个样本X1,...,XN～q(x)
- 近似解：1N∑Nif(xi)π(xi)q(xi)
注：q(x)的形状应与f(x)π(x)的形状足够近似，且q(x)的尾部比π(x)的尾部厚，估计的方差才不为无限大。
在高维空间中合适的q(x)很难找到。
MCMC采样方法
重要性采样和接受-拒绝采样都在q(x)与π(x)很相似的时候才表现好。而且在高维空间问题中，标准的采样方法会失败：
- 接受-拒绝采样：维数增高，拒绝率→100%
- 重要性采样：大多数的样本权重→0
而且上述两种采样方法都是独立采样的，效率较低，MCMC采样方法是关联采样，即下一个样本与这个样本有关系。在蒙特卡洛模拟中，我们在后验分布中抽取样本，当这些样本独立时，利用大数定律样本均值会收敛到期望值。如果得到的样本是不独立的，那么就要借助于马尔科夫链进行抽样。MCMC（Markov Chain Monte Carlo)方法就是为了这个目的而诞生的。
MCMC方法的基本思想是：通过构建一个markov chain使得该markov chain的稳定分布是我们所要采样的分布f(x)。如果这个markov chain达到稳定状态，那么来自这个chain的每个样本都是f(x)的样本，从而实现抽样的目的。这里存在一个核心问题，如何构建满足要求的markov chain？
马尔科夫链基础
皮姐姐的博客写过这部分，写的很清楚了，直接放在这里吧
http://blog.csdn.net/pipisorry/article/details/46618991
参考文献【1】中也有讲得很详细。
总之这部分的重点是马氏链的平稳分布，这是MCMC算法的核心基础
马氏链的收敛性质主要由转移矩阵P 决定, 所以基于马氏链做采样的关键问题是如何构造转移矩阵P,使得平稳分布恰好是我们要的分布π(x)。
这里的关键是细致平稳条件（具体内容在上面两个参考链接里都有）。细致平稳条件的好处，就是我们能控制马尔科夫链收敛到我们指定的分布。
构造满足条件的马氏链
假设我们已经有一个转移矩阵为Q的马氏链，q(i,j)表示从状态i到j的概率，通常情况下p(i)q(i,j)≠p(j)q(j,i)也就是细致平稳条件不成立，所以p(x)不太可能是这个马氏链的平稳分布。因此，我们就要对这个马氏链进行一些小改造，使他满足细致平滑条件。例如我们引入一个α(i,j)使得
p(i)q(i,j)α(i,j)=p(j)q(j,i)α(j,i)
最简单的我们按照对称性取
α(i,j)=p(j)q(j,i),α(j)=p(i)q(i,j)
就可以使得上式成立，可以看做具有新的转移矩阵Q′的马氏链，其中
Q′(i,j)=q(i,j)α(i,j),Q′(j,i)=q(j,i)α(j,i).
因此马氏链Q′满足细致平稳条件，且平稳分布是p(x)。
在改造 Q的过程中引入的 α(i,j)称为接受率，物理意义可以理解为在原来的马氏链上，从状态 i 以q(i,j) 的概率转跳转到状态j 的时候，我们以α(i,j)的概率接受这个转移，于是得到新的马氏链Q′的转移概率为q(i,j)α(i,j)。
注：当按照上面介绍的构造方法把Q–>Q′后，就不能保证Q′是一个转移矩阵了，即Q′的每一行加和为1。这时应该在当 j != i 的时候概率Q′(i, j) 就如上处理， 当j = i 的时候， Q′(i, i) 应该设置Q′(i, i) = 1- 其它概率之和，归一化概率转移矩阵。
MCMC采样算法(参考文献1)
上述过程中 p(x),q(x|y) 说的都是离散的情形，事实上即便这两个分布是连续的，以上算法仍然是有效，于是就得到更一般的连续概率分布 p(x)的采样算法，而 q(x|y) 就是任意一个连续二元概率分布对应的条件分布。
以上的 MCMC 采样算法已经能很漂亮的工作了，不过它有一个小的问题：马氏链Q在转移的过程中的接受率α(i,j)可能偏小，这样采样过程中马氏链容易原地踏步，拒绝大量的跳转，这使得马氏链遍历所有的状态空间要花费太长的时间，收敛到平稳分布p(x)的速度太慢。有没有办法提升一些接受率呢?
假设 α(i,j)=0.1,α(j,i)=0.2
, 此时满足细致平稳条件，于是
p(i)q(i,j)×0.1=p(j)q(j,i)×0.2
上式两边扩大5倍，我们改写为
p(i)q(i,j)×0.5=p(j)q(j,i)×1
看，我们提高了接受率，而细致平稳条件并没有打破！这启发我们可以把细致平稳条件式中的α(i,j),α(j,i) 同比例放大，使得两数中最大的一个放大到1，这样我们就提高了采样中的跳转接受率。所以我们可以取
α(i,j)=min{p(j)q(j,i)p(i)q(i,j),1}
于是，经过对上述MCMC 采样算法中接受率的微小改造，我们就得到了如下教科书中最常见的 Metropolis-Hastings 算法。
Metropolis-Hasting算法
对于分布p(x)，我们构造转移矩阵Q′使其满足细致平稳条件
p(x)Q′(x→y)=p(y)Q′(y→x)
此处x并不要求是一维的，对于高维空间的p(X)，如果满足细致平稳条件
p(X)Q′(X→Y)=p(Y)Q′(Y→X)
那么以上的M-H算法一样有效
Gibbs Sampling
对于高维的情况，由于接受率α（通常小于1）的存在，M-H算法效率不够高，能否找到一个转移矩阵Q使得接受率α=1呢？
先看二维的情况，假设有一个概率分布p(x,y)，考察x坐标相同的两个点A(x1,y1),B(x1,y2)，我们发现
p(x1,y1)p(y2|x1)=p(x1)p(y1|x1)p(y2|x1)
p(x1,y2)p(y1|x1)=p(x1)p(y1|x2)p(y1|x1)
所以得到
p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)
即
p(A)p(y2|x1)=p(B)p(y1|x1)
基于以上等式，我们发现，在x=x1这条平行于y轴的直线上，如果使用条件分布p(y|x_1)作为任何两个点之间的转移概率，那么任何两个点之间的转移满足细致平稳条件。同样的，如果我们在y=y1这条直线上任意取两个点A，C则有
p(A)p(x2|y1)=p(C)p(x1|y1)
于是我们可以如下构造平面上任意两点之间的概率转移矩阵Q
有了如上的转移矩阵Ｑ，我们很容易验证对平面上的任意两点Ｘ，Ｙ，满足细致平稳条件
p(X)Q(X→Y)=p(Y)Q(Y→X)
于是这个二维空间上的马氏链将收敛到平稳分布p(x,y)。而这个算法就称为 Gibbs Sampling 算法,是 Stuart Geman 和Donald Geman 这两兄弟于1984年提出来的，之所以叫做Gibbs Sampling 是因为他们研究了Gibbs random field, 这个算法在现代贝叶斯分析中占据重要位置。
二维Gibbs Sampling算法
Gibbs Sampling算法中的马氏链转移
以上采样过程中，如图所示，马氏链的转移只是轮换的沿着坐标轴 x轴和y轴做转移，于是得到样本 (x0,y0),(x0,y1),(x1,y1),(x1,y2),(x2,y2),? 马氏链收敛后，最终得到的样本就是 p(x,y) 的样本，而收敛之前的阶段称为 burn-in period。额外说明一下，我们看到教科书上的 Gibbs Sampling 算法大都是坐标轴轮换采样的，但是这其实是不强制要求的。最一般的情形可以是，在t时刻，可以在x轴和y轴之间随机的选一个坐标轴，然后按条件概率做转移，马氏链也是一样收敛的。轮换两个坐标轴只是一种方便的形式。
n 维Gibbs Sampling算法
以上的过程我们很容易推广到高维的情形，如果x1 变为多维情形x1，可以看出推导过程不变，所以细致平稳条件同样是成立的
p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)
此时转移矩阵Q由条件分布p(y|x1)定义。上式知识说明了一根坐标轴的情形，二维情形类似。所以n维空间中对于概率分布p(x1,x2,...,xn)可以如下定义转移矩阵
- 如果当前状态为(x1,x2,...,xn)，马氏链转移的过程中，只能沿着坐标轴做转移，沿着xi这根坐标轴转移时，转移概率由条件概率p(xi|x1,...,xi?1,xi+1,...,xn)定义
- 其它无法沿着单根坐标轴进行的跳转，转移概率都设置为 0
Gibbs抽样就是在这m个条件分布中迭代产生样本。于是我们可以把吉布斯采样由二维推广到n维：
[图片]
以上算法收敛后，得到的就是概率分布p(x1,x2,...,xn)的样本，当然这些样本并不独立，但是我们此处要求的是采样得到的样本符合给定的概率分布，并不要求独立。
同样的，在以上算法中，坐标轴轮换采样不是必须的，可以在坐标轴轮换中引入随机性，这时候转移矩阵 Q 中任何两个点的转移概率中就会包含坐标轴选择的概率，而在通常的 Gibbs Sampling 算法中，坐标轴轮换是一个确定性的过程，也就是在给定时刻t，在一根固定的坐标轴上转移的概率是1。
吉布斯采样可以看做是M-H算法的一个特例，即接受率α=1的情况。证明如下：
考虑一个M-H采样的步骤，它涉及到变量zk，剩余变量z?k保持不变。同时，从z到z?的转移概率为qk(z?|z)=p(z?k|z?k)。我们注意到z??k=z?k，因为再采样步骤中，向量的各个元素都不变。又p(z)=p(zk|z?k)p(z?k)，因此，确定M-H算法中的接受概率为
A(z?,z)=p(z?)qk(z|z?)p(z)qk(z?|z)=p(z?k|z??k)p(z??k)p(zk|z??k)p(zk|z?k)p(z?k)p(z?k|z?k)=1
收敛性判断
当然无论是metropolis-hasting算法还是gibbs算法，都有一个burn in的过程，所谓burn in的过程就是因为这个两个算法本身都是markov chain的算法，要达到稳定状态需要一定的步骤才能达到，所以需要一个burn in过程，只有在达到平衡状态时候得到的样本才能是平衡状态时候的目标分布的样本，因此，在burn in过程中产生的样本都需要被舍弃。如何判断一个过程是否达到了平衡状态还没有一个成熟的方法来解决，目前常见的方法是看是否状态已经平稳（例如画一个图，如果在较长的过程中，变化已经不大，说明很有可能已经平衡）当然这个方法并不能肯定一个状态是否平衡，你可以举出反例，但是却是实际中没有办法的办法。
MCMC方法依赖于产生的马氏链在t足够大时要收敛。
关于链的收敛有这样一些检验方法。
（1）图形方法 这是简单直观的方法。我们可以利用这样一些图形：
（a）迹图（trace plot）：将所产生的样本对迭代次数作图，生成马氏链的一条样本路径。如果当t足够大时，路径表现出稳定性没有明显的周期和趋势，就可以认为是收敛了。
(b)自相关图（Autocorrelation plot）：如果产生的样本序列自相关程度很高，用迹图检验的效果会比较差。一般自相关随迭代步长的增加而减小，如果没有表现出这种现象，说明链的收敛性有问题。
（c）遍历均值图（ergodic mean plot）：MCMC的理论基础是马尔科夫链的遍历定理。因此可以用累积均值对迭代步骤作图，观察遍历均值是否收敛。
（2）蒙特卡洛误差
（3）Gelman-Rubin方法
From：
http://blog.csdn.net/sherrylml/article/details/51434549
Ref:
[1] http://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/
[2] http://blog.csdn.net/xianlingmao/article/details/7768833
[3]http://www.cnblogs.com/xbinworl/p/4266146.html?utm_source=tuicool&utm_medium=referral
[4] https://site.douban.com/182577/widget/notes/10567181/note/292072927/
[5]Pattern Recognition and Machine Learning. Chapter11 Sampling Method