MCMC(一)蒙特卡罗方法（随机模拟）
如何用蒙特卡罗方法来随机模拟求解一些复杂的连续积分或者离散求和。
例子：正方形、内切圆中撒点
rejection sampling
https://en.wikipedia.org/wiki/Rejection_sampling#Theory
使用接受-拒绝采样，我们可以解决一些概率分布不是常见的分布的时候，得到其采样集并用蒙特卡罗方法求和的目的。但是接受-拒绝采样也只能部分满足我们的需求，在很多时候我们还是很难得到我们的概率分布的样本集。比如：
| （1）对于一些二维分布$p(x,y)$，有时候我们只能得到条件分布$p(x | y)$和$p(y | x)$,却很难得到二维分布$p(x,y)$一般形式，这时我们无法用接受-拒绝采样得到其样本集。 |
（2）对于一些高维的复杂非常见分布$p(x_1,x_2,…,x_n)$，我们要找到一个合适的$q(x)$和$k$非常困难。
从上面可以看出，要想将蒙特卡罗方法作为一个通用的采样模拟求和的方法，必须解决如何方便得到各种复杂概率分布的对应的采样样本集的问题。而我们下一篇要讲到的马尔科夫链就是帮助找到这些复杂概率分布的对应的采样样本集的白衣骑士。
MCMC(二)马尔科夫链
马尔科夫链在很多时间序列模型中得到广泛的应用，比如循环神经网络RNN，隐式马尔科夫模型HMM等，当然MCMC也需要它。
细致平稳条件：如果非周期马尔科夫链的状态转移矩阵P和概率分布π(x)对于所有的$i,j$满足： \(\pi(i)P(i,j)=\pi(j)P(j,i)\) 则称概率分布π(x)是状态转移矩阵P的平稳分布。
如果我们得到了某个平稳分布所对应的马尔科夫链状态转移矩阵，我们就很容易采用出这个平稳分布的样本集。
假设我们任意初始的概率分布是$π_0(x)$, 经过第一轮马尔科夫链状态转移后的概率分布是$π_1(x)$，。。。第i轮的概率分布是$π_i(x)$。假设经过n轮后马尔科夫链收敛到我们的平稳分布π(x)，即：
\(\pi_n(x)=\pi_{n+1}(x)=\pi_{n+2}(x)=\cdots=\pi(x)\) 对于每个分布πi(x)，我们有： \(\pi_i(x)=\pi_{i-1}(x)P=\pi_{i-2}(x)P^2=\pi_0(x)P^i\) 现在我们可以开始采样了，首先，基于初始任意简单概率分布比如高斯分布$π_0(x)$采样得到状态值$x_0$，基于条件概率分布$P(x|x_0)$采样状态值$x_1$，一直进行下去，当状态转移进行到一定的次数时，比如到n次时，我们认为此时的采样集$(x_n,x_{n+1},x_{n+2},\cdots)$即是符合我们的平稳分布的对应样本集，可以用来做蒙特卡罗模拟求和了。
基于马尔科夫链的采样过程
（1）输入马尔科夫链状态转移矩阵P，设定状态转移次数阈值$n_1$，需要的样本个数$n_2$
（2）从任意简单概率分布采样得到初始状态值$x_0$
（3）for t=0 to $n_1+n_2−1$: 从条件概率分布$P(x|x_t)$中采样得到样本$x_{t+1}$
样本集$(x_{n1},x_{n_1+1},…,x_{n_1+n_2−1})$即为我们需要的平稳分布对应的样本集。
如果假定我们可以得到我们需要采样样本的平稳分布所对应的马尔科夫链状态转移矩阵，那么我们就可以用马尔科夫链采样得到我们需要的样本集，进而进行蒙特卡罗模拟。但是一个重要的问题是，随意给定一个平稳分布π,如何得到它所对应的马尔科夫链状态转移矩阵P呢？这是个大问题。我们绕了一圈似乎还是没有解决任意概率分布采样样本集的问题。
幸运的是，MCMC采样通过迂回的方式解决了上面这个大问题，我们在下一篇来讨论MCMC的采样，以及它的使用改进版采样: M-H采样和Gibbs采样.
MCMC(三)MCMC采样和M-H采样
在M-H采样中通过引入接受率使细致平稳条件满足。
目标矩阵P可以通过任意一个马尔科夫链状态转移矩阵Q乘以α(i,j)得到。
α(i,j)我们有一般称之为接受率。取值在[0,1]之间，可以理解为一个概率值。即目标矩阵P可以通过任意一个马尔科夫链状态转移矩阵Q以一定的接受率获得。这个很像接受-拒绝采样，那里是以一个常用分布通过一定的接受-拒绝概率得到一个非常见分布，这里是以一个常见的马尔科夫链状态转移矩阵Q通过一定的接受-拒绝概率得到目标转移矩阵P,两者的解决问题思路是类似的。
MCMC采样过程
（1）输入我们任意选定的马尔科夫链状态转移矩阵Q，平稳分布π(x)，设定状态转移次数阈值$n_1$，需要的样本个数$n_2$ （2）从任意简单概率分布采样得到初始状态值$x_0$ （3）for t=0 to $n_1+n_2−1$:
(a) 从条件概率分布$Q(x|x_t)$中采样得到样本$x^*$ (b) 从均匀分布采样u∼uniform[0,1] (c) 如果$u<\alpha(x_t,x^∗)=\pi(x^∗)Q(x^∗,x_t)$, 则接受转移$x_t→x^∗$，即$x_{t+1}=x^∗$ (d) 否则不接受转移，即$x_{t+1}=x_t$ 样本集$(x_{n_1},x_{n_1+1},…,x_{n_1+n_2−1})$即为我们需要的平稳分布对应的样本集。
MCMC采样比较难在实际中应用，为什么呢？问题在上面第三步的c步骤，接受率这儿。由于$α(x_t,x^∗)$可能非常的小，比如0.1，导致我们大部分的采样值都被拒绝转移，采样效率很低。有可能我们采样了上百万次马尔可夫链还没有收敛，也就是上面这个n1要非常非常的大，这让人难以接受，怎么办呢？这时就轮到M-H采样出场了。
M-H(Metropolis-Hastings)采样
\[\alpha(i,j)=min\left\{\frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\right\}\]（1）输入我们任意选定的马尔科夫链状态转移矩阵Q，平稳分布π(x)，设定状态转移次数阈值$n_1$，需要的样本个数$n_2$ （2）从任意简单概率分布采样得到初始状态值$x_0$ （3）for t=0 to $n_1+n_2−1$:
(a) 从条件概率分布$Q(x|x_t)$中采样得到样本$x^*$ (b) 从均匀分布采样u∼uniform[0,1] (c) 如果$u<\alpha(x_t,x^∗)=min\left{\frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)},1\right}$, 则接受转移$x_t→x^∗$，即$x_{t+1}=x^∗$ (d) 否则不接受转移，即$x_{t+1}=x_t$ 样本集$(x_{n_1},x_{n_1+1},…,x_{n_1+n_2−1})$即为我们需要的平稳分布对应的样本集。
很多时候，我们选择的马尔科夫链状态转移矩阵$Q$如果是对称的，即满足
$Q(i,j)=Q(j,i)$,这时我们的接受率可以进一步简化为： \(\alpha(i,j)=min\left\{\frac{\pi(j)}{\pi(i)},1\right\}\) M-H采样完整解决了使用蒙特卡罗方法需要的任意概率分布样本集的问题，因此在实际生产环境得到了广泛的应用。
但是在大数据时代，M-H采样面临着两大难题：
（1） 我们的数据特征非常的多，M-H采样由于接受率计算式$\frac{\pi(j)Q(j,i)}{\pi(i)Q(i,j)}$的存在，在高维时需要的计算时间非常的可观，算法效率很低。同时α(i,j)一般小于1，有时候辛苦计算出来却被拒绝了。能不能做到不拒绝转移呢？
（2） 由于特征维度大，很多时候我们甚至很难求出目标的各特征维度联合分布，但是可以方便求出各个特征之间的条件概率分布。这时候我们能不能只有各维度之间条件概率分布的情况下方便的采样呢？
Gibbs采样解决了上面两个问题，因此在大数据时代，MCMC采样基本是Gibbs采样的天下
MCMC(四)Gibbs采样
二维Gibbs采样
二维Gibbs采样需要两个维度之间的条件概率。具体过程如下：
（1）输入平稳分布$π(x_1,x_2)$，设定状态转移次数阈值$n_1$，需要的样本个数$n_2$
（2）随机初始化初始状态值$x_1^{(0)}$和$x_2^{(0)}$
（3）for t=0 to $n_1+n_2−1$:
| (a) 从条件概率分布$P(x_2 | x_1^{(t)})$中采样得到样本$x_2^{t+1}$ |
| (b) 从条件概率分布$P(x_1 | x_2^{(t+1)})$中采样得到样本$x_1^{t+1}$ |
样本集$\left{(x_1^{(n_1)},x_2^{(n_1)}),(x_1^{(n_1+1)},x_2^{(n_1+1)}),…,(x_1^{(n_1+n_2−1)},x_2^{(n_1+n_2−1)})\right}$即为我们需要的平稳分布对应的样本集。
整个采样过程中，我们通过轮换坐标轴，采样的过程为：
$(x_1^{(1)},x_2^{(1)})→(x_1^{(1)},x_2^{(2)})→(x_1^{(2)},x_2^{(2)})→…→(x_1^{(n_1+n_2−1)},x_2^{(n_1+n_2−1)})$
多维Gibbs采样
| 上面的这个算法推广到多维的时候也是成立的。比如一个n维的概率分布$\pi(x_1,x_2,…x_n)$，我们可以通过在n个坐标轴上轮换采样，来得到新的样本。对于轮换到的任意一个坐标轴$x_i$上的转移，马尔科夫链的状态转移概率为$P(x_i | x_1,x_2,…,x_{i−1},x_{i+1},…,x_n)$，即固定n−1个坐标轴，在某一个坐标轴上移动。 |
具体的算法过程如下：
（1）输入平稳分布$π(x_1,x_2，…,x_n)$或者对应的所有特征的条件概率分布，设定状态转移次数阈值$n_1$，需要的样本个数$n_2$
（2）随机初始化初始状态值$(x_1^{(0)},x_2^{(0)},…,x_n^{(0)})$
（3）for t=0 to $n_1+n_2−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_4^{(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}$ |
样本集$\left{(x_1^{(n_1)},x_2(n_1),…,x_n(n_1)),…,(x_1^{(n_1+n_2−1)},x_2^{(n_1+n_2−1)},…,x_n^{(n_1+n_2−1)})\right}$即为我们需要的平稳分布对应的样本集。
整个采样过程和Lasso回归的坐标轴下降法算法非常类似，只不过Lasso回归是固定n−1个特征，对某一个特征求极值。而Gibbs采样是固定n−1个特征在某一个特征采样。
同样的，轮换坐标轴不是必须的，我们可以随机选择某一个坐标轴进行状态转移，只不过常用的Gibbs采样的实现都是基于坐标轴轮换的。
由于Gibbs采样在高维特征时的优势，目前我们通常意义上的MCMC采样都是用的Gibbs采样。当然Gibbs采样是从M-H采样的基础上的进化而来的，同时Gibbs采样要求数据至少有两个维度，一维概率分布的采样是没法用Gibbs采样的,这时M-H采样仍然成立。
有了Gibbs采样来获取概率分布的样本集，有了蒙特卡罗方法来用样本集模拟求和，他们一起就奠定了MCMC算法在大数据时代高维数据模拟求和时的作用。
References
- https://www.cnblogs.com/pinard/p/6625739.html 刘建平
- https://cosx.org/2013/01/lda-math-mcmc-and-gibbs-samplingLDA-math-MCMC 和 Gibbs Sampling 靳志辉
- https://zhuanlan.zhihu.com/p/37121528 马尔可夫链蒙特卡罗算法（MCMC）
