MCMC算法

马尔科夫链蒙特卡罗方法(Markov Chain Monte Carlo,以下简称MCMC),是用于得到服从某一特定分布的随机采样的方法。很多复杂算法求解的基础,都使用到了 MCMC 算法。

本文内容主要参考张勤老师的《动物遗传育种中的计算方法》和刘建平的博客MCMC(一)蒙特卡罗方法

随机数产生方法

服从某一特定分布的随机变量的抽样值称为随机数,我们通常采用数学方法来产生随机数。

这种方法是利用一定的数学递推公式来产生随机数,即在给定一个任意的初值后,就可由此递推公式产生任意多的随机数,这种方法非常容易在计算机上实现。这种方法最大的缺陷在于给定初值后,以后所有的随机数便被唯一地确定下来了,而且这些随机数还存在周期现象,即随机序列达到一定长度后会出现重复,因此严格来说这些随机数并不是真正的随机数,故一般称之为伪随机数 (pseudo random number)。但由于初值可随机确定,因此其后的一系列数也可看作是随机的。

[0,1]均匀随机数的产生

[0,1]均匀随机数是在[0,1]区间上均匀分布随机变量的抽样值,其是产生其他一切分布随机数的基础。也就是说,任何其他分布的随机数都可通过对[0,1]均匀随机数进行某种转换得到。

目前产生[0,1]均匀随机数的最广泛的方法是线性同余法,也称为线性同余发生器 (linear congruential generator) 。

线性同余法

Lehmer (1951) 提出这种方法,它利用数论中的同余运算原理产生随机数,其递推公式为

其中 是事先给定的参数,均为非负整数; 是在 [0,m] 区间取值的一个随机整数; 是在 [0,1] 区间取值的一个随机实数。

在给定一个任意非负整数的初值 后,由上式求出整数序列 和实数序列

,则称为乘同余法;若 ,则称为混同余法

其他分布随机数的产生

任何其他分布的随机数都可以通过[0,1]均匀进行某种转换得到,例如正态分布可以通过著名的 Box-Muller 变换得到,根据下面公式通过两个[0,1]均匀随机数,我们可以得到两个彼此独立的标准正态分布随机数(缺证明,下式中的 为两个独立的[0,1]均匀随机数)。

下面介绍两种基本方法,反函数法与接受-拒绝采样法。

反函数法

假设我们已知某个分布的累积分布函数 ,那么抽样值可以通过下式得到

其中, 为[0,1]均匀随机数。

证明:因为 服从[0,1]均匀分布,所以 ,因此,我们有

所以 的累积分布函数就是 ,因此 就是 的一个抽样值。

这种方法非常简单,但是要求随机变量累积分布函数的反函数存在,而且容易计算。在很多情况下,累积分布函数的反函数不存在,或者虽然存在但其计算却非常困难,此时这种方法并不适用。

接受-拒绝采样法

对于概率分布不常见的分布,假设我们已知其概率密度分布为 ,此时我们用一个程序可采样的分布 ,比如正态分布。具体过程如下:

  1. 首先我们找到一个常数 , 使得 恒成立,即 恒在 下方。

  2. 采样得到 的一个样本 ,然后从均匀分布 中采样得到一个值 。如果 ,则拒绝这一次抽样,否则则接受这次抽样。

  3. 重复以上过程,直到得到 个抽样值。

整个过程从直观理解来看,就是从 曲线下方的二维图形中随机抽取一个点,如果这个点正好落在 的曲线下方,则接受这个抽样,否则就拒绝这个抽样。

因此使用接受-拒绝采样法时,只有当 曲线图形越接近时,抽样接受率越高,抽样效率越高。

蒙特卡罗方法

蒙特卡罗方法 (Monte Carlo) 是指使用随机数来解决很多计算问题的方法。

总的来说,蒙特卡罗方法处理的问题可以分为两大类:

第一类是确定性的数学问题。如计算定积分、解线性和非线性方程组,解积分方程和某些偏积分方程等。在这类问题中,因不含时间因素,故常将相应的模拟模型称为静态模型。

第二类是随机性问题,例如中子在介质中的扩散问题,运筹学中的库存问题等。

用蒙特卡罗方法计算定积分

用蒙特卡罗方法计算定积分是其最典型的应用,这里有两种算法,一是随机投点法,二是平均值法。下面以单重积分为例,分别进行介绍。

随机投点法

假设所求定积分为

因此,如果我们向 (0,0) (0,1) (1,1) (1,0) 这四个点构成的单位正方形均匀地投点,落在曲线 下方的概率就是我们要求的定积分,证明如下:

因此,我们可以按照以下步骤计算这个积分:

  1. 独立地产生两个 [0,1]均匀随机数 ,则 就构成了在单位正方形中地一个均匀分布的随机点
  2. 检验该随机点是否落在曲线 下方

重复上面两个步骤 次,设有 次随机点落在曲线 下方,则当 足够大时,我们有

我们可以将上面的方法一般化,假设我们所求的定积分为

我们需要先对这个积分进行变换。由 可得 ,因此我们有

,有

其中 ,它满足 ,因此积分 $ \int_{0}^{1} g^{*}(z) d z$ 可以按照上面的方法进行计算。

平均值法

假设我们所求的定积分为

是 [a,b] 区间上的均匀分布随机变量,则

于是,我们有

因此,我们可以用以下方法计算积分的近似解

  1. 产生 [a,b] 区间上的均匀分布随机数
  2. 计算

足够大时,我们有

其实这种做法就是在 [a,b] 区间上随机采样 个点,用这 个点的函数值均值代表 [a,b] 区间上的所有函数值,即用一个矩形的面积代表原来的曲线下方的面积。

这里我们有一个隐藏的假设是 在 [a,b] 区间是均匀分布的,而大多数时候这个假设均不成立。此时我们进一步假设 在 [a,b] 区间上的概率密度函数为 ,此时我们有蒙特卡罗方法的一般形式 :

因此,平均值法计算积分的一般过程如下:

  1. 从概率密度函数为 的分布中生成随机数 $ x_{i}$ ,
  2. 计算

足够大时,我们有

可以看出当我们假设 在 [a,b] 区间满足均匀分布时,我们将 带入上式,得到

我们得到了均匀分布的估计公式,即均匀分布是上式的一个特列。

小结

蒙特卡罗方法的核心在于生成特定分布的随机数,常用的分布可以通过对于 [0,1] 均匀随机数进行转换得到。通过采用接受-拒绝采样,我们也可以得到一些不常见分布的随机数。

但是在很多时候,有一些分布我们还是很难得到其随机数,比如:

1)对于一些二维分布,这个函数本身计算困难,但是我们能得到其条件分布 ,这时我们无法用接受-拒绝采样得到其样本集。

2)对于一些高维的分布 ,我们要找到一个合适的 非常困难。

此时我们需要马尔科夫链的方法帮助我们进行采样。

马尔科夫链

马尔科夫链 (Markov chain) 是一种最简单,应用最多的随机过程。随机过程是依赖于时间的一组随机变量。马尔科夫链要求具备“无记忆”的性质:下一状态的概率分布只能由当前状态决定,在时间序列中它前面的事件均与之无关,表示为

我们将系统随着时间的变化而发生的状态改变称为系统状态的转移,将已知在时刻 系统处于状态 的条件下,在下一时刻处于状态 的概率称为转移概率,表示为 ,即

如果系统由状态 在下一时刻转移为状态 的概率与系统当前所处的时刻无关,则称其为齐次的 (homogeneous) 。对于齐次的马尔科夫链,我们可以将转移概率简单地表示为 。一般来说,我们涉及到的马尔科夫链均为齐次的,因此下文均默认为齐次的马尔科夫链。

此时我们将所有的转移概率放到一个矩阵中,称为转移概率矩阵,如下。

转移概率矩阵的每一行之和为1,即

举个例子,假设天气变化是一个随机过程,假设天气的状态空间为 (有雨,晴天,多云),再假设天气变化是一个马尔科夫链过程,即明天出现何种天气的概率仅与今天的天气有关。设状态转移矩阵为

其中,第一行表示如果今天有雨,则明天为有雨,晴天,多云的概率分别为 0.5,0.25 和 0.25 。

假设系统初始的概率分布向量为 ,那么第一天的天气概率分布向量为 ,第二天的分布向量为 ,依次类推,第 天的分布向量为 .

我们发现随着天数的增加, 会逐渐变成一个稳定不变的分布,此时我们称马尔科夫链收敛,收敛时系统状态空间的概率分布称为平稳分布 (stationary distribution),表示为

在本例中为 ,而且与初始值无关。进一步发现其原因是由于 之后就变成了一个稳定不变的矩阵,即

但是不是所有的马尔科夫链都能达到平稳分布。

性质

不可约:如果从马尔科夫链中的任一状态出发,都能以大于零的概率经过有限次转移变为其他任一状态,则称该马尔科夫链是不可约的 (irreducible) 。假设 为从状态 经过 k 次 () 转移变成状态 的概率,不可约是说对于任意的 ,必定存在一个正整数 ,使得 。换句话说,如果马尔科夫链是不可约的,那么它的任何两个状态是连通的(两个状态彼此均可以通过有限次转移变为另一种状态,记为 ),即可以经过一次或多次转移从某一种状态转移到任意一种状态。

非周期:如果所有满足 中没有大于1的公因子,则称状态 非周期的 (aperiodic)。举个例子,假设存在一个马尔科夫链,其状态空间为 (1,2,3), 而状态之间的转移只能是 1-2-3-1-2-3 ,即系统状态每经过三次转移就会发生重复,假设我们从状态 2 开始,则经过 3,6,9… 次转移之后会重复出现状态 2 ,因此我们说状态 2 是周期性的,其周期为3。由周期性的定义存在一下推论:

  1. 推论:若状态A与状态B连通,则A与B周期相同。
  2. 推论:若不可约的马尔科夫链有周期性状态A,则该马尔科夫链的所有状态为周期性状态

只有不可约非周期的马尔科夫链才能收敛于平衡状态,并且存在以下性质

  1. 平稳分布 的唯一非负解。

基于马尔科夫链采样

假如我们想得到某个马尔科夫链的平稳分布的随机数,如果我们得到了其转移概率矩阵,那么我们就可以轻松做到这一点。

我们可以从任意一个初始概率分布 ,假设我们认为经过 轮后马尔科夫链收敛于平稳分布,那么之后采样得到的随机数就是我们需要服从平稳分布的随机数,其过程如下:

  1. 输入转移概率矩阵 ,设定状态转移次数阈值为 ,需要的样本数为
  2. 从任意简单概率分布采样得到初始状态值
  3. 从条件概率分布 采用得到

最终,样本集合 ( ) 即为我们需要的满足平稳分布的随机数。

但是随意给定一个平稳分布 ,我们如何得到其相应的转移概率矩阵 呢?此时我们就需要将马尔科夫链和蒙特卡洛方法结合起来,即MCMC 方法。

MCMC 算法

在解决如果得到某个平稳分布相应的转移概率矩阵 这个问题之前,我们需要先看马尔科夫链的细致平稳条件

马尔科夫链的细致平稳条件

如果非周期马尔科夫链的转移概率矩阵 和概率分布 对于所有的 ,均满足

就是转移概率矩阵 的平稳分布。

证明如下,由细致平稳条件有

其中我们有

因此我们得到

因此, 就是转移概率矩阵 的平稳分布,证明完毕。

MCMC 采样

一般情况下,我们的目标平稳分布 和我们随机找的某一个马尔科夫链转移概率矩阵 不满足细致平稳条件,即

我们可以对上式做一个改造,使得细致平稳条件成立,方法式引入 ,使得

那么什么样的 可以使得上式成立呢?最简单的就是只要满足下面两个式子即可:

这样,我们就得到了我们的目标平稳分布 对应的转移概率矩阵 ,其元素满足

也就是说,我们的目标矩阵 可以通过任意一个转移概率矩阵 构建得到,这里的 我们一般称为接受率,这个思路很像上面的接受-拒绝采样。

MCMC 采样的步骤如下:

  1. 输入任意的转移概率矩阵 ,平稳分布 ,设定状态转移次数阈值 ,需要的样本个数
  2. 从任意简单概率分布采样得到初始状态值
    1. 从条件概率分布 采用得到
    2. 从[0,1]均匀分布采样
    3. 如果 ,则接受转移,即 ;否则不接受转移,即

最终,样本集合 ( ) 即为我们需要的满足平稳分布的随机数。

但是这个采样算法很难在实践中应用,因为这里的接受率 可能非常的小,这样导致大部分的采样值都被拒绝转移,采样效率很低。

我们采用M-H 采样方法可以提高采样效率。

M-H 采样

M-H采样是Metropolis-Hastings采样的简称,这个算法首先由Metropolis提出,被Hasting扩展,其解决了接受率过低的问题。

我们回到MCMC采样的细致平稳条件:

我们的问题是接受率 太小了,比如 ,此时如果等式的两边同时扩大 5 倍,此时接受率提高到了 0.5,但是细致平稳条件依然是满足的。

因此,我们的接受率可以改为

其实这种做法可以换种方式解释,我们回到一开始我们说 不满足细致平稳条件,即

这可以分为两种情况

对于第一种情况,我们可以在左侧乘以概率 ,使得细致平稳条件成立

对于第二种情况,我们可以在右侧乘以概率 ,使得细致平稳条件成立

第一种情况的接受率为 ,第二种情况的接受率为 ,因此我们得到

通过这个微小的改造,M-H采样算法过程如下

  1. 输入任意的转移概率矩阵 ,平稳分布 ,设定状态转移次数阈值 ,需要的样本个数
  2. 从任意简单概率分布采样得到初始状态值
    1. 从条件概率分布 采用得到
    2. 从[0,1]均匀分布采样
    3. 如果 ,则接受转移,即 ;否则不接受转移,即

最终,样本集合 ( ) 即为我们需要的满足平稳分布的随机数。

很多时候,我们选择的转移概率矩阵 都是对称矩阵,即 ,此时我们的接受率可以进一步简化为

小结

M-H算法有下面几个问题:

  1. 对于高维度数据,由于每一次抽样均需要计算 ,计算效率低
  2. 由于 ,仍然有可能拒绝一部分抽样。
  3. 有时,对于高纬度数据,我们可能不知道其联合分布,但是可以方便计算各个特征之间的条件概率分布,此时无法使用 M-H 抽样。

吉布斯抽样 (Gibbs) 解决了上面的几个问题,其接受概率为1,即不拒绝任何采样,而且其抽样过程中只计算条件分布。

Gibbs 抽样

二维抽样

我们回到细致平稳条件,我们有

我们假设我们的数据是二维数据,平稳分布为 ,观测第一个特征相同的两个点 ,根据概率论的知识,容易发现下面两个式子成立

由于两个等式的右手项相等,因此我们有

也就是说

联合上面的细致平稳条件,我们发现在 $ x_{1}=x_{1}^{(1)}$ 这条直线上,我们只需要选择 ,则任意两个点均满足细致平稳条件。

同样的道理,在在 $ x_{2}=x_{2}^{(1)}$ 这条直线上,我们只需要选择 ,则任意两个点也满足细致平稳条件。

基于上面的发现,我们可以构建 的转移概率矩阵为

我们可以轻松验证对于二维平面上的任意两点 ,细致平稳条件均成立,即:

于是这个二维空间的马尔科夫链将收敛到平稳分布

二维 Gibbs 抽样具体过程如下:

  1. 输入平稳分布 ,设定状态转移次数阈值 ,需要的样本个数
  2. 从任意简单概率分布采样得到初始状态值
    1. 从条件概率分布 采样得到
    2. 从条件概率分布 采样得到

最终,样本集合 即为我们需要的满足平稳分布的随机数。

整个采样过程中,我们不停轮换坐标轴,即

多维抽样

我们从二维扩展到多维,假设我们有一个 维的平稳分布 ,我们可以通过对 个坐标轴轮换采样,来得到新的样本。具体的抽样过程如下:

  1. 输入平稳分布 或者对应的所有特征的条件概率分布,设定状态转移次数阈值 ,需要的样本个数
  2. 随机初始化初始状态值
  3. for to :
    1. 从条件概率分布 中采样得到样本
    2. 从条件概率分布 中采样得到样本
    3. 从条件概率分布 中采样得到样本
    4. 从条件概率分布 中采样得到样本

样本集 即为我们需要的平稳分布对应的样本集。

小结

由于 Gibbs 抽样对于高维度特征具有优势,因此一般所说的 MCMC 采样都是用的 Gibbs 抽样。但是 Gibbs 抽样要求至少有两个维度,因此一维的数据采样只能采用 M-H 采样。

收敛性判断

在MCMC 算法中,无论是 M-H 算法还是 Gibbs 算法,我们均认为随着马尔科夫链的延长,它就越接近其平稳分布,我们称其为收敛。而收敛前的所有迭代是要被舍弃的,这些被舍弃的抽样称为 burn-in 期。但是我们如何判断一条链是否收敛呢?也就是说,上面的 应该设定为多少了?

这似乎没有一个统一的答案,目前常用的有两种方法。

图式法

我们将 MCMC 链上的轨迹点进行可视化,横坐标是迭代次数,纵坐标是抽样值。如果经过一段时间的剧烈波动之后逐渐趋于平稳,则意味着这条链可能收敛了。

更为合理的做法是,用不同的参数初值构建两个甚至多个独立的链,观察不同链的轨迹点,随着链的延长,不同链开始发生重叠,甚至几乎完全重合,此时可以认为收敛。

需要说明的是,图式法不能保证所观察的链一定收敛了,为稳妥起见, burn-in 期稍微设大一点较好。

Geweke 检验

首先假定一个burn-in 期,在其后的抽样值按一定间隔取2个样本,如前 10% 和 后 10% 的抽样值的平均数之间不存在显著差异,则可认为链收敛了。否则,需要延长burn-in 期,再重新进行检验。

参考文献

  1. 张勤. 动物遗传育种中的计算方法[M]. 科学出版社, 2007.

  2. MCMC(一)蒙特卡罗方法

  3. Brémaud P. Markov chains: Gibbs fields, Monte Carlo simulation, and queues[M]. Springer Science & Business Media, 2013.

  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2019-2022 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信