MCMC方法整理与理解

首先感谢大神的整理
http://www.cnblogs.com/xbinworld/p/4266146.html
这篇博客从这几点引入，并介绍了MCMC方法

1. 蒙特卡洛数值积分
2. 均匀分布，Box-Muller 变换
3. Monte Carlo principle
4. 接受-拒绝抽样（Acceptance-Rejection sampling)
5. 重要性抽样(Importance sampling)
6. 马尔科夫链，马尔科夫稳态
7. MCMC——Metropolis-Hasting算法
8. MCMC——Gibbs Sampling算法

下面是我对这方面知识的学习过程

1, 蒙特卡洛数值积分

如果我们要求f(x)的积分，如而f(x)的形式比较复杂积分不好求，则可以通过数值解法来求近似的结果。常用的方法是蒙特卡洛积分：
$${\int }_a^b\frac{f(x)}{q(x)}q(x)\mathrm{d}x$$
这样把q(x)看做是x在区间内的概率分布，而把前面的分数部门看做一个函数，然后在q(x)下抽取n个样本，当n足够大时，可以用采用均值来近似：
因此只要q（x）比较容易采到数据样本就行了。随机模拟方法的核心就是如何对一个概率分布得到样本，即抽样（sampling）。下面我们将介绍常用的抽样方法。

[注]关于上述公式的推导和理解

刚开始很难理解，为什么$${\int }_a^{b}f(x)\mathrm{d}x={\int }_a^b\frac{f(x)}{q(x)}q(x)\mathrm{d}x$$
而且可以用q(x)的概率分布去近似f(x)的概率分布，这里其实运用到了蒙特卡洛求定积分的期望法，蒙特卡洛求定积分比较好理解的是投点法，期望法的推导如下：
首先我们设$g^{\ast }(x)=\frac{f(x)}{q(x)}$，其中$q(x)$为任意一组相互独立且同分布的随机变量 { X i } \lbrace X_i\rbrace {Xi​}，$X_i$在$[a,b]$上服从分布律$q(x)$，那么$q(x)$是随机变量$X$的PDF（或PMF），$g^{\ast }(x)$是关于$x$的函数，根据LOTUS（law of the unconscious statistician）那么

命题（随机变量的函数的期望规则）设随机变量$X$的PMF为$p_x(x_k)$，又设$g(x)$是$x$的函数，则$g(X)$的期望
$$E[g(x)]=\sum _{k}g(x_k)p_x(x_k)$$
推论1 期望和方差的性质
$$E[aX+b]=aE[x]+b$$
$$D[aX+b]=a^{2}D[x]$$
推论2 方差的计算公式
$$D[X]=E[X^2]-E[X{]}^2$$
推论3 矩的计算公式
$$E[(X-E[X]{)}^n]=\sum _k(x-E[X]{)}^n{p}_x(x_k)$$

PDF：概率密度函数（probability density function）, 在数学中，连续型随机变量的概率密度函数（在不至于混淆时可以简称为密度函数）是一个描述这个随机变量的输出值，在某个确定的取值点附近的可能性的函数。
PMF : 概率质量函数（probability mass function), 在概率论中，概率质量函数是离散随机变量在各特定取值上的概率。
CDF : 累积分布函数 (cumulative distribution function)，又叫分布函数，是概率密度函数的积分，能完整描述一个实随机变量X的概率分布。

$$E[g^{\ast }(x)]={\int }_a^b{g}^{\ast }(x)q(x)\mathrm{d}x={\int }_a^{b}f(x)\mathrm{d}x=I$$
由强大数定理
$$P(\underset{N\to +\infty }{\lim }\frac{1}{N}\sum _{i=1}^n{g}^{\ast }(X_i)=I)=1$$
若选
$$\overline{I}=\frac{1}{N}\sum _{i=1}^n{g}^{\ast }(X_i)$$
则$\overline{I}$依概率1收敛到I。平均值法就用$\overline{I}$作为I的近似值。

---

上面是理论部分，下面是实现的过程
假设要计算的积分有如下形式
$$I={\int }_a^{b}f(x)\mathrm{d}x$$
$f(x)$在$[a,b]$内可积，任意选择一个有简便办法可以进行抽样的概率密度函数$q(x)$，使其满足以下条件

1. 当$f(x)̸=0$时，$q(x)̸=0$，$(a\le x\le b)$
2. ${\int }_a^{b}q(x)\mathrm{d}x=1$

如果记
$$g^{\ast }(x)=\{\begin{array}{ll}\frac{f(x)}{q(x)},& \text{q(x)̸=0}\\ 0,& \text{q(x)=0}\end{array}$$
那么原积分式可以写成
$$I={\int }_a^b{g}^{\ast }(x)q(x)\mathrm{d}x$$
因而求积分的步骤是：

1. 产生服从分布律$q(x)$的随机变量$X_i(i=1,2,\cdots ,N)$
2. 计算均值
   $$\overline{I}=\frac{1}{N}\sum _{i=1}^n{g}^{\ast }(X_i)$$
   并用它作为I的近似值，即$I\approx \overline{I}$。

如果$a,b$为有限值，那么$q(x)$可取作为均匀分布：
$$q(x)=\{\begin{array}{ll}\frac{1}{b-a}& ,\text{q(x)̸=0}\\ 0& ,\text{otherwise}\end{array}$$
此时原来的积分式变为
$$I=(b-a){\int }_a^{b}g(x)\frac{1}{b-a}\mathrm{d}x$$
具体步骤如下：

1. 产生服从分布律$q(x)$的随机变量$X_i(i=1,2,\cdots ,N)$
2. 计算均值
   $$\overline{I}=\frac{b-a}{N}\sum _{i=1}^N{g}^{\ast }(X_i)$$
   并用它作为I的近似值，即$I\approx \overline{I}$。

自此，结束对蒙特卡洛数值积分的理解

---

---

2. 均匀分布，Box-Muller 变换

均匀分布很好理解，一般可以表示为X~Uniform[a,b]，计算机中的random方法就可以生成[0,1]区间内的随机数，这些随机数总体上服从均匀分布，Box-Muller 变换是通过均匀分布生成正态分布的随机变量。

Box-Muller 变换
如果随机变量 U1,U2 独立且U1,U2∼Uniform[0,1]
$Z_0=\sqrt{-2\ln U_1}\cos (2\pi U_2)$
$Z_1=\sqrt{-2\ln U_1}\sin (2\pi U_2)$
则$Z_0,Z_1$服从标准正态分布，即$Z_0$，$Z_1$都服从标准正态分布，$(Z_0,Z_1{)}^T$服从二元标准正态分布

[注]Box-Muller 变换的推导证明过程

PS:刚看完上面蒙特卡洛积分的推导过程，先歇一歇，看看代码去，明天补上0 0

[注]随机数生成算法

先来看一下对随机数的定义

在计算物理学中，随机数被准确地分成了三类：真随机数、准随机数、伪随机数。其中书上的定义为
1）真随机数：产生的数不可预计，也不可能重复产生两个相同的真随机数序列。真随机数只能通过某些随机的物理过程来产生，如放射性衰变、电子设备的热噪声等。
2）准随机数：其随机数序列不具备随机性质，仅仅是用它来处理问题能够得到正确的结果。
3）伪随机数：通过某种数学公式或者算法产生的数值序列。虽然在数学意义上伪随机数是不随机的，但是如果能够通过统计检验，可以当成真随机数使用。

PRNG(Pseudo-random number generator)，伪随机数发生器

1. 梅森旋转算法
2. 线性同余法

QRNG(Quasi-random number generator)，准随机数发生器

1. Halton序列
2. Sobol序列
3. Latin超立方体序列

伪随机数和准随机数的区别，伪随机数在小样本的情况下均匀分布的差异性很大，同样，维数灾难（curse of dimensionality）是另一个造成差异性的原因。准随机数就是生成一个不具备随机性质但是满足均匀分布的样本，用准随机数的蒙特卡罗方法叫做拟蒙特卡罗，在tensorflow的tfp.mcmc.sample_halton_sequence方法中，警告当样本维数不确定的时候，准随机数是完全确定的并且通常具有显着的负自相关。

Warning: The sequence elements take values only between 0 and 1. Care must be taken to appropriately transform the domain of a function if it differs from the unit cube before evaluating integrals using Halton samples. It is also important to remember that quasi-random numbers without randomization are not a replacement for pseudo-random numbers in every context. Quasi random numbers are completely deterministic and typically have significant negative autocorrelation unless randomization is used.

3. Monte Carlo principle（蒙特卡洛原则）

Monte Carlo principle就是上文中对蒙特卡洛数值积分的过程。

4. 接受-拒绝抽样（Acceptance-Rejection sampling)

4、5的内容可以在我的这篇博客中找到（拒绝采样、重要性采样）
但是，在高维的情况下，Rejection Sampling 会出现两个问题，第一是合适的 q 分布比较难以找到，第二是很难确定一个合理的 k 值。这两个问题会导致拒绝率很高，无用计算增加。

5. 重要性抽样(Importance sampling)

同样，可惜的是，在高维空间里找到一个这样合适的 q 非常难。即使有 Adaptive importance sampling 和 Sampling-Importance-Resampling(SIR) 的出现，要找到一个同时满足容易采样并且与原样本近似的提议分布，往往是不存在的！

6. 马尔科夫链，马尔科夫稳态

马尔科夫链的剋荣可以在我这篇博客中找到马尔科夫决策过程
在马尔科夫的转移概率矩阵中，任意初始概率分布开始都会收敛到稳定的结果，这个就是马尔科夫收敛定理。

马氏链的收敛定理
如果一个非周期的马尔科夫链具有状态转移矩阵$P$,且他的任意两个状态是连通的，那么${\lim }_{N\to +\infty }{P}_{ij}^n$存在且与$i$无关，记${\lim }_{N\to +\infty }{P}_{ij}^n=\pi (j)$，我们有

1. $$\underset{N\to +\infty }{\lim }{P}^n=\left[\begin{array}{ccccc}\pi (1)& \pi (2)& \dots & \pi (j)& \dots \\ \pi (1)& \pi (2)& \dots & \pi (j)& \dots \\ \dots & \dots & \dots & \dots & \dots \\ \pi (1)& \pi (2)& \dots & \pi (j)& \dots \\ \dots & \dots & \dots & \dots & \dots \end{array}\right]$$
2. $\pi (j)={\sum }_{i=0}^{\infty }\pi (i)P_{ij}$
3. $\pi$ 是 $\pi P=\pi$的唯一非负解，其中
   $$\pi =\left[\begin{array}{ccccc}\pi (1)& \pi (2)& \dots & \pi (j)& \dots \end{array}\right]，\sum _{i=0}^{\infty }\pi (i)=1$$
   $\pi$称为马尔科夫链的平稳分布

[注1]要求图是联通的（没有孤立点），同时不存在一个联通的子图是没有对外的出边的（就像黑洞一样）。
这个马氏链的收敛定理非常重要，所有的 MCMC(Markov Chain Monte Carlo) 方法都是以这个定理作为理论基础的。
[注2]马尔可夫链的极限分布性质决定了MCMC是无偏估计（unbiased estimation），即采样数趋于无穷时会得到求解目标数学期望的真实值，这将MCMC与其替代方法，例如变分贝叶斯估计（variational Bayesian inference）相区分，后者计算量通常小于MCMC但无法得到无偏估计。

7. MCMC——Metropolis-Hasting算法

对于给定的概率分布p(x),我们希望能有便捷的方式生成它对应的样本。由于马氏链能收敛到平稳分布， 于是一个很的漂亮想法是：如果我们能构造一个转移矩阵为P的马氏链，使得该马氏链的平稳分布恰好是p(x), 那么我们从任何一个初始状态x0出发沿着马氏链转移, 得到一个转移序列 x0,x1,x2,⋯xn,xn+1⋯,， 如果马氏链在第n步已经收敛了，于是我们就得到了 π(x) 的样本xn,xn+1⋯。
这个绝妙的想法在1953年被 Metropolis想到了，为了研究粒子系统的平稳性质， Metropolis 考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法，即Metropolis算法，并在最早的计算机上编程实现。Metropolis 算法是首个普适的采样方法，并启发了一系列 MCMC方法，所以人们把它视为随机模拟技术腾飞的起点。 Metropolis的这篇论文被收录在《统计学中的重大突破》中， Metropolis算法也被遴选为二十世纪的十个最重要的算法之一。
我们接下来介绍的MCMC 算法是 Metropolis 算法的一个改进变种，即常用的 Metropolis-Hastings 算法。由上一节的例子和定理我们看到了，马氏链的收敛性质主要由转移矩阵P 决定, 所以基于马氏链做采样的关键问题是如何构造转移矩阵P,使得平稳分布恰好是我们要的分布p(x)。如何能做到这一点呢？我们主要使用如下的定理。

定理：[细致平稳条件]
如果非周期的马氏链的转移矩阵$P$和分布$\pi (x)$满足
$$\pi (i)P_{ij}=\pi (j)P_{ji}，\text{for all i,j}$$
则$\pi (x)$是马氏链的平稳分布，上式被成为细致平稳条件（detailed balance condition）

假设我们已经有一个状态转移矩阵Q马氏链$(q(i,j))$表示从状态$i$转移到状态$j$的概率，也可以写成$q(j|i)$或者$q(i\to j)$，显然，通常情况下
$$p(i)q(i,j)̸=p(j)q(j,i)$$
也就是细致平稳条件不成立，所以$p(x)$不太可能是这个马氏链的平稳分布。我们可否对马氏链做一个改造，使得细致平稳条件成立呢？譬如，我们引入一个$\alpha (i,j)$，我们希望
$$p(i)q(i,j)\alpha (i,j)=p(j)q(j,i)\alpha (j,i)$$
取什么样的$\alpha (i,j)$以上等式能成立呢？最简单的，按照对称矩阵，我们可以取
$$\alpha (i,j)=p(j)q(j,i)$$
$$\alpha (j,i)=p(i)q(i,j)$$
于是，$p(i)q(i,j)\alpha (i,j)=p(j)q(j,i)\alpha (j,i)$就成立了，所以有
$$p(i)q(i,j)\alpha (i,j)=p(j)q(j,i)\alpha (j,i)$$
其中$$Q\prime (i,j)=q(i,j)\alpha (i,j)，Q\prime (j,i)=q(j,i)\alpha (j,i)$$
于是我们把原来具有转移矩阵$Q$的一个很普通的马氏链，改造为了具有转移矩阵$Q\prime$的马氏链，而$Q\prime$恰好满足细致平稳条件，由此马氏链$Q\prime$的平稳分布就是$p(x)$。
在改造$Q$的过程中引入的$a(i,j)$称为接受率，物理意义可以理解为在原来的马氏链上，从状态$i$以$q(i,j)$的概率转跳到状态$j$的时候，我们以$\alpha (i,j)$的概率接受这个转移，于是得到新的马氏链$Q\prime$的转移概率为$q(i,j)\alpha (i,j)$。
终于，通过上述的推导，对Metropolis-Hasting算法的跳转接受率有了正确的认识。

8. MCMC——Gibbs Sampling算法