前言:本文主要是关于MCMC, 这是一类对概率进行采样的方法, 主要参考了Coursera上国立高等经济大学Advanced Machine Learning系列课程Course3: Bayesian Methods for Machine Learning Week4.
本文讲一讲Markov chain Monte Carlo(MCMC), MCMC是一类对概率\(q(x)\)进行采样的方法, 其本质是构建稳态分布恰为所需分布\(q(x)\)的马尔可夫链, 这里主要介绍两种MCMC方法: Gibbs采样和Metropolis-Hastings算法.
首先要明确的是我们期望抽样的目标分布通常是高维分布, 而对高维分布抽样的方法都基于一维分布的抽样, 因此接下来会介绍一些一维抽样方法. 接着我们会分别回顾MCMC的前半和后半, 即马尔可夫链和Monte Carlo蒙特卡洛的相关知识, 最后介绍Gibbs采样和Metropolis-Hastings算法, 并结合两者的优点给出一种能够进行并行计算的抽样方法.
准备知识
一维分布的抽样
一维离散分布的抽样:
- 假设一维离散分布的密度函数为\[q(x)=q(x_1,...,x_k)=\prod_{i=1}^kp_i^{x_i}, \sum_{i=1}^{k}p_k=1\]其中\(x_i\in\{0,1\}, \sum_{i=1}^{k}x_k=1\). 只需取\(y\sim U(0,1)\), 令\(x=e_k\)当\(y\in\left[\sum_{i=1}^{k-1}p_i, \sum_{i=1}^{k}p_i\right)\)
一维连续分布的抽样:
- 假设一维连续分布的密度函数为\(q(x)\), 其对应的分布函数为\(F(y)=P(x\le y)=\int_{-\infty}^{y}q(x)dx\). 若\(F^{-1}(y)\)易求, 只需取\(y\sim U(0,1)\), 令\(x=F^{-1}(y)\), 可以得到\(x\sim q(x)\). 证明如下\[P(x\le a)=P(F^{-1}(y)\le a)=P(y\le F(a))=F(a)\]对a求导可得\(P(x=a)=q(a)\)
- 一般\(F^{-1}(y)\)难求, 这里介绍一种拒绝采样方法:
- 取\(y\sim N(a,b)=r(x)\), 找到常数\(M>0\), 满足\(q(x)\le Mr(x)\)
- 取\(z\sim U(0,Mq(y))\), 令\(x=y\)当\(z\le q(y)\)时, 否则拒绝\(y\)
- 其他的采样方法如Adaptive rejection sampling、Importance sampling及Sampling-importance-resampling可参考Bishop的PRML的11.1
马氏链
马尔可夫过程, 又称马氏过程, 且当其时间状态离散时, 马尔可夫过程又称为马尔可夫链, 又称马氏链.
- 马氏链由结构\(G\)和转移概率矩阵\(M\)构成. 结构\(G\)是一个有向图, 每个节点对应一个状态;转移概率矩阵\(M\)定义在节点上, \(M_{ij}\)表示第\(i\)个状态转移到第\(j\)个状态的概率. 我们说由转移概率矩阵\(M\)可以构造一个马氏链, 若\(M\)满足:
- \(\sum_{j=1}M_{ij}=1, \forall i\), 即\(M\)的行和为1
- K-C方程(自然满足)
- 称\(\pi\)为稳态分布, 如果满足如下平稳条件\(\pi_j=\sum_i\pi_i M_{ij}\). 可以说, 一个马氏链从任意状态出发最终都能到达稳态分布, 若能构造稳态分布即为期望分布的马氏链, 我们的目的也就达成了.
p.s. 这里假定本文提到的马氏链均为时间离散且状态空间可数的时齐马氏链. 状态空间可数很好理解, 对于每一台机器来说都有一个机器精度的存在, 所以实数上取值集合实际上是可数的; 这里的时齐可以简单理解为不管从哪个时间点出发, 其转移概率矩阵都是一样的. KC方程说的是对任意的\(s\le r\le t\), 从时间\(s\)转移到时间\(t\)的转移概率矩阵等于\(s\)到\(r\)再到\(t\)的三维转移张量对\(r\)时刻的所有状态求和, 乍一看这是很自然的事情, 而也正因为时齐及状态空间可数, K-C方程才能得以自然满足.
蒙特卡洛
蒙特卡洛是通过重复随机抽样来获得数值结果的一种计算方法, 假定目标为获得\(\mathbb{E}_{q(x)}f(x)\)的估计值, 蒙特卡洛做的就是:
- \(x_{s}\sim q(x)\)
- \(\mathbb{E}_{q(x)}f(x)\approx \frac{1}{M}\sum_{s=1}^Mf(x_s)\)
下面举个例子来说明MCMC的用处, 回忆EM算法中的求期望的E步, 可见博文EM算法
- 在\(Q\)中找到\(q(T)\approx P(T\mid X, \theta)\)
- \(\mathbb{E}_{P(T\mid X, \theta)}\log P(X,T\mid\theta)\approx \mathbb{E}_{q(T)}\log P(X,T\mid\theta)\)
两者同样是对某期望的近似, 一个自然的想法是把MCMC替换掉EM算法中的E步, 即
- \(T_s\sim P(T\mid X, \theta)\)
- \(\mathbb{E}_{P(T\mid X, \theta)}\log P(X,T\mid\theta)\approx \sum_{s=1}^M\log P(X,T_s\mid\theta)\)
我们要做的就是从\(P(T\mid X, \theta)\)中抽样, 然后得到期望的估计值. 下面以EM算法为例介绍两种MCMC方法:Gibbs采样和Metropolis-Hastings算法
Gibbs采样
Gibbs采样的想法是化多维抽样为一维抽样, 假定采样目标分布为\(P(T\mid X, \theta)\), 假定\(T=\{T_1,...,T_d\}\), 第\(k\)次采样得到的数据为\(T^k=\{t_1^k,...,t_d^k\}\), 第\(k+1\)次采样策略如下\[\begin{aligned}t_1^{k+1}&\sim P(T_1\mid T_2=t_2^k, T_3=t_3^k,..., T_d=t_d^k, X, \theta) \\\ t_2^{k+1}&\sim P(T_2\mid T_1=t_1^{k+1}, T_3=t_3^k,..., T_d=t_d^k, X, \theta) \\\ &...\\\ t_d^{k+1}&\sim P(T_d\mid T_1=t_1^{k+1}, T_2=t_2^{k+1},..., T_{d-1}=t_{d-1}^{k+1}, X, \theta)\end{aligned}\]首先验证上述策略定义的矩阵\(M\)满足转移概率矩阵的条件, 即验证转移阵的行和均为1:\[\begin{aligned}\forall i, \sum_j M_{ij}&=\sum_jP(T_1=j_1\mid T_2=i_2, T_3=i_3,..., T_d=i_d, X, \theta) \\\ &\times P(T_2=j_2\mid T_1=j_1, T_3=i_3,..., T_d=i_d, X, \theta) \\\ &...\\\ &\times P(T_d=j_d\mid T_1=j_1, T_2=j_2,..., T_{d-1}=j_{d-1}, X, \theta)\\\ &=\sum_{j_{-d}}P(T_1=j_1\mid T_2=i_2, T_3=i_3,..., T_d=i_d, X, \theta) \\\ &\times P(T_2=j_2\mid T_1=j_1, T_3=i_3,..., T_d=i_d, X, \theta) \\\ &...\\\ &\times \sum_{j_{d}}P(T_d=j_d\mid T_1=j_1, T_2=j_2,..., T_{d-1}=j_{d-1}, X, \theta)\end{aligned}\]容易看到\(\sum_{j_{d}}P(T_d=j_d\mid T_1=j_1, T_2=j_2,..., T_{d-1}=j_{d-1}, X, \theta)=1\), 因而以此类推能得到结论, 所以可以说由上述策略定义的转移阵\(M\)能构造出一个概率空间及其上的马氏链.
下面证明我们的目标分布恰为该马氏链的稳态分布, 即验证平稳条件: \[P(T=j\mid X, \theta)=\sum_i P(T=i\mid X, \theta)M_{ij}\]式子太长这里就不再写出来了, 有兴趣的朋友可以尝试去验证一下:\[P(T_1=j_1, T_2=j_2,..., T_d=j_d\mid X, \theta)=\sum_i P(T_1=i_1, T_2=i_2, ..., T_d=i_d\mid X, \theta) M_{ij}\]你需要做的就是仿照上面验证行和为1, 将\(M_{ij}\)带入验证 ,很容易发现, 右边带有\(T_1=i_1\)的只有一项, 因此可以先对\(i\)的第一个分量\(i_1\)求和, 剩下的以此类推, 不在赘述.
最后由马氏链的性质就可以推出, 在转移足够长的次数\(K\)后, 其后的数据\(T_K, T_{K+1},...\)恰来自我们目标分布\(P(T\mid X, \theta)\)
Metropolis-Hastings算法
Metropolis-Hastings算法的想法是先让抽样按照随机选取的\(Q\)非常随意的进行, 然后再根据稳态分布应该满足的条件去制定一个拒绝策略\(A\). 假定第\(k\)次采样数据为\(x_k\), 其第\(k+1\)次采样思路如下:
- \(x'\sim Q(x^k\to x)\)
- 以\(A(x^k\to x')\)的概率接受\(x'\), 令\(x^{k+1}=x'\);否则令\(x^{k+1}=x^k\)
写成数学表达式即为:\[\begin{aligned}M_{ij}&= Q_{ij}A_{ij}, j\neq i \\\ M_{ii}&= Q_{ii}A_{ii}+\sum_{j\neq i}Q_{ij}(1-A_{ij})\end{aligned}\]下面根据\(M\)应该满足的性质来构造拒绝策略\(A\), 首先应当满足\(\sum_j M_{ij}=1\) \[\begin{aligned}\sum_j M_{ij}&= \sum_{j\neq i}Q_{ij}A_{ij}+Q_{ii}A_{ii}+\sum_{j\neq i}Q_{ij}(1-A_{ij})\\\ &= \sum_{j\neq i}Q_{ij}+Q_{ii}A_{ii}\end{aligned}\]若\(A_{ii}=1\), 那么上式\(=1\), 即M为转移概率阵. \(\pi_j M_{ji}=\pi_i M_{ij}\)被称为是细平稳条件, 因为两边对\(i\)求和能推得平稳条件\[\sum_i\pi_i M_{ij}=\sum_i\pi_j M_{ji}=\pi_j\]
若直接假定细平稳条件成立, 当\(i\neq j\)时即为\[\pi_j Q_{ji}A_{ji}=\pi_i Q_{ij}A_{ij}\]变形得\(\frac{A_{ij}}{A_{ji}}=\frac{\pi_j Q_{ji}}{\pi_i Q_{ij}}, i\neq j\). 最后令\(A_{ij}=\min\{1, \frac{\pi_j Q_{ji}}{\pi_i Q_{ij}}\}, i\neq j\)即可. 事实上, \(A_{ii}=1\)恰能统一进上式故\[A_{ij}=\min\left\{1, \frac{\pi_j Q_{ji}}{\pi_i Q_{ij}}\right\}\]
举个例子, \(x'\sim Q(x\to x')=N(x,1)\), 此时拒绝策略为\[A(x\to x')=\min\left\{1, \frac{\pi(x') }{\pi(x)}\right\}\]
Gibbs+MH
若Metropolis-Hastings算法的Q直接选用稍作修改版的Gibbs采样策略, 这样的修改使得并行计算能得以进行\[\begin{aligned}t_1^{k+1}&\sim P(T_1\mid T_2=t_2^k, T_3=t_3^k,..., T_d=t_d^k, X, \theta) \\\ t_2^{k+1}&\sim P(T_2\mid T_1=t_1^{k}, T_3=t_3^k,..., T_d=t_d^k, X, \theta) \\\ &......\\\ t_d^{k+1}&\sim P(T_d\mid T_1=t_1^{k}, T_2=t_2^{k},..., T_{d-1}=t_{d-1}^{k}, X, \theta)\end{aligned}\]此时并不能保证此种策略能收敛到目标分布\(P(T\mid X, \theta)\), 但加上拒绝策略就能保证\[A(T\to T')=\min\left\{1, \frac{P(T'\mid X, \theta) Q(T'\to T)}{P(T\mid X, \theta) Q(T\to T')}\right\}\]