0. Main Takeaway
- Sec 1. 介绍了蒙特卡洛采样 (MC) 和拒绝-接受采样.
- Sec 2. 介绍了马尔可夫链(MC)和基于马尔可夫链的采样.
- Sec 3. 介绍了马尔可夫链蒙特卡洛采样(MCMC)和它的一个改进版本Metropolis-Hastings (MH).
- Sec 4. 介绍了Gibbs采样, 它用来处理MCMC采样在高维空间低效的问题.
Note:
- 本文章主要基于刘建平的博客整理写成, 补充了一些细节.
1. Monte Carlo (MC) 采样
核心: $\mathbb{E}_{x \sim p(x)}[\cdots] = \frac{1}{n} \sum\limits_{i=1}^n[\cdots]$. 故我们可以用数值采样的方法近似期望, 而期望又可以表示为求和/积分的形式, 故提供了一种使用数值采样计算求和/积分的方法.
Monte Carlo是一类使用随机模拟去解决问题的方法的统称, 例如, 当我们希望估计某些不太好求解的积分问题:
我们可以只在$[a,b]$间随机选取一个点$x_0$, 然后将$f(x_0)*(b-a)$作为阴影部分面积的近似. 我们也可以更精细一点, 在$[a,b]$间随机采$N$个点$\left \{ x_i\right \}$, 然后将$(b-a) * \frac{\sum\limits_{i=1}^N x_i}{n}$的值作为面积的近似. 我们可以这么做的理论依据是: 若$x\sim U[a,b]$, 则$p(x) = \frac{1}{b-a}$, 则:
$$ \int_x f(x) = \int_x \frac{f(x)}{p(x)} p(x) = \mathbb{E}_{x\sim p(x)}[\frac{f(x)}{p(x)}] = \mathbb{E}_{x\sim p(x)}[(b-a) f(x)] \approx \frac{1}{n}\sum_{i=1}^N (b-a) f(x_i) $$
故我们实际上不用局限于均匀分布, 对于任意分布$p(x)$, 我们都可以使用$\frac{1}{n}\sum\limits_{i=1}^N \frac{f(x_i)}{p(x_i)} , \text{ where } x_i \sim p(x)$来近似求解$\int_x f(x)$.
1.1. Box-Muller
即便我们知道概率密度函数$p(x)$, 如何从$p(x)$中采样仍然是一个non-trivial的问题 (特别是从高维空间中), 使用Probability Integral Theorem是一个选择, 但是它需要我们求出累计分布函数同时求逆, 这仍是十分困难的. 为此一些trick常常需要被使用.
Theorem: Let $X$ be a continuous random variable with c.d.f. $F(x)$, then $F(X) \sim \mathcal{U}(0,1)$
Proof: Let $Y = F(X)$ and has a c.d.f. $\mathcal{G}(y)$, then:
$$ \mathcal{G}(y) = P(Y \leq y) = P(F(X) \leq y) = P(x \leq F^{-1}(y)) = F(F^{-1}(y)) = y $$
Then: Let $U \sim \mathcal{U}(0,1)$, given an arbitrary distribution with c.d.f. $F(x)$, then assign $X = F^{-1}(U)$, $X$ will have the distribution $F(x)$.
Box-Muller变换:
Input: $U_1,U_2 \sim U[0,1]$
Output: $X,Y \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$
$$ X = \cos(2\pi U_1)\sqrt{-2\ln U_2}\\ Y = \sin(2\pi U_1)\sqrt{-2\ln U_2} $$
Proof:
假设$X,Y \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$, 则$p(x,y) = \frac{1}{2\pi} \exp(- \frac{x^2 + y^2}{2})$, 若我们记$X=R\cos \Theta, Y = R \sin \Theta$, 由概率密度换元可得:
$$ p(r,\theta) = \frac{1}{2\pi}\exp (-\frac{r^2}{2})\lvert\det \frac{\partial (x,y)}{\partial (r,\theta)} \rvert = \frac{r}{2\pi}\exp (-\frac{r^2}{2}) $$
则$p(\theta) = \int_0^\infty p(r,\theta)\mathrm{d}r = \frac{1}{2\pi}$, 可知$\Theta \sim U[0,2\pi]$.
而$F(R\leq r) = \int_0^r \int_0^{2\pi}p(r,\theta) \mathrm{d}\theta \mathrm{d}r = 1 - \exp (-\frac{r^2}{2})$, 故由Probability Integral Theorem, $R= F^{-1}(U) =\sqrt {-2\ln (1-U)}$, 其中$U\sim U[0,1]$.
故我们如果取$U_1 = \frac{\Theta}{2\pi}$, $U_2 = 1-U$, 则我们有$U_1, U_2 \sim U[0,1]$, 且$X = R\cos\Theta =\cos(2\pi U_1)\sqrt{-2\ln U_2}, Y = R\sin\Theta = \sin(2\pi U_1)\sqrt{-2\ln U_2}$, 且$X,Y \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$.
其他一些常见的连续分布, 比如t分布, F分布, Beta分布, Gamma分布等, 都可以通过类似的方式从$U[0,1]$得到的采样样本转化得到. 在python的numpy, scikit-learn等类库中, 都有生成这些常用分布样本的函数可以使用.
1.2. Rejection Sampling
假设我们有一个简单易采样的分布$q(z)$, 我们希望从复杂的目标分布$\tilde p(z)$中采样 (我们记归一化后的分布为$p(z) = \frac{\tilde p(z)}{\int_z \tilde p(z)} = \frac{\tilde p(z)}{Z_p}$).
拒绝接受采样的工作过程如下:
首先选取$k$, 使得$\forall z, kq(z)\geq \tilde p(z)$.
Sample $z_0 \sim q(z)$.
以$\frac{\tilde p(z)}{kq(z)}$的概率接受这个样本, 否则拒绝.
- Sample $u \sim U[0,kq(z_0)]$, 若$u > \tilde p(z_0)$则拒绝这个样本 (落在图中灰色部分), 反之接受.
首先我们证明接受拒绝采样可以做到从$p(z)$中采样:
我们记二值随机变量$A$, $A= 1$表示接受, 反之拒绝, $p(A=1) = \int_z p(A=1|z) q(z) = \int_z\frac{\tilde p(z)}{kq(z)}q(z)= \frac{Z_p}{k}$则我们有:
- $$ p(z_0|A=1) = \frac{p(A=1|z_0)p(z_0)}{p(A=1)} = \frac{\frac{\tilde p(z_0)}{kq(z_0)}q(z_0)}{p(A=1)} = \frac{\tilde p(z_0)}{Z_p} = p(z_0) $$
接下来我们讨论$k$应该如何取值:
- 由$p(A=1)$可知, 接受率和$k$成反比, 故为了接受率更高, $k$应该越小越好.
Summary
作为MC方法中的一种, 接受-拒绝采样可以在一定程度上帮我们处理目标分布过于复杂难以采样的问题, 但是它仍然存在很多缺陷:
对于某些多维分布, 条件概率分布$p(X_i|X_{-i})$可能容易得到, 但是联合分布$p(X)$很难求.
对于高维度问题存在挑战:
- 假设$X\in R^{1000}$, $p(X) = \mathcal{N}(0, \Sigma_1), q(X) = \mathcal{N}(0,\Sigma_2)$, 且$\Sigma_2 = 1.001 \Sigma_1$.
- 则为了求$k$, 考虑零点处$\frac{p(0)}{q(0)} = \sqrt \frac{|\Sigma_2|}{|\Sigma_1|} = 1.001^{500}\approx 496,984,196,731,226,689,628.694$ (Curse of dimensionality), 会导致接受率极低.
总的来说接受-拒绝采样只在$q(z)$能较好拟合$\tilde p(z)$的情况下work的比较好.
2. Markov Chain (MC)
2.1. Markov Chain and Stationary Distribution
虽然接下来的推导都是在离散情况下, 但是可以被扩展到连续情况下.
对于随机序列$\left \{ X_i\right \}$, 如果对于任意$X_i$, 我们有
$$ p(X_i| X_{i-1}, \cdots, X_1) = p(X_i|X_{i-1}) $$
我们称其为马尔可夫链. 同时我们记转移矩阵$P_{ij} = P(i,j) = p(x_j |x_i)$, 在$t$时刻$X_t$取$x_i$的概率为: $\pi_{t}(i) = p(X_t = x_i)$.
我们记录 $\pi_t = \left \{ \pi_t(x_1), \cdots, \pi_t(x_i)\cdots\right \}$为当前链条在$t$时刻的概率分布, 则我们易得$\pi_{t+1} = \pi_t P$. 我们记当前马氏链的极限分布为$\lim\limits_{t\rightarrow\infty}\pi_t = \pi^{*}$ (如果存在), 我们即当前马氏链的平稳分布为 $\pi$ (如果存在):
$$ \pi = \pi P $$
由马氏链的性质我们可得:
对于非周期遍历链, 马氏链的极限分布存在且等于平稳分布.
- 非周期: 对于任意状态$x_i$, 记$d$为$\left \{ n|n\geq 1, p^n_{ii}>0\right \}$的最大公约数, 若$d=1$, 则该马氏链为非周期的 ($p^n_{ii}$是从$x_i$出发, 经过$n$步恰好回到$x_i$的概率). 直觉上理解, 对于周期马氏链会存在循环, 自然就不会有平稳分布.
- 遍历链: 正常返非周期链被称为遍历链条. 正常返指的是任意状态, 在无限长的马氏链中都会被访问无限次. 故遍历链可以简单认为就是任意状态之间两两可达. (不会有状态在$t\rightarrow \infty$的情况下被无法访问).
2.2. 基于Markov Chain采样
对于任意给定分布$\pi$, 如果我们能找到对应的(非周期遍历)马氏链 (对应的状态转移矩阵$P$), 使得$\pi$是该马氏链的平稳分布,那我们就可以通过这条马氏链采样得到服从分布$\pi$的样本.
假设$n$轮后马氏链收敛到平稳分布, 即:
$$ \pi_n = \pi_{n+1} = \cdots = \pi $$
则我们从$\pi_n$开始采样就可以得到符合我们目标分布$\pi$的样本了. 具体操作如下:
输入马尔科夫链状态转移矩阵 $P$ , 设定状态转移次数阈值 $n_{1}$ , 需要的样本个数 $n_{2}$.
从任意简单概率分布采样得到初始状态值 $x_{0}$
for $t=0$ to $n_{1}+n_{2}-1$ :
- 从条件概率分布 $P\left(x \mid x_{t}\right)$ 中采样得到样本 $x_{t+1}$
样本集 $\left(x_{n_{1}}, x_{n_{1}+1}, \ldots, x_{n_{1}+n_{2}-1}\right)$ 即为我们需要的平稳分布对应的样本集.
3. Markov Chain Monte Carlo (MCMC) and Metropolis-Hastings (MH)
3.1. Detailed Balance
在基于MC采样中, 最大的问题是给定分布$\pi$, 如何找到对应的马氏链 ($P$)使之满足$P$的平稳分布是$\pi$. MCMC用一种巧妙的迂回的方式 (类似接受-拒绝采样) 解决了这个问题, 同时MH是一种提升MCMC采样效率的简单改进.
在MCMC中, 主要解决的是, 给定分布$\pi$, 如何找到对应的马氏链 ($P$)使之满足$P$的平稳分布是$\pi$. 首先我们定义马氏链的细致平稳条件(detailed balance), 我们称$\pi$满足细致平稳条件, 如果它满足:
$$ \forall x_i, x_j, \pi(x_i)P(i,j) = x_j P(j,i) $$
不难证明, 如果$\pi$满足了细致平稳条件, 则它是马氏链的平稳分布:
$$ \sum_{i=1}^{\infty} \pi(x_i) P(i, j)=\sum_{i=1}^{\infty} \pi(x_j) P(j, i)=\pi(x_j) \sum_{i=1}^{\infty} P(j, i)=\pi(x_j) $$
即:
$$ \pi = \pi P $$
因此现在问题变成了, 给定分布$\pi$, 如何寻找转移方程$Q$, 满足 $\pi(x_i)Q(i,j) = \pi(x_j) Q(j,i)$.
3.2. MCMC
对于一般的$Q(i,j)$, 细致平稳条件是不满足的:
$$ \pi(x_i)Q(i,j) \neq \pi(x_j) Q(j,i) $$
如果我们对上式做一个改造, 引入一个$\alpha(i,j)$来使上式取等, 即:
$$ \pi(x_i)Q(i,j)\alpha(i,j) = \pi(x_j) Q(j,i)\alpha(j,i)\\ \text{where } \alpha(i,j) = \pi(j) Q(j,i) $$
这样我们就得到了满足细致平稳条件的转移矩阵$P$, 满足:
$$ P(i,j) =Q(i,j)\alpha(i,j) $$
同时$P$对应的平稳分布就是目标分布$\pi$.
这里我们的$P$是由某一个马氏链转移矩阵$Q(i,j)$乘以$\alpha(i,j)$得到, 其中$\alpha(i,j)$我们一般称为接受率, 它的取值在$[0,1]$之间, 可以理解为是一个概率值, 代表了我们接受转移的概率. 它和第1节中提到的拒绝-接受采样的思路类似, 前者是通过拒绝-接受概率拟合一个复杂分布, 这里是通过拒绝-接受概率得到一个满足细致平稳条件的转移矩阵.
注意, 对于$P(i,j) =Q(i,j)\alpha(i,j)$, $\sum\limits_{j} P(i,j)$不一定等于$1$, 这是正常的, 因为我们在使用MCMC做采样的时候, 下一步采样是有概率被拒绝的, 所以你可以在马氏链中额外加一个指向自己的自环, 概率值是$1- \sum\limits_{j} P(i,j)$.
我们得到MCMC的采样过程:
输入我们任意选定的马尔科夫链状态转移矩阵 $Q$ , 平稳分布 $\pi(x)$ , 设定状态转移次数阈值 $n_{1}$ , 需要的样本个数 $n_{2}$
从任意简单概率分布采样得到初始状态值 $x_{0}$.
for $t=0$ to $n_{1}+n_{2}-1$ :
从条件概率分布 $Q\left(x \mid x_{t}\right)$ 中采样得到样本 $x_{*}$.
以概率$\alpha\left(x_{t}, x_{*}\right)$接受样本$x_{*}$:
- 从均匀分布采样 $u\sim U [0,1]$
- 如果 $u<\alpha\left(x_{t}, x_{*}\right)$, $=\pi\left(x_{*}\right)$, $Q\left(x_{*}, x_{t}\right)$ 则接受转移 $x_{t} \rightarrow x_{*}$ , 即 $x_{t+1}=x_{*}$
- 否则不接受转移, 即 $x_{t+1}=x_{t}$
样本集 $\left(x_{n_{1}}, x_{n_{1}+1}, \ldots, x_{n_{1}+n_{2}-1}\right)$ 即为我们需要的平稳分布对应的样本集。
3.3. Metropolis-Hastings Sampling
Metropolis-Hastings采样简称M-H采样，这个算法首先由Metropolis提出，被Hastings改进，因此被称之为Metropolis-Hastings采样.
在上述MCMC采样中, 一个问题是, $\alpha(i,j)$作为接受概率, 如果$\alpha(i,j)$过低, 会导致MCMC采样的效率很低下. 回顾平稳条件:
$$ \pi(x_i)Q(i,j)\alpha(i,j) = \pi(x_j) Q(j,i)\alpha(j,i) $$
采样低下的根本原因是$\alpha(i,j)$太小了, 我们假设有:
$$ \pi(x_i)Q(i,j)\times 0.1 = \pi(x_j) Q(j,i)\times 0.2 $$
如果我们能同时把等式两边扩大$5$倍, 可以看到细致平稳条件仍然满足:
$$ \pi(x_i)Q(i,j)\times 0.5 = \pi(x_j) Q(j,i)\times 1 $$
因此如果我们令 $\alpha(i,j) = \min\left \{ \frac{\alpha(i,j)}{\alpha(j,i)}, 1\right \} = \min\left \{ \frac{\pi(x_j)Q(j,i)}{\pi(x_i)Q(i,j)},1\right \}$, 就得到了M-H采样, 同时易证细致平稳条件仍然满足:
$$ \pi(x_i)Q(i,j)\times \min\left \{ \frac{\alpha(i,j)}{\alpha(j,i)}, 1\right \}= \pi(x_j) Q(j,i)\times\min\left \{ \frac{\alpha(j,i)}{\alpha(i,j)}, 1\right \} $$
由此M-H采样过程为:
输入我们任意选定的马尔科夫链状态转移矩阵 $Q$ , 平稳分布 $\pi(x)$ , 设定状态转移次数阈值 $n_{1}$ , 需要的样本个数 $n_{2}$
从任意简单概率分布采样得到初始状态值 $x_{0}$.
for $t=0$ to $n_{1}+n_{2}-1$ :
从条件概率分布 $Q\left(x \mid x_{t}\right)$ 中采样得到样本 $x_{*}$.
以概率$\min\left \{ \frac{\alpha(x_t,x^{*})}{\alpha(x^{*},x_t)}, 1\right \}$接受样本$x_{*}$:
- 从均匀分布采样 $u\sim U [0,1]$
- 如果 $u< \min\left \{ \frac{\alpha(x_t,x^{*})}{\alpha(x^{*},x_t)}, 1\right \}$ , 则接受转移 $x_{t} \rightarrow x_{*}$ , 即 $x_{t+1}=x_{*}$
- 否则不接受转移, 即 $x_{t+1}=x_{t}$
如果我们选择的转移概率矩阵$Q$是对称的, 即满足$Q(i,j)=Q(j,i)$, 此时$\alpha(i,j)$的表示可以进一步简化为:
$$ \alpha(i,j) = \min\left \{ \frac{\pi(x_j)}{\pi(x_i)},1\right \} $$
在PRML中指出, 通过修改接受率$\alpha(i,j)$, M-H算法能在很大程度上提升MCMC的采样效率.
4. Gibbs Sampling
在第3节中, MCMC (M-H)算法已经可以做到从任意样本中采样的问题, 但是它仍然存在两个问题:
- 是需要计算接受率, 在高维时计算量大, 并且由于接受率的原因导致算法收敛时间变长.
- 对于高维数据, 往往数据的条件概率分布易得, 而联合概率分布不易得.
4.1. 重新寻找细致平稳条件
在3.2. MCMC中, 我们通过引入接受概率$\alpha(i,j)$来满足平稳分布条件, 这里我们尝试换一个思路.
考虑二维随机变量$(X,Y)$, 我们假设它们服从联合分布$\pi(x,y)$. 我们考虑两个第一维特征相同的两个点$A(x^1, y^1)$和$B(x^1, y^2)$, 我们不难发现:
$$ \pi(x^1, y^1) \pi(y^2|x^1) = \pi(y^1|x^1)\pi(x^1)\pi(y^2|x^1)\\ \pi(x^1, y^2) \pi(y^1|x^1)= \pi(y^2|x^1)\pi(x^1)\pi(y^1|x^1) $$
即$\pi(x^1, y^1) \pi(y^2|x^1) = \pi(x^1, y^2) \pi(y^1|x^1)$.
换言之, 在$X=x^1$这条直线上, 如果我们取$\pi(y^2|x^1)$作为从$A$到$B$的状态转移概率$P(A,B)$, 则在这条直线上的任意两点满足细致平稳条件. (保持$Y= y^1$可以得到类似的结果)
基于这个发现, 我们可以构造满足细致平稳条件的状态转移矩阵$P$:
$$ \left \{ \begin{aligned} P(A,B) &= \pi(Y^{(B)}|X), \text{ if }X^{(A)} = X^{(B)} =X\\ P(A,C) &= \pi(X^{(C)}|Y), \text{ if }Y^{(A)} = Y^{(C)} =Y\\ P(A,D) &= 0, \text{ otherwise} \end{aligned} \right . $$
由此我们得到了$P$, 它对应的平稳分布就是$\pi(x,y)$, 我们给出二维Gibbs采样的流程:
- 输入平稳分布 $\pi\left(x, y\right)$ , 设定状态转移次数阈值 $n_{1}$ , 需要的样本个数 $n_{2}$.
- 随机初始化初始状态值 $x^{(0)}$ 和 $y^{(0)}$
- for $t=0$ to $n_{1}+n_{2}-1$ :
- 从条件概率分布 $P\left(y \mid x^{(t)}\right)$ 中采样得到样本 $y^{t+1}$
- 从条件概率分布 $P\left(x \mid y^{(t+1)}\right)$ 中采样得到样本 $x^{t+1}$
样本集 $\left\{\left(x^{\left(n_{1}\right)}, y^{\left(n_{1}\right)}\right),\left(x^{\left(n_{1}+1\right)}, y^{\left(n_{1}+1\right)}\right), \ldots,\left(x^{\left(n_{1}+n_{2}-1\right)},y^{\left(n_{1}+n_{2}-1\right)}\right)\right\}$ 即为我们需要的平稳分布对应的样本集.
整个采样过程中，我们通过轮换坐标轴，采样的过程为:
$$ \left(x_{1}^{(1)}, x_{2}^{(1)}\right) \rightarrow\left(x_{1}^{(1)}, x_{2}^{(2)}\right) \rightarrow\left(x_{1}^{(2)}, x_{2}^{(2)}\right) \rightarrow \ldots \rightarrow\left(x_{1}^{\left(n_{1}+n_{2}-1\right)}, x_{2}^{\left(n_{1}+n_{2}-1\right)}\right) $$
用下图可以很直观的看出, 采样是在两个坐标轴上不停的轮换的. (坐标轴轮换不是必须的, 我们也可以每次随机选择一个坐标轴进行采样. 不过常用的Gibbs采样的实现都是基于坐标轴轮换的.)
4.2. 多维Gibbs采样
Gibbs采样至少需要两个维度 (1维可以使用M-H), 于是它可以被推广到高维. 在高维度中, 我们选取的马氏链的状态转移概率为$\pi(x_i|x_1,x_2,…,x_{i-1},x_{i+1},…,x_n) = \pi(x_i|x_{-i})$, 即我们固定$n-1$个坐标轴, 在剩余的一个坐标轴上进行移动(采样). 具体采样过程和上述二维中的类似.