什么是采样

我们知道了一个变量的分布，要生成一批样本服从这个分布，这个过程就叫采样。比如在信号处理领域，采样是将信号从连续时间域上的模拟信号转换到离散时间域上的离散信号的过程。

为什么要采样

有许多原因需要我们从某个分布中采样。比如获取离散信号，比如数据挖掘中进行模型匹配，比如机器学习中估计某个函数在某个分布上的期望，比如以较小的代价近似许多项的和或者某个积分。
在许多情况下，采样的目标就是抽样，比如对大数据集中获取合适的子集。

采样的方法

蒙特卡罗采样

蒙特卡罗方法主要是找到满足分布的随机变量，对该随机变量随意取值生成若干样本。当样本足够多时，样本分布近似于目标分布。
当我们无法精确计算和或积分时，可以采用蒙特卡罗方法来近似。这种想法是把和或积分视作某个分布下的期望，然后通过估计对应的平均值来近似这个期望。令
s=∑xp(x)f(x)=Ep[f(x)]s=\sum\limits_xp(x)f(x)=E_p[f(x)]s=x∑p(x)f(x)=Ep[f(x)]
或者
s=∫p(x)f(x)dx=Ep[f(x)]s=\int p(x)f(x)dx=E_p[f(x)]s=∫p(x)f(x)dx=Ep[f(x)]
其中p(x)p(x)p(x)为概率分布函数(求和)或者概率密度函数(积分)
我们可以通过从ppp中抽取n个样本x(1),...,x(n)x^{(1)},...,x^{(n)}x(1),...,x(n)来近似s并得到一个经验平均值
sn=1n∑i=1nf(x(i))s_n=\frac{1}{n}\sum\limits_{i=1}^nf(x^{(i)})sn=n1i=1∑nf(x(i))
显然这是一个无偏估计，且由大数定律，lim⁡n→∞sn=s\lim\limits_{n\rightarrow\infty}s_n=sn→∞limsn=s.
由中心极限定理有sn∼N(s,Var[f(x)]n)s_n\sim\mathcal{N}(s,\frac{Var[f(x)]}{n})sn∼N(s,nVar[f(x)])，由此可以估计sns_nsn的置信区间。

重要采样

蒙特卡罗采样中假设了可以从基准分布p(x)p(x)p(x)中轻易采样，当无法轻易采样时，可以换一个能够容易采样的分布，如上例中，p(x)f(x)=q(x)p(x)f(x)q(x)p(x)f(x)=q(x)\frac{p(x)f(x)}{q(x)}p(x)f(x)=q(x)q(x)p(x)f(x)，由此可以从qqq分布中采样，然后估计pfq\frac{pf}{q}qpf的期望。这样的转化显然不唯一，因此我们可以寻找最优的分布q∗q^*q∗。这样的分布可以容易的推导出来。令
sq=1n∑i=1,x(i)∼qnp(x(i))f(x(i))q(x(i))s_q=\frac{1}{n}\sum\limits_{i=1,x^{(i)}\sim q}^n\frac{p(x^{(i)})f(x^{(i)})}{q(x^{(i)})}sq=n1i=1,x(i)∼q∑nq(x(i))p(x(i))f(x(i))
计算得：
E[sq]=s,Var[sq]=Var[p(x)f(x)q(x)]/n\mathbb{E}[s_q]=s,Var[s_q]=Var[\frac{p(x)f(x)}{q(x)}]/nE[sq]=s,Var[sq]=Var[q(x)p(x)f(x)]/n
可以发现估计的期望与qqq分布无关。我们希望方差取最小值，则qqq需要满足：
q∗(x)=p(x)∣f(x)∣Zq^*(x)=\frac{p(x)|f(x)|}{Z}q∗(x)=Zp(x)∣f(x)∣
这里ZZZ是归一化常数，此时Var[sq]=0Var[s_q]=0Var[sq]=0.

一个好的qqq分布的选择可以显著地提高蒙特卡罗估计的效率，而一个糟糕的qqq分布则会使效率更糟糕。q分布通常会取一些简单常用的分布使得我们能够从qqq分布中容易的采样，但当xxx是高维数据时，qqq分布的简单性使得它很难和ppp或p∣f∣p|f|p∣f∣相匹配。

马尔科夫链蒙特卡罗方法（MCMC）

马尔科夫蒙特卡罗方法的基本思想是构造一个马氏链{Xn:n⩾0}\{X_n:n\geqslant0\}{Xn:n⩾0}，使目标分布作为该马氏链的不变分布，再利用马氏链的极限理论（遍历定理），用XnX_nXn代替目标分布。

马氏过程相关概念和定理

定义1 转移矩阵
设状态空间EEE上的矩阵P=(p(i,j):i,j∈E)P=(p(i,j):i,j\in E)P=(p(i,j):i,j∈E)为转移矩阵，满足一下性质：
(1) 对任何i,j∈Ei,j\in Ei,j∈E，有p(i,j)⩾0p(i,j)\geqslant 0p(i,j)⩾0
(2) 对任何i∈Ei\in Ei∈E，有∑j∈Ep(i,j)=1\sum\limits_{j\in E}p(i,j)=1j∈E∑p(i,j)=1

定义2 马氏链
对于随机过程X={X0,X1,X2,...,Xn,...}X=\{X_0,X_1,X_2,...,X_n,...\}X={X0,X1,X2,...,Xn,...}。我们称XXX是以PPP为转移矩阵的马氏链，当且仅当P(Xn+1=j∣X0=i0,X1=i1,...,Xn=in)=p(in,j)P(X_{n+1}=j|X_0=i_0,X_1=i_1,...,X_n=i_n)=p(i_n,j)P(Xn+1=j∣X0=i0,X1=i1,...,Xn=in)=p(in,j)
由此可以得到P(Xn+1=j∣X0=i0,X1=i1,...,Xn=in)=p(in,j)=P(Xn+1=j∣Xn=in)P(X_{n+1}=j|X_0=i_0,X_1=i_1,...,X_n=i_n)=p(i_n,j)=P(X_{n+1}=j|X_n=i_n)P(Xn+1=j∣X0=i0,X1=i1,...,Xn=in)=p(in,j)=P(Xn+1=j∣Xn=in)
特别地，利用转移矩阵和初始分布可以给出马氏链的有限维分布的表示。

定义3 平稳分布（不变分布）
对于以PPP为转移矩阵的马氏链XXX，如果概率分布η=(ηi:i∈E)\eta=(\eta_i:i\in E)η=(ηi:i∈E)满足∑i∈Eηipi,j=ηj,j∈E\sum\limits_{i\in E}\eta_ip_{i,j}=\eta_j,j\in Ei∈E∑ηipi,j=ηj,j∈E，则称η\etaη为PPP或XXX的平稳分布（不变分布）。
定义4 可逆分布
对于以PPP为转移矩阵的马氏链XXX，如果概率分布η=(ηi:i∈E)\eta=(\eta_i:i\in E)η=(ηi:i∈E)满足ηipi,j=ηjpj,i,i,j∈E\eta_ip_{i,j}=\eta_jp_{j,i},i,j\in Eηipi,j=ηjpj,i,i,j∈E，则称η\etaη为PPP或XXX的可逆分布。
显然可逆分布一定是平稳分布。

如果马氏链X以PPP为转移矩阵，以平稳分布η\etaη为初始分布，则称马氏链XXX是平稳的。
如果马氏链X以PPP为转移矩阵，以可逆分布η\etaη为初始分布，则称马氏链XXX是可逆的，也叫细致平衡的。

定理1（弱遍历定理）
设转移矩阵为PPP，令L(n)=1n∑n=0nP(n)L^{(n)}=\frac{1}{n}\sum\limits_{n=0}^{n}P^{(n)}L(n)=n1n=0∑nP(n)。则lim⁡n→∞L(n)=L\lim\limits_{n\to\infty}L^{(n)}=Ln→∞limL(n)=L存在且满足PL=LP=L=L2PL=LP=L=L^2PL=LP=L=L2

定理2
设PPP是马氏链的转移矩阵，假定PPP不可约正常返，则PPP有唯一不变分布，并且转移概率的平均极限分布是马尔可夫链的平稳分布。

事实上
lim⁡n→∞L(n)=L=(ππ...)\lim\limits_{n\to\infty}L^{(n)}=L=\begin{pmatrix}\pi\\\pi\\.\\.\\.\end{pmatrix}n→∞limL(n)=L=⎝⎜⎜⎜⎜⎛ππ...⎠⎟⎟⎟⎟⎞
π\piπ为行向量且任意元素都大于0，此时π\piπ为PPP的唯一不变分布。

定理3
若PPP不可约，正常返，非周期，则有lim⁡n→∞=P(n)=L\lim\limits_{n\to\infty}=P^{(n)}=Ln→∞lim=P(n)=L

MCMC算法

有了这些理论就可以给出马尔科夫蒙特卡罗法的采样过程了。先假设存在转移矩阵PPP，目标分布就是PPP的不变分布。
先从初始分布从采一个样本x(0)x^{(0)}x(0)，这个初始分布可以是状态空间上的任意分布。接着从利用转移矩阵得到在x(0)x^{(0)}x(0)条件下的分布P(⋅∣X0=x(0))\mathbf P(\cdot|X_0=x^{(0)})P(⋅∣X0=x(0))，在此分布下采样得到x(1)x^{(1)}x(1)，不断重复上述过程，可以发现每个时刻在这个马尔可夫链上进行随机游走一次，就可以得到一个样本。

# 算法1：生成一条马氏链轨迹
Proceduce MCMC_Sample(E             #状态空间P_0        #初始分布P          #马尔科夫转移模型n          #执行步数
)
1 Sample x^(0)~P_0
2 for t=1,2,...,n
3     Sample x^(t)~P(·|x^(t-1))
4 return x^(0),x^(1),...,x^(n)

在算法执行m(<n)m(<n)m(<n)步后，可以认为得到的样本集合{x(m+1),x(m+2),...x(n)}\{x^{(m+1)},x^{(m+2)},...x^{(n)}\}{x(m+1),x(m+2),...x(n)}为所需要的样本。前mmm步称为混合期。
现在我们可以得到对sss的估计sM=1n−m∑i=m+1nf(x(i))s_M=\frac{1}{n-m}\sum\limits_{i=m+1}^nf(x^{(i)})sM=n−m1i=m+1∑nf(x(i))

Metropolis-Hastings算法

在MCMC算法中，我们假设了已经有这样的转移矩阵PPP，但事实上这样的矩阵很难直接找到。而Metropolis-Hastings算法就是用一种简单的方法构造出来转移矩阵PPP，从而能够应用MCMC算法采样。

注意到细致平衡的分布也是不变分布，由此可以构造一个符合条件的可逆过程的转移矩阵。设目标分布为π\piπ，将目标分布设为转移矩阵PPP的不变分布。首先令QQQ是任意一个取正整数的特定不可约马尔科夫转移概率矩阵，以q(i,j)q(i,j)q(i,j)表示QQQ的iii行jjj列的元素。设目标分布为π\piπ现按照细致平衡的条件，对任意i≠ji\neq ji=j，有
π(i)q(i,j)a(i,j)=π(j)q(j,i)a(j,i)\pi(i)q(i,j)a(i,j)=\pi(j)q(j,i)a(j,i)π(i)q(i,j)a(i,j)=π(j)q(j,i)a(j,i)
其中a(i,j)a(i,j)a(i,j)为配比出来的系数，于是可令a(i,j)=min⁡(π(i)q(i,j)π(j)q(j,i),1)a(i,j)=\min(\frac{\pi(i)q(i,j)}{\pi(j)q(j,i)},1)a(i,j)=min(π(j)q(j,i)π(i)q(i,j),1)
所以可设p(i,j)=q(i,j)a(i,j)p(i,j)=q(i,j)a(i,j)p(i,j)=q(i,j)a(i,j)，为保持转移矩阵的性质，得到p(i,i)=q(i,i)+∑k≠iq(i,k)(1−a(i,k))p(i,i)=q(i,i)+\sum\limits_{k\neq i}q(i,k)(1-a(i,k))p(i,i)=q(i,i)+k=i∑q(i,k)(1−a(i,k))，最后设P=(p(i,j))P=(p(i,j))P=(p(i,j))，则PPP即为所需的转移矩阵。

# 算法2：M-H算法
Proceduce MCMC_Sample(E, P_0, Q, a, n)
1 Sample x^(0)~P_0
2 for t=1,2,...,n
3     Sample x*~Q(·|x^(t-1))
4     Sample u~U[0,1]
5     if u<a(x^(t-1),x*) then
6         x^(t) = x*
7     else
8         x(t) = x(t-1)
9 return x^(0),x^(1),...,x^(n)

MCMC方法的问题

采集的相邻样本可能含有相关性
一种解决办法是每隔n个样本返回一个样本，但这样不仅提高了采样时的计算代价，而且也没有完全解除掉相关性；另一种解决办法是同时并行多个马氏链，每隔马氏链只采集一个样本。这两种方法是两个极端，更多从业者选择马氏链和小批量中的样本数接近，然后再用这些马氏链抽取所需要的样本。
无法预测混合期
检测一个马尔科夫链是否达到了平衡很困难，只能粗略的估计这段时间是否是足够的。由此可见，此算法的关键就是能否针对具体的π\piπ找到一个能够快速收敛的马氏链，收敛速度的快慢决定了算法的优劣。
由不变分布方程πP=π\pi P = \piπP=π以及不变分布的唯一性知，PPP的特征值中有且只有一个111，且π\piπ是对应的唯一特征向量。注意到PPP是不可约的，则由Perron-Frobenius定理可以保证这个矩阵最大特征值为1.这就导致了其他特征值都随时间变化趋向于0，那么第二大特征值将决定了收敛速度，也就是混合期的大小。
注：在M-H算法中，预设分布的标准差也影响着混合时间。如果标准差过小，则采集的样本多集中在概率大的周围；若标准差过大，则被拒绝的概率将变大，混合时间将变长。

吉布斯采样（Gibbs）

在计算过程中常常会遇到复杂的联合分布，这增加了计算的复杂度。事实上不同分量之间可能存在相关性，利用概率图模型，将联合分布分解成多个条件分布，可以减少状态空间的参数。在采样过程中，将高维总体的采样，化为一系列一维分布的取样，这就是Gibbs采样法的要义。

# 算法3：Gibbs算法
Proceduce MCMC_Sample(E, P_0, P, T)
1 Sample x^(0)~P_0
2 for t=1,2,...,T
3     Sample x_1^(t) ~ P(x_1|x_2^(t-1), x_3^(t-1),...,x_n^(t-1))
4     Sample x_2^(t) ~ P(x_2|x_1^(t), x_3^(t-1),...,x_n^(t-1))
5     ...
6     Sample x_j^(t) ~ P(x_j|x_1^(t),..., x_(j-1)^(t),x_(j+1)^(t-1),...,x_n^(t-1))
7     ...
8     Sample x_n^(t) ~ P(x_n|x_1^(t), x_2^(t),...,x_(n-1)^(t))
9 return x^(0),x^(1),...,x^(n)

参考文献

古德费洛, 本吉奥, 库维尔, 赵申剑, 黎彧君, and 符天凡. 深度学习. 北京: 人民邮电出版社, 2017.
钱敏平. 应用随机过程. 北京: 高等教育出版社, 2011.
罗斯, and 龚光鲁. 应用随机过程概率模型导论. 北京: 人民邮电出版社, 2011.
龚光鲁. 应用随机过程教程及在算法和智能计算中的随机模型. 北京: 清华大学出版社, 2004.
Andrieu, Christophe, De Freitas, Nando, Doucet, Arnaud, and Jordan, Michael I. “An Introduction to MCMC for Machine Learning.” Machine Learning 50.1 (2003): 5-43.
科勒, 弗里德曼, 王飞跃, and 韩素青. 概率图模型原理与技术. 北京: 清华大学出版社, 2015.
Perron, Oskar. “Zur Theorie Der Matrices.” Mathematische Annalen 64.2 (1907): 248-63.

马尔科夫蒙特卡罗方法相关推荐

马尔可夫蒙特卡罗 MCMC 原理及经典实现
我们在做机器学习.深度学习或自然语言处理等项目时,经常采用什么方法采样呢?大家马上会想到吉布斯 Gibbs 采样,今天我们来分享一种比较实用的采样方法:马尔可夫蒙特卡罗方法,吉布斯采样是其中的一种. ...
用Python中的马尔科夫链进行营销渠道归因
用Python中的马尔科夫链进行营销渠道归因 --第一部分:"更简单 "的方法任何积极开展营销活动的企业都应该对确定哪些营销渠道推动了实际转化率感兴趣.投资回报率(ROI)是一个 ...
通过Python实现马尔科夫链蒙特卡罗方法的入门级应用
通过把马尔科夫链蒙特卡罗(MCMC)应用于一个具体问题,本文介绍了 Python 中 MCMC 的入门级应用. GitHub 地址:https://github.com/WillKoehrsen/ai ...
MCMC+马尔科夫链蒙特卡罗
MCMC+马尔科夫链蒙特卡罗为了解决什么问题,所以出现了这一种方法? 后来又因为出现了什么新情况,所以产生了对应的变种?
MCMC原理解析(马尔科夫链蒙特卡洛方法)
马尔科夫链蒙特卡洛方法(Markov Chain Monte Carlo),简称MCMC,MCMC算法的核心思想是我们已知一个概率密度函数,需要从这个概率分布中采样,来分析这个分布的一些统计特性,然而 ...
蒙特卡罗马尔科夫链与Gibbs采样
这几个概念看了挺多遍都还是含混不清,最近看了几篇博客,才算大致理解了一点点皮毛,所以来总结一下. MCMC概述从名字我们可以看出,MCMC由两个MC组成,即蒙特卡罗方法(Monte Carlo Si ...
第十五课.马尔科夫链蒙特卡洛方法
目录 M-H采样 Metropolis-Hastings采样原理 M-H采样步骤 Gibbs方法 Gibbs核心流程 Gibbs采样的合理性证明 Gibbs采样实验在第十四课中讲述了马尔科夫链与其 ...
深度强化学习-马尔科夫决策过程和表格型方法
深度强化学习-马尔科夫决策过程和表格型方法-笔记(二)_wield_jjz的博客-CSDN博客深度强化学习2--马尔科夫决策过程(MDP)_谢宜廷的博客-CSDN博客 (零基础可以看懂)强化学习中的 ...
马尔科夫链和马尔科夫链蒙特卡洛方法
前言译自:<Training Restricted Boltzmann Machines: An Introduction > 马尔科夫链在RBM的训练中占据重要地位,因为它提供了从复杂 ...
基于张量的多元多阶马尔科夫多模态预测方法
本博客整理自研读的论文,文末会附上出处. 基于张量的多元多阶马尔科夫多模态预测方法一.问题背景二.多元多阶马尔科夫模型 1.张量连接和张量统一乘 2.多元多阶马尔科夫转移模型 3.多元多阶马尔 ...

马尔科夫蒙特卡罗方法

文章目录