在复杂函数求期望、非线性函数近似等问题中，需要从一些特定的函数中采样。而不同于高斯分布、均匀分布和Gamma分布等较为简单的分布，这些分布往往难以直接采样，因此需要从其他角度设计采样方法。这里介绍几种常用的方法。
假设 \(p(x)\) 难以直接采样，\(q(x)\)是一个比较容易采样的分布，如高斯、均匀分布，且正整数\(M\)使得\(p(x)/(Mq(x))<1\)。
则接收-拒绝采样的流程：
（从上步骤可以看出，得到一个有效粒子需要平均\(M\)次采样）
可以理解为，\(p(x)\)表示一个超多面体包围的区域（如圆形面积），\(Mq(x)\)为另一个完全包含\(p(x)\)的超多面体区域（如一个包含圆形的正方形），则当随机粒子\(\mu<p(x)/(Mq(x))\)，则该粒子位于\(p(x)\)所包围区域中。为了使得接受率尽可能高，\(M\)应该在满足\(p(x)/(Mq(x))<1\)的条件下尽量小。该方法的缺点是：找到合适的\(q(x)\)很难；且接受概率可能很低，使得采样效率低下。
假设\(p(x)\)难以直接采样，\(q(x)\)是一个比较容易采样的分布（称为proposal function），则对\(p(x)\)的采样可以转换成从\(q(x)\)中采样的粒子\(x_i\)的权重和,
其中 \(w_i = \frac{p({x_i})}{q({x_i})}\)为重要性权重（Importance Weight）。重要性采样的关键在于寻找合适的proposal function，通常\(q(x)\)与\(p(x)\)越相似越好。在贝叶斯估计中，一般可以将预测分布作为proposal function。其他还有各种对该函数设计的研究，这里不展开。粒子滤波（Particle filters or Sequential Monte Carlo (SMC)）被用于在给定观测的情况下，估计状态的后验密度，其可以看作是重要性采样的递归版本（如序贯重要性重采样）。
MCMC算法是基于马尔科夫链（Markov chain）的从一个复杂的分布中随机采样的方法。当下比较流行的MCMC方法有Metropolis-Hastings算法，Gibbs采样。他们都是将所要采样的分布建模为马尔科夫链的平稳分布（equilibrium distribution）。
定义（Markov chains）：下一状态的概率分布只由当前状态决定，与时间序列中它前面的事件无关，即
假设初始的分布为\(p^{(0)}(x)= [p(x_1^{(0)}),p(x_2^{(0)}),...]\)，系统的马尔科夫转移函数\(f(x|x‘)\)表示从\(x‘\)转移至\(x\)的概率。经\(n \to \infty\)次转移之后的系统分布为\(p(x)=p^{(n)}(x)\)\(=p^{(0)}(x){f^n}(x|x‘)\)，其中\(p(x) = [{p(x_1)},{p(x_2)},...]\)。
如果对于任意的状态\(x_i\)和\(x_j\)，概率分布满足\(p(x_j) = \sum\limits_i p(x_i)f(x_j|x_i)\)，则\(p(x)\)为该马尔科夫链的平稳分布。且若\(f(x|x‘)\)为离散的，则
这一性质满足的条件是该马尔科夫链是可约且非周期的，平常遇到的大部分马尔科夫链均满足该条件。
回到抽样问题，若希望从一个难以直接采样的分布中进行采样，则MCMC的思想是：构造一个转移函数为\(f(x|x‘)\)的马尔科夫链，使得其平稳分布就是所要采样的分布，则可以从任意的初始分布出发，经过多次Markov转移，马氏链收敛之后的采样就是源自其平稳分布（即要采样的分布）的采样。
（这个想法在1953年被 Metropolis想到了，为了研究粒子系统的平稳性质，Metropolis 考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法，即Metropolis算法，并在最早的计算机上编程实现。Metropolis 算法是首个普适的采样方法，并启发了一系列 MCMC方法，所以人们把它视为随机模拟技术腾飞的起点。 Metropolis的这篇论文被收录在《统计学中的重大突破》中， Metropolis算法也被遴选为二十世纪的十个最重要的算法之一。）
但是，在已知需要采样的分布（此处即，平稳分布），如何构造转移转移函数\(f(\cdot|\cdot)\)？可以从细致平稳条件（detailed balance condition）出发来计算\(f(\cdot|\cdot)\)。
定理: [细致平稳条件] 如果非周期马尔科夫链的状态转移函数\(f\)和概率分布\(p(x)\)对于所有的\(x_i\)和\(x_j\)，满足：\(p(x_i)f(x_j|x_i) = p(x_j)f(x_i|x_j)\)，则\(p(x)\)是该马尔科夫链的平稳分布。
所以只要找到使得概率分布\(p(x)\)满足细致平稳条件的函数\(f\)即可。假设对于任意的马尔科夫转移函数\(g(\cdot|\cdot)\)，由于他一般不满足细致平稳条件，则
则引入因子\(\alpha (i,j) < 1\)，使得
满足以上的\(\alpha\)可取为
因此转移矩阵\(f(x_j|x_i) = g(x_j|x_i)\alpha(i,j)\) 。这样，将任意的转移函数\(g(\cdot|\cdot)\)转换成了所希望得到的函数\(f(\cdot|\cdot)\)。以上\(\alpha(i,j)\)可理解为接受率，即以概率\(\alpha(i,j)\)接受“从状态\(x_i\)以概率\(g(x_j|x_i)\)转移至状态\({x_j}\)”这一事件。
综合以上，MCMC算法流程为：
需要采样的函数为\(p(x)\)，任意选定马尔科夫转移密度\(g(\cdot|\cdot)\)，设定转移次数阈值\(M\)，需要样本数为\(N\)。
初始化状态为\(x^0\)，
\(For \quad t=0:N+M-1\)
a) 采样得\(y \sim g(x|{x^t})\)
b) 从均匀分布\(U(0,1)\)采样得\(\mu\)
c) 如果\(\mu < \alpha (x^t,y) = p(y)g(x^t|y)\)，则接受转移\(x^{t + 1}= y\)，否则不接受转移且返回步骤a）。
样本集\(({x^M},{x^{M + 1}},...,{x^{N + M - 1}})\)为采样的样本集。
注：第三步其实相当于接受-拒绝采样：实际要从\(f(x_j|x_i) = g(x_j|x_i)\alpha(i,j)\)，对于\(\alpha(\cdot) < 1\),根据接受-拒绝采样方法，构造函数\(\frac{g(x_j|x_i)\alpha(i,j)}{g(x_j|x_i)} = \alpha(i,j)\)，则当\(\mu < \alpha(i,j)\)，接受。
以上算法的一个问题就是接受率可能非常小，导致大部分采样值被拒绝，采样效率十分低下。
M-H采样解决了MCMC采样中的采样效率低下的问题。
根据以上的细致平稳条件
当\(\alpha\)较小，对上式两端同比例放大，使得\(p(x_i)g(x_j|x_i)A(i,j)=p(x_j)g(x_i|x_j)\times 1\)，此时依然满足细致平稳条件，且接受概率为\(A(i,j) = \frac{p(x_j)g(x_i|x_j)}{p(x_i)g(x_j|x_i)}\)。又接受概率应当小于1，则可得到
M-H采样步骤如下:
需要采样的函数为\(p(x)\)，任意选定马尔科夫转移密度\(g(\cdot|\cdot)\)，设定转移次数阈值\(M\)，需要样本数为\(N\)。
初始化状态为\(x_0\)，
\(For \quad t=0:N+M-1\)
a) 采样得\(y \sim g(x|{x_t})\)
b) 从均匀分布\(U(0,1)\)采样得\(\mu\)
c) 如果 \(\mu<A(i,j)=\min \left( {1,\frac{p({x_j})g({x_i}|{x_j})}{p({x_i})g({x_j}|{x_i})}} \right)\)，则接受转移\({x_{t + 1}} = y\)，否则不接受转移且返回步骤a）。
样本集\(({x_M},{x_{M + 1}},...,{x_{N + M - 1}})\)为采样的样本集。
M-H采样的缺点在于1、在高维时，接受概率计算较为复杂；2、还是存在拒绝的可能性，导致算法收敛时间长；3、通常的贝叶斯问题中，M-H所需的特征的联合分布难以求得，而条件分布相对于好求。Gibbs采样可以解决这些问题。
Gibbs采样是Metropolis-Hastings采样算法的特殊形式，即找到一个已知的分布，使得接受率1，这样，每次的采样都会被接受，可以提高MCMC的收敛速度。
这里从二维分布开始介绍，并推广到高维分布。假设\(p(X)=p(x,y)\)是一个二维联合分布并希望从中采样，对于两点\(A:({x_A},{y_A})\)和\(B:({x_B},{y_B})\)，当第一维度相同\({x_A} = {x_B} = {x_0}\)，容易得到（对条件密度进行贝叶斯展开即可）
即\(p(A)p(y_B|x_0) = p(B)p(y_A|x_0)\)。所以当位于\(x = {x_0}\)这一条直线上的，使用\(p( \cdot |{x_0})\)作为状态转移函数，则该直线上的任意两点满足细致平稳条件。同理若两点的第二维度相同的点\({y_A} = {y_B} = {y_0}\)，则在\(y = {y_0}\)这条直线上的任意两点，满足转移函数为\(p( \cdot |{y_0})\)的细致平稳条件。
综上所述，对于任意两点\(A:({x_A},{y_A})\)和\(B:({x_B},{y_B})\)，构造分布\(p(X)\)的概率分布对应的马尔科夫转移密度为，
则它们满足细致平稳条件\(P(A)f(B|A) = P(B)f(B|A)\)，于是该二维空间上的马氏链将收敛到平稳分布\(p(X)\)。利用该状态转移密度，可以得到二维Gibbs采样方法，这个采样需要两个维度之间的条件概率。具体过程如下：
确定所需的采样的分布\(p(X)=p(x,y)\)，设定状态转移次数阈值\(M\)，以及需要的样本数\(N\)
随机化初始状态值\(x^0\)和\(y^0\).
\(for \quad t=0:N+M-1\)
a) 从条件概率分布\(f(y|{x^t})\)中采样得到\({y^{t + 1}}\)
b) 从条件概率分布\(P(x|{y^{t + 1}})\)中采样得到\({x^{t + 1}}\)
样本集\(\{ ({x^M},{y^M}),({x^{M + 1}},{y^{M + 1}}),...,({x^{M + N - 1}},{y^{M + N - 1}})\}\)为所需的样本集。
以上过程实际通过轮换坐标轴实现，即采样顺序在X,Y轴之间相互轮换。但是实际上轮换坐标轴不是必须的，我们可以随机选择某一个坐标轴进行状态转移，只不过常用的Gibbs采样的实现都是基于坐标轴轮换的。多维Gibbs采样是二维的推广，可以假设在n维坐标的轮换采样得到。\(N\)维空间函数采样具体过程如下：
假设需要从联合分布\(p(X)\)中获得\(N\)个粒子，其中\(X = \left\{ {x_1},...,{x_n}\right\}\)，其中将第\(i\)个粒子表示为\({X^{i}}\)。
随机化初始状态值\(X^0\)
\(for \quad t=0:N+M-1\)
a) 从条件分布\(P({x_1}|x_2^t,x_3^t,...,x_n^t)\)中采样得到\(x_1^{t + 1}\)
b) 从条件分布\(P({x_2}|x_1^{t + 1},x_3^t,...,x_n^t)\)中采样得到\(x_2^{t + 1}\)
c) ……
d) 从条件分布\(P({x_j}|x_1^{t + 1},x_2^{t + 1},...,x_{j - 1}^{t + 1},x_{j + 1}^t,...,x_n^t)\)中采样得到\(x_j^{t + 1}\)
e) ……
f) 从条件分布\(P({x_n}|x_1^{t + 1},x_2^{t + 1},....,x_{n - 1}^{t + 1})\)中采样得到\(x_n^{t + 1}\)
最终样本集为\({X^M},X^{M + 1},...,X^{M + N - 1}\)，其中\({X^j} =\left\{ x_1^j,...,x_n^j\right\}\)。
由于Gibbs采样在高维特征时的优势，目前我们通常意义上的MCMC采样都是用的Gibbs采样。当然Gibbs采样是从M-H采样的基础上的进化而来的，同时Gibbs采样要求数据至少有两个维度，一维概率分布的采样是没法用Gibbs采样的,这时M-H采样仍然成立。
原文：https://www.cnblogs.com/tobeingoodmood/p/13277014.html