MCMC

最近在处理文本主题分割的项目,感觉还是应该仔细整理、推导一遍马尔可夫理论相关的知识。该笔记主要摘自Dustin Stansbury 的博客。


蒙特卡洛方法


强大数定律:若有一系列独立同分布随机变量序列 ${X_{1},...,X_{n}}$ 且其期望 $E(X_{i}) = \mu < \infty$,则有:

$P(lim_{n \to \infty}\frac{\sum_{i=1}^{n} X_{i}}{n}=\mu)=1$

无意识统计学家法则,又称LOTUS:已知随机变量 $X$ 的密度函数 $f_{X}$,而 $g(X)$ 的密度函数较难获得,则可通过如下方法计算其期望:

$E(g(X))=\int_{-\infty}^{\infty}g(x)f_{X}(x)dx$

由上述两前提可知,若要求函数 $g(x)$ 的积分,采样蒙特卡洛的方法如下式:

$\int_{a}^{b}g(x)dx=\int_{a}^{b}\frac{g(x)}{f_{X}(x)}f_{X}(x)dx=E(\frac{g(x)}{f_{X}(x)}) \approx \frac{1}{N}\sum_{i=1}^{N}\frac{g(X_{i})}{f_{X}(X_{i})}$

对于 $f_{X}$ 的选取,除了尽量选择可以简便抽样的概率密度函数外,还应满足当 $g(x) \neq 0$ 时,$f_{X}(x) \neq 0, (a \leq x \leq b)$


由上述描述可知,通过蒙特卡洛方法处理问题时,需要给出符合 $f_{X}$ 分布的样本序列 ${x_{1},...,x_{N}}$,即该序列满足:

$\forall k,\frac{\left \| \{x_{i}|x_{i} = k\} \right \|}{\left \| \{x_{i}\} \right \|} \approx f_{X}(k)$

所以,蒙特卡洛方法的核心问题在于如何从一个分布上采样

马尔可夫蒙特卡洛方法,又称MCMC,就是专门来解决该类问题的方法,其主题思想为通过马尔可夫链来生成所需要的样本。


马尔可夫链


马尔可夫链是一个按顺序执行的随机过程

$x^{(0)} \to x^{(1)} \to ... \to x^{t} \to ...$

一个马尔可夫链具有以下三个特征:

  • 状态空间 $x$:该马尔可夫链上所有状态可取的值的集合
  • 转移操作 $p(x^{(t+1)}|x^{(t)})$:从状态 $x^{(t)}$ 转移到 $x^{(t+1)}$ 的概率
  • 初始条件分布 $\pi^{(0)}$:描述了 $t=0$ 时所有可能状态的概率

马尔可夫链通过从 $\pi^{(0)}$ 种采样获得初始的状态,然后通过转移操作变换到下一个状态。

若存在转移操作的值不再随着转移操作的进行而改变的性质时,则称该马尔可夫链是时间同质的。当 $t \to \infty$ 时,时间同质马尔可夫链将会达到一个链的平稳分布,即

$p(x^{(t+1)}|x^{(t)})=p(x^{(t)}|x^{(t-1)})$


若一时间同质马尔可夫链的状态空间由有限个不同值所构成,则其转移操作可被定义为矩阵 $P$,其中元素为

$p_{ij}=P(x(t+1)=j|x^{(t)}=i)$

若马尔可夫链的状态空间是连续的,比如 $x \in \mathbb{R^{N}}$,则此时转移操作不再是一个矩阵,而是一些实数上的连续函数。同样连续状态空间的马尔可夫链也存在平稳分布。


Metropolis采样


由上述讨论可知,若直接从目标概率分布 $p(x)$ 采样比较困难的话,可以考虑构造一个平稳分布是该目标分布的马尔可夫链,通过从初始分布中采样再经过足够次数的转移操作来获取 $p(x)$ 的采样,下面介绍Metropolis抽样

  • 从某种初始分布 $\pi^{(0)}$ 中随机采样出初始状态 $x^{(0)}$
  • 循环进行以下步骤直至获得足够的样本数量
    • 根据上一步的状态生成新的建议分布 $q(x|x^{(t-1)})$
    • 从该建议分布中抽样一个候选样本 $x^{\ast}$
    • 计算该候选样本作为候选状态的接受概率 $\alpha=min(1,\frac{p(x^{\ast})}{p(x^{(t-1)})})$
    • 从均匀分布 $Unif(0, 1)$ 中采样出数值 $u$,若 $u \leq \alpha$,则 $x^{(t)}=x^{\ast}$,否则 $x^{(t)}=x^{(t-1)}$

需要注意的是,Metropolis采样并不需要 $p(x)$ 被归一化,这点从 $\alpha$ 的表达式中可以看出。


再次回顾马尔可夫链在平稳分布时的性质可知,此时 $x^{(t)} \to x^{(t-1)}$ 和 $x^{(t-1)} \to x^{(t)}$ 的转移概率应该是相等的,我们称这种性质为可逆性细致平衡

若建议分布是对称的,则Metropolis采样具有可逆性。正态分布、柯西分布、学生分布和均匀分布作为建议分布都是对称的。但有时候目标分布并不存在合适的对称的建议分布,由此引出下节相对应的解决方案。


Metropolis-Hastings采样


为解决非对称目标分布时Metropolis采样无法选用具有细致平衡性质的建议分布的问题,Metropolis-Hastings采样引入了计算校正因子的步骤,具体如下:

  • 从某种初始分布 $\pi^{(0)}$ 中随机采样出初始状态 $x^{(0)}$
  • 循环进行以下步骤直至获得足够的样本数量
    • 根据上一步的状态生成新的建议分布 $q(x|x^{(t-1)})$
    • 从该建议分布中抽样一个候选样本 $x^{\ast}$
    • 计算该候选样本的校正因子 $c=\frac{q(x^{(t-1)}|x^{\ast})}{q(x^{\ast}|x^{(t-1)})}$
    • 计算该候选样本作为候选状态的接受概率 $\alpha=min(1,\frac{p(x^{\ast})}{p(x^{(t-1)})} \times c)$
    • 从均匀分布 $Unif(0, 1)$ 中采样出数值 $u$,若 $u \leq \alpha$,则 $x^{(t)}=x^{\ast}$,否则 $x^{(t)}=x^{(t-1)}$

当状态 $x$ 为多维时,称以上更新方式为分块更新


可以预见,随着状态维数的增加,候选样本被拒绝的概率大大增加,一种解决办法是依次单独更新状态中每一维的数据,称为分量更新,具体步骤如下:

  • 从某种初始分布 $\pi^{(0)}$ 中随机采样出初始状态 $\mathbb{x}^{(0)}$
  • 循环进行以下步骤直至获得足够的样本数量
    • 分布对每一维 $i = 1,...,D$ 进行更新,具体步骤如下:
    • 根据上一步的状态生成新的建议分布 $q(x_{i}|x_{i}^{(t-1)})$
    • 从该建议分布中抽样一个候选样本 $x_{i}^{\ast}$
    • 计算该候选样本的校正因子 $c=\frac{q(x_{i}^{(t-1)}|x_{i}^{\ast})}{q(x_{i}^{\ast}|x_{i}^{(t-1)})}$
    • 记 $\mathbb{x}_{o} = {x_{j}|j \neq i}$
    • 计算该候选样本作为候选状态的接受概率 $\alpha=min(1,\frac{p(x_{i}^{\ast},\mathbb{x}_{o}^{(t-1)})}{p(x^{(t-1)}, \mathbb{x}_{o}^{(t-1)})} \times c)$
    • 从均匀分布 $Unif(0, 1)$ 中采样出数值 $u$,若 $u \leq \alpha$,则 $x_{i}^{(t)}=x_{i}^{\ast}$,否则 $x_{i}^{(t)}=x_{i}^{(t-1)}$

Gibbs采样


分量更新虽然比分块更新要有效率,但仍会存在拒绝的现象。Gibbs采样沿用了分量更新的模式,但接受所有的候选样本,从而避免计算的浪费。

要使用Gibbs采样,状态 $\mathbb{x}$ 必须满足两个标准:

  • 条件概率分布 $p(x_{i}|\mathbb{x}_{o})$ 存在数学书的解析表达式
  • 条件概率分布可被采样

具体Gibbs采样的步骤如下:

  • 从某种初始分布 $\pi^{(0)}$ 中随机采样出初始状态 $\mathbb{x}^{(0)}$
  • 循环进行以下步骤直至获得足够的样本数量
    • 分布对每一维 $i = 1,...,D$,根据 $p(x^{(t)}_{i}|x_{1}^{(t)},...,x_{i-1}^{(t)},x_{i+1}^{(t-1)},...,x_{D}^{(t-1)})$ 采样 $x_{i}^{(t)}$

Gibbs采样常用于贝叶斯方法中,而在于一些无法得出条件概率解析表达式的情况下很难被运用。


哈密顿蒙特卡洛


哈密顿蒙特卡洛,又称混合蒙特卡洛,记为HMC,是一种借鉴物理学上哈密顿系统的思路来预测马尔可夫链下一个状态的方法。与Metropolis-Hastings采样和Gibbs采样相比,HMC能避免二者在高维模型中因随机游走而造成的寻找目标分布低效的问题。


哈密顿动力学是一种描述目标在系统中运动的方法,基本的变量有:

  • $t$:观测的时间。
  • $\mathbb{x}$:目标的位置,对应于统计学中的参数 $\theta$ 的值。
  • $\mathbb{p}$:目标的动量,可以视为人为添加的辅助变量。
  • $U(\mathbb{x}),K(\mathbb{p})$:势能和动能。
  • $H(\mathbb{x},\mathbb{p}) = U(\mathbb{X}) + K(\mathbb{p})$:总能量。

在哈密顿动力学,动能和势能间存在着一组转换关系式:

  • $\frac{\partial x_{i}}{\partial t}=\frac{\partial H}{\partial p_{i}}=\frac{\partial K(\mathbb{p})}{\partial p_{i}}$
  • $\frac{\partial p_{i}}{\partial t}=-\frac{\partial H}{\partial x_{i}}=-\frac{\partial U(\mathbb{x})}{\partial x_{i}}$

由上述关系式可知,若已知 $\frac{\partial U(\mathbb{x})}{\partial x_{i}}$ 和 $\frac{\partial K(\mathbb{p})}{\partial p_{i}}$ 的表达式以及初始值 $\mathbb{x}^{(t_{0})},\mathbb{p}^{(t_{0})}$,则可预测任意时刻 $t_{n}$ 目标所对应的位置和动量。

而为了在计算机中计算哈密顿方程式,需要将连续的等式采取近似的方式离散化。其主体思路是通过将一段时间分割成一系列间隔为 $\delta$ 的时间序列来达到离散化的目的。

xiehao

Read more posts by this author.