马尔科夫蒙特卡罗方法

什么是采样

我们知道了一个变量的分布，要生成一批样本服从这个分布，这个过程就叫采样。比如在信号处理领域，采样是将信号从连续时间域上的模拟信号转换到离散时间域上的离散信号的过程。

为什么要采样

有许多原因需要我们从某个分布中采样。比如获取离散信号，比如数据挖掘中进行模型匹配，比如机器学习中估计某个函数在某个分布上的期望，比如以较小的代价近似许多项的和或者某个积分。
在许多情况下，采样的目标就是抽样，比如对大数据集中获取合适的子集。

采样的方法

蒙特卡罗采样

蒙特卡罗方法主要是找到满足分布的随机变量，对该随机变量随意取值生成若干样本。当样本足够多时，样本分布近似于目标分布。
当我们无法精确计算和或积分时，可以采用蒙特卡罗方法来近似。这种想法是把和或积分视作某个分布下的期望，然后通过估计对应的平均值来近似这个期望。令
$$s=\sum _{x}p(x)f(x)=E_p[f(x)]$$
或者
$$s=\int p(x)f(x)dx=E_p[f(x)]$$
其中$p(x)$为概率分布函数(求和)或者概率密度函数(积分)
我们可以通过从$p$中抽取n个样本$x^{(1)},...,x^{(n)}$来近似s并得到一个经验平均值
$$s_n=\frac{1}{n}\sum _{i=1}^{n}f(x^{(i)})$$
显然这是一个无偏估计，且由大数定律，$\underset{n\to \infty }{\lim }{s}_n=s$.
由中心极限定理有$s_n\sim \mathcal{N}(s,\frac{Var[f(x)]}{n})$，由此可以估计$s_n$的置信区间。

重要采样

蒙特卡罗采样中假设了可以从基准分布$p(x)$中轻易采样，当无法轻易采样时，可以换一个能够容易采样的分布，如上例中，$p(x)f(x)=q(x)\frac{p(x)f(x)}{q(x)}$，由此可以从$q$分布中采样，然后估计$\frac{pf}{q}$的期望。这样的转化显然不唯一，因此我们可以寻找最优的分布$q^{\ast }$。这样的分布可以容易的推导出来。令
$$s_q=\frac{1}{n}\sum _{i=1,x^{(i)}\sim q}^n\frac{p(x^{(i)})f(x^{(i)})}{q(x^{(i)})}$$
计算得：
$$\mathbb{E}[s_q]=s,Var[s_q]=Var[\frac{p(x)f(x)}{q(x)}]/n$$
可以发现估计的期望与$q$分布无关。我们希望方差取最小值，则$q$需要满足：
$$q^{\ast }(x)=\frac{p(x)|f(x)|}{Z}$$
这里$Z$是归一化常数，此时$Var[s_q]=0$.

一个好的$q$分布的选择可以显著地提高蒙特卡罗估计的效率，而一个糟糕的$q$分布则会使效率更糟糕。q分布通常会取一些简单常用的分布使得我们能够从$q$分布中容易的采样，但当$x$是高维数据时，$q$分布的简单性使得它很难和$p$或$p|f|$相匹配。

马尔科夫链蒙特卡罗方法（MCMC）

马尔科夫蒙特卡罗方法的基本思想是构造一个马氏链 { X n : n ⩾ 0 } \{X_n:n\geqslant0\} {Xn​:n⩾0}，使目标分布作为该马氏链的不变分布，再利用马氏链的极限理论（遍历定理），用$X_n$代替目标分布。

马氏过程相关概念和定理

定义1 转移矩阵
设状态空间$E$上的矩阵$P=(p(i,j):i,j\in E)$为转移矩阵，满足一下性质：
(1) 对任何$i,j\in E$，有$p(i,j)⩾0$
(2) 对任何$i\in E$，有$\sum _{j\in E}p(i,j)=1$

定义2 马氏链
对于随机过程 X = { X 0 , X 1 , X 2 , . . . , X n , . . . } X=\{X_0,X_1,X_2,...,X_n,...\} X={X0​,X1​,X2​,...,Xn​,...}。我们称$X$是以$P$为转移矩阵的马氏链，当且仅当$P(X_{n+1}=j|X_0=i_0,X_1=i_1,...,X_n=i_n)=p(i_n,j)$
由此可以得到$P(X_{n+1}=j|X_0=i_0,X_1=i_1,...,X_n=i_n)=p(i_n,j)=P(X_{n+1}=j|X_n=i_n)$
特别地，利用转移矩阵和初始分布可以给出马氏链的有限维分布的表示。

定义3 平稳分布（不变分布）
对于以$P$为转移矩阵的马氏链$X$，如果概率分布$\eta =({\eta }_i:i\in E)$满足$\sum _{i\in E}{\eta }_i{p}_{i,j}={\eta }_j,j\in E$，则称$\eta$为$P$或$X$的平稳分布（不变分布）。
定义4 可逆分布
对于以$P$为转移矩阵的马氏链$X$，如果概率分布$\eta =({\eta }_i:i\in E)$满足${\eta }_i{p}_{i,j}={\eta }_j{p}_{j,i},i,j\in E$，则称$\eta$为$P$或$X$的可逆分布。
显然可逆分布一定是平稳分布。

如果马氏链X以$P$为转移矩阵，以平稳分布$\eta$为初始分布，则称马氏链$X$是平稳的。
如果马氏链X以$P$为转移矩阵，以可逆分布$\eta$为初始分布，则称马氏链$X$是可逆的，也叫细致平衡的。

定理1（弱遍历定理）
设转移矩阵为$P$，令$L^{(n)}=\frac{1}{n}\sum _{n=0}^n{P}^{(n)}$。则$\underset{n\to \infty }{\lim }{L}^{(n)}=L$存在且满足$PL=LP=L=L^2$

定理2
设$P$是马氏链的转移矩阵，假定$P$不可约正常返，则$P$有唯一不变分布，并且转移概率的平均极限分布是马尔可夫链的平稳分布。

事实上
$$\underset{n\to \infty }{\lim }{L}^{(n)}=L=\left(\begin{array}{c}\pi \\ \pi \\ .\\ .\\ .\end{array}\right)$$
$\pi$为行向量且任意元素都大于0，此时$\pi$为$P$的唯一不变分布。

定理3
若$P$不可约，正常返，非周期，则有$\underset{n\to \infty }{\lim }=P^{(n)}=L$

MCMC算法

有了这些理论就可以给出马尔科夫蒙特卡罗法的采样过程了。先假设存在转移矩阵$P$，目标分布就是$P$的不变分布。
先从初始分布从采一个样本$x^{(0)}$，这个初始分布可以是状态空间上的任意分布。接着从利用转移矩阵得到在$x^{(0)}$条件下的分布$\mathbf{P}(\cdot |X_0=x^{(0)})$，在此分布下采样得到$x^{(1)}$，不断重复上述过程，可以发现每个时刻在这个马尔可夫链上进行随机游走一次，就可以得到一个样本。

# 算法1：生成一条马氏链轨迹
Proceduce MCMC_Sample(
	E 			#状态空间
	P_0 		#初始分布
	P	 		#马尔科夫转移模型
	n			#执行步数
)
1 Sample x^(0)~P_0
2 for t=1,2,...,n
3     Sample x^(t)~P(·|x^(t-1))
4 return x^(0),x^(1),...,x^(n)

在算法执行$m(<n)$步后，可以认为得到的样本集合 { x ( m + 1 ) , x ( m + 2 ) , . . . x ( n ) } \{x^{(m+1)},x^{(m+2)},...x^{(n)}\} {x(m+1),x(m+2),...x(n)}为所需要的样本。前$m$步称为混合期。
现在我们可以得到对$s$的估计$s_M=\frac{1}{n-m}\sum _{i=m+1}^{n}f(x^{(i)})$

Metropolis-Hastings算法

在MCMC算法中，我们假设了已经有这样的转移矩阵$P$，但事实上这样的矩阵很难直接找到。而Metropolis-Hastings算法就是用一种简单的方法构造出来转移矩阵$P$，从而能够应用MCMC算法采样。

注意到细致平衡的分布也是不变分布，由此可以构造一个符合条件的可逆过程的转移矩阵。设目标分布为$\pi$，将目标分布设为转移矩阵$P$的不变分布。首先令$Q$是任意一个取正整数的特定不可约马尔科夫转移概率矩阵，以$q(i,j)$表示$Q$的$i$行$j$列的元素。设目标分布为$\pi$现按照细致平衡的条件，对任意$i\ne j$，有
$$\pi (i)q(i,j)a(i,j)=\pi (j)q(j,i)a(j,i)$$
其中$a(i,j)$为配比出来的系数，于是可令$a(i,j)=\min (\frac{\pi (i)q(i,j)}{\pi (j)q(j,i)},1)$
所以可设$p(i,j)=q(i,j)a(i,j)$，为保持转移矩阵的性质，得到$p(i,i)=q(i,i)+\sum _{k\ne i}q(i,k)(1-a(i,k))$，最后设$P=(p(i,j))$，则$P$即为所需的转移矩阵。

# 算法2：M-H算法
Proceduce MCMC_Sample(E, P_0, Q, a, n)
1 Sample x^(0)~P_0
2 for t=1,2,...,n
3     Sample x*~Q(·|x^(t-1))
4     Sample u~U[0,1]
5     if u<a(x^(t-1),x*) then
6         x^(t) = x*
7     else
8         x(t) = x(t-1)
9 return x^(0),x^(1),...,x^(n)

MCMC方法的问题

1. 采集的相邻样本可能含有相关性
   一种解决办法是每隔n个样本返回一个样本，但这样不仅提高了采样时的计算代价，而且也没有完全解除掉相关性；另一种解决办法是同时并行多个马氏链，每隔马氏链只采集一个样本。这两种方法是两个极端，更多从业者选择马氏链和小批量中的样本数接近，然后再用这些马氏链抽取所需要的样本。
2. 无法预测混合期
   检测一个马尔科夫链是否达到了平衡很困难，只能粗略的估计这段时间是否是足够的。由此可见，此算法的关键就是能否针对具体的$\pi$找到一个能够快速收敛的马氏链，收敛速度的快慢决定了算法的优劣。
   由不变分布方程$\pi P=\pi$以及不变分布的唯一性知，$P$的特征值中有且只有一个$1$，且$\pi$是对应的唯一特征向量。注意到$P$是不可约的，则由Perron-Frobenius定理可以保证这个矩阵最大特征值为1.这就导致了其他特征值都随时间变化趋向于0，那么第二大特征值将决定了收敛速度，也就是混合期的大小。
   注： 在M-H算法中，预设分布的标准差也影响着混合时间。如果标准差过小，则采集的样本多集中在概率大的周围；若标准差过大，则被拒绝的概率将变大，混合时间将变长。
   

吉布斯采样（Gibbs）

在计算过程中常常会遇到复杂的联合分布，这增加了计算的复杂度。事实上不同分量之间可能存在相关性，利用概率图模型，将联合分布分解成多个条件分布，可以减少状态空间的参数。在采样过程中，将高维总体的采样，化为一系列一维分布的取样，这就是Gibbs采样法的要义。

# 算法3：Gibbs算法
Proceduce MCMC_Sample(E, P_0, P, T)
1 Sample x^(0)~P_0
2 for t=1,2,...,T
3     Sample x_1^(t) ~ P(x_1|x_2^(t-1), x_3^(t-1),...,x_n^(t-1))
4     Sample x_2^(t) ~ P(x_2|x_1^(t), x_3^(t-1),...,x_n^(t-1))
5     ...
6     Sample x_j^(t) ~ P(x_j|x_1^(t),..., x_(j-1)^(t),x_(j+1)^(t-1),...,x_n^(t-1))
7     ...
8     Sample x_n^(t) ~ P(x_n|x_1^(t), x_2^(t),...,x_(n-1)^(t))
9 return x^(0),x^(1),...,x^(n)

参考文献

1. 古德费洛, 本吉奥, 库维尔, 赵申剑, 黎彧君, and 符天凡. 深度学习. 北京: 人民邮电出版社, 2017.
2. 钱敏平. 应用随机过程. 北京: 高等教育出版社, 2011.
3. 罗斯, and 龚光鲁. 应用随机过程 概率模型导论. 北京: 人民邮电出版社, 2011.
4. 龚光鲁. 应用随机过程教程 及在算法和智能计算中的随机模型. 北京: 清华大学出版社, 2004.
5. Andrieu, Christophe, De Freitas, Nando, Doucet, Arnaud, and Jordan, Michael I. “An Introduction to MCMC for Machine Learning.” Machine Learning 50.1 (2003): 5-43.
6. 科勒, 弗里德曼, 王飞跃, and 韩素青. 概率图模型 原理与技术. 北京: 清华大学出版社, 2015.
7. Perron, Oskar. “Zur Theorie Der Matrices.” Mathematische Annalen 64.2 (1907): 248-63.