刘浚嘉：PR：机器人学的概率方法学习路径​zhuanlan.zhihu.com

PR 采样分章 第二节：马尔可夫链蒙特卡洛 MCMC

Paper | 1970

在许多实例中，我们希望采用蒙特卡罗方法，然而往往又不存在一种简单的方法可以直接从目标分布

中精确采样或者一个好的（方差较小的）重要采样分布 q(x)。

在这种情况下，为了从分布

中近似采样，我们引入了一种称为

马尔可夫链（Markov Chain）的数学工具。

为了避免处理棘手的后验分布计算，我们从分布中采样，并用这些样本来计算分布的统计参数。

马尔科夫链

Markov Chain最重要的两个特性：

1. 无后效性，
   ；
2. 平稳性，详见 马尔可夫链，作者给出了一个关于社会阶层的很好例子。简要来说，不论初始概率分布如何，转移矩阵相同的情况下，最终状态的概率分布都会趋于同一个稳定的概率分布。即马尔科夫链模型的状态转移矩阵收敛到的稳定概率分布与初始状态概率分布无关。

如果非周期（状态转化不是循环的，否则永远不收敛）马尔可夫链的转移矩阵

和某个分布

满足：

则

是马尔可夫链的平稳分布。 这被称作马尔可夫链的细致平稳条件

detailed balance condition ，其证明如下：

基于马尔科夫链采样

马尔科夫链的平稳性意味着：

如果可以得到了一个稳定概率分布对应的马尔科夫链模型的状态转移矩阵，则我们可以用任意的概率分布样本开始，带入马尔科夫链模型的状态转移矩阵，这样经过一些序列的转换，最终就可以得到符合对应稳定概率分布的样本。

假设我们任意初始的概率分布是

, 经过第 1 轮马尔科夫链状态转移后的概率分布是

，...... 第 i 轮的概率分布是

。假设经过 n 轮后马尔科夫链收敛到我们的平稳分布

，即：

对于每个分布

，我们有：

采样步骤：

1. 从任意简单概率分布（比如高斯）
   开始采样得到状态值
   ；
2. 基于条件概率分布
   采样状态值
   ，一直进行下去，直至第 n 次 产生平稳分布；
3. 此时的采样集
   符合目标的平稳分布，可以用来做蒙特卡洛模拟求和了。

那么MCMC到底是什么思路呢？

1. 建立一个马尔科夫链，使其收敛到平稳分布恰好为 p(x)；
2. 从马尔可夫链中模拟随机的状态序列，该序列足够长，能够（几乎）达到稳态；
3. 保留生成的一些状态作为样本

好像和刚刚没什么区别。是不是还是一脸懵，别急，慢慢来，我们一步步分析。

我们已知在有状态转移矩阵转移矩阵的前提下，可以得到对应平稳分布的采样，但是问题是怎么得到状态转移矩阵啊？？

根据马尔可夫链的细致平稳条件，只要我们找到了可以使概率分布

满足

的矩阵

即可。但这往往不容易。

MCMC采样

一般情况下，矩阵 Q 不满足细致平稳条件：

对上式进行改造，使其等式成立：

这个

又该怎么得到？？其实满足下两式即可：

即可得到分布

对应的马尔科夫链状态转移矩阵

:

这就很像上一章重要性采样的思想，其中

可以视为权重/接受率。

MCMC采样过程：

1. 选定任意的马尔科夫链状态转移矩阵 Q，平稳分布
   ，设定状态转移次数阈值
   ，需要的样本个数
   ；
2. 从任意简单概率分布开始，采样得到初始状态
   ；
3. for
   to
   :

   - 从条件概率分布

   中采样得到样本

   - 从均匀分布采样u∼uniform[0,1]
   - 如果

   则接受转移
   ，即
   - 否则不接受转移，即

样本集

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

Implementation

还是我们在

Machine-Learning-is-ALL-You-Need​github.com

中的作风，实现形式包括纯python from scratch 的 self-implement 形式 + 常用标准库形式。

首先是 From scratch，自己复现对原理的理解会更加透彻：

TODO

概率机器学习标准库 GPflow：

GP regression​gpflow.readthedocs.io

我们想从以

为参数的后验分布采样:

num_burnin_steps = ci_niter(300)  # n_1
num_samples = ci_niter(500)       # n_2# Note that here we need model.trainable_parameters, not trainable_variables - only parameters can have priors!
hmc_helper = gpflow.optimizers.SamplingHelper(model.log_posterior_density, model.trainable_parameters
)hmc = tfp.mcmc.HamiltonianMonteCarlo(target_log_prob_fn=hmc_helper.target_log_prob_fn, num_leapfrog_steps=10, step_size=0.01
)
adaptive_hmc = tfp.mcmc.SimpleStepSizeAdaptation(hmc, num_adaptation_steps=10, target_accept_prob=f64(0.75), adaptation_rate=0.1
)@tf.function
def run_chain_fn():return tfp.mcmc.sample_chain(num_results=num_samples,num_burnin_steps=num_burnin_steps,current_state=hmc_helper.current_state,kernel=adaptive_hmc,trace_fn=lambda _, pkr: pkr.inner_results.is_accepted,)samples, traces = run_chain_fn()
parameter_samples = hmc_helper.convert_to_constrained_values(samples)param_to_name = {param: name for name, param in gpflow.utilities.parameter_dict(model).items()}

贝叶斯推断标准库 PyMC3：

Getting started with PyMC3​docs.pymc.io

Bayesian Regression in PYMC3 using MCMC & Variational Inference

import pymc3 as pmwith pm.Model() as basic_model:alpha_prior = pm.HalfNormal('alpha', sd=2, shape=2)beta_prior = pm.Normal('beta', mu=0, sd=2, shape=2)sigma_prior = pm.HalfNormal('sigma', sd=2, shape=1)mu_likelihood = alpha_prior[cat_tensor] + beta_prior[cat_tensor] * x_tensory_likelihood = pm.Normal('y', mu=mu_likelihood, sd=sigma_prior, observed=y_tensor)# HMC: Hamiltonian Monte Carlo, default MCMC method in PYMC3hmc_trace = pm.sample(draws=5000, tune=1000, cores=2)# 画图 + 输出pm.traceplot(hmc_trace)pm.summary(hmc_trace)

各种MCMC的标准库汇总

pythonMCMC​github.com

关联阅读

刘浚嘉：PR Sampling Ⅰ: 蒙特卡洛采样、重要性采样及python实现​zhuanlan.zhihu.com

刘浚嘉：PR Sampling Ⅲ: M-H and Gibbs 采样​zhuanlan.zhihu.com

刘浚嘉：Reinforcement-Learning-in-Robotics 专栏目录​zhuanlan.zhihu.com

Reference

1. Bayesian inference problem, MCMC and variational inference
2. 采样方法(Sampling Method)
3. MCMC sampling for dummies
4. MCMC(一)蒙特卡罗方法
5. MCMC(二)马尔科夫链
6. MCMC(三)MCMC采样和M-H采样 - 刘建平Pinard - 博客园