蒙特卡洛采样_生信小课堂 | MCMC（三）：MCMC采样方法

上期我们介绍了蒙特卡洛模拟(Monte Carlo sampling)和马尔科夫链(Markov Chain)。基于马尔科夫链对复杂高维概率分布进行采样，即给定初始概率分布及其对应的状态转移矩阵，通过很多次抽样使得马尔科夫链收敛到平稳分布(Equilibrium distribution)，状态概率和转移概率收敛，这时我们便通过采样得到了用以描述采样事件分布的参数。那如何找到平稳分布对应的状态转移矩阵呢？这就需要马尔科夫链蒙特卡洛模拟(MCMC)算法采样。本期我们将对MCMC算法模型进行简单介绍。

要基于MCMC方法找到马尔科夫链的平稳分布及其状态转移矩阵，首先需要了解马尔科夫链的性质：细致平稳条件(detailed balance condition)。给定事件概率分布π(x)，马尔科夫链状态转移矩阵P，对于所有的x, y满足：

则称π(x)是状态转移矩阵P的平稳分布。即两状态各自的状态概率与转移到对方状态的转移概率相乘相等。基于马尔科夫链的细致平稳条件，只要找到可以使得概率分布π(x)满足细致平稳条件的矩阵P即可。

给定目标分布π(x)和随机选定的状态转移矩阵P，如何求得满足条件的目标状态转移矩阵Q呢？在完全随机抽样情况下，π(x)和任意选定矩阵P很难满足细致平稳条件：

这时我们可以引入接受率(Acceptance rate)的概念。在我们更新概率矩阵的过程中，使得接受抽样并更新矩阵的概率同抽样得到的转移概率成正比。引入接受率α(x,y),可使其满足细致平稳条件。即：

由上式可以看出，目标矩阵Q可以通过任意状态转移矩阵P乘以接受率α得到。接受率α(x,y)，取值范围[0,1]。也就是说，在状态y，依赖于随机选定状态转移矩阵P转移到状态x 时，以一定的接受率α(x,y)接受状态转移。

上述过程即MCMC采样的基本原理，但这个采样过程很难在实际采样过程中得到应用。因其接受率往往很小，导致大部分采样值都被拒绝，采样效率极低。所以，Metropolis-Hastings(M-H)采样和Gibbs采样，解决MCMC采样过程接受率低的问题。 

Metropolis-Hastings(M-H) 采样

M-H采样即满足细致平稳条件的前提下，同时成比例扩大α(x,y),α(x,y)，使其中一个值扩大到1，这样同比例扩大接受率，但确没有破坏细节平衡条件。进而提高了采样的接受率。如公式3，假设α(y,x)=1，那么接受率α(x,y)：

相较于MCMC采样方法，M-H算法通过提高接受率，使得采样效率提高。下图为基于M-H采样的分布与真实分布比较示例：

由此可见，低维数据内，采样分布与真实分布拟合度较高。但在高维数据内，由于接受率的存在，有部分采样样本还是会被拒绝。那么，能否改进算法使α=1，从而使采得的样本达到100%的接受。这就需要Gibbs采样，Gibbs采样仅对高维数据情况有效。

Gibbs采样

Gibbs采样模型优势在于应用到高维数据，通过轮换n维坐标数据采集样本。对于高维数据，N维概率分布π(x1，x2，…xn)，假设对xi维坐标上进行采样，基于条件概率固定N-1维坐标轴在xi坐标轴进行转移采样。Gibbs采样分别对N维坐标轮换采样，使得在每个坐标维度采样均满足细节平衡条件，保证每次采得样本均以100%概率接受。以N维数据为例进行简要介绍：

鉴于前面关于贝叶斯模型、蒙特卡洛模拟、马尔科夫链及MCMC算法背景的介绍，下期将介绍MCMC算法在体细胞突变检测中的应用。

推荐阅读

1	生信小课堂第一期：概率分布与假设检验基础
2	生信小课堂第二期：概率分布在NGS中的常见应用
3	生信小课堂第三期：通俗理解贝叶斯统计
4	生信小课堂第四‍‍期：MCMC(二)：蒙特卡洛抽样与马尔科夫链

泛生子(纳斯达克代码：GTH)是中国领先的癌症精准医疗公司，专注于癌症基因组学研究和应用，并致力依托先进的分子生物学及大数据分析能力改变癌症诊疗方式。泛生子已打造了全面的产品及服务管线，覆盖从癌症早筛到诊断及治疗建议，再到监测及预后管理的癌症全周期。

4000-996-336

www.genetronhealth.com