标准采样方法

先从最基本的方法说起，可能很多人都知道只要可以对分布函数F(x)\displaystyle F( x)F(x)求逆，并从均匀分布中采样u,并将u代进逆函数中就能得到x的样本，即x=F−1(u),u∼U(0,1)\displaystyle x=F^{-1}( u) ,u\sim U( 0,1)x=F−1(u),u∼U(0,1)。他的原理是什么？其实他的出发点是找到一个从均匀分布到目标分布的可逆变换g\displaystyle gg：

x=g(u)p(x)=p(u)∣dudx∣=pu(g−1(x))∣dg−1(x)dx∣x=g( u)\\ p( x) =p( u)\left| \frac{du}{dx}\right| =p_{u}\left( g^{-1}( x)\right)\left| \frac{dg^{-1}( x)}{dx}\right| x=g(u)p(x)=p(u)∣∣∣∣dxdu∣∣∣∣=pu(g−1(x))∣∣∣∣dxdg−1(x)∣∣∣∣

因为我们要保证x是我们的目标分布，那么他一定要满足∫−∞xp(x)dx=F(x)\displaystyle \int ^{x}_{-\infty } p( x) dx=F( x)∫−∞xp(x)dx=F(x)，因为pu(f−1(x))\displaystyle p_{u}\left( f^{-1}( x)\right)pu(f−1(x))是0到1的均匀分布所以他等于1，于是必须有∫−∞x∣dg−1(x)dx∣dx=F(x)⟹g−1(x)=F(x)⟹g(x)=F−1(x)\displaystyle \int ^{x}_{-\infty }\left| \frac{dg^{-1}( x)}{dx}\right| dx=F( x) \Longrightarrow g^{-1}( x) =F( x) \Longrightarrow g( x) =F^{-1}( x)∫−∞x∣∣∣∣dxdg−1(x)∣∣∣∣dx=F(x)⟹g−1(x)=F(x)⟹g(x)=F−1(x)。
然而该方法的问题在于，并不是每个概率分布的都可以写出来并且可以求逆的。对于那些特别复杂的分布这样方法就无能为力了。

拒绝采样

现在我们的目标是要采样p(x)\displaystyle p( x)p(x)的分布。为了采样这个分布，实际上，我们只需要知道p~(x)\displaystyle \tilde{p}( x)p~(x)就足够了，p~(x)\displaystyle \tilde{p}( x)p~(x)是没有标准化的p(x)=p~(x)/Z\displaystyle p( x) =\tilde{p}( x) /Zp(x)=p~(x)/Z, 其中Z=∫p~(x)dx\displaystyle Z=\int \tilde{p}( x) dxZ=∫p~(x)dx 是不需要知道的。那么拒绝采样的流程是这样的，首先随便找一个proposal distributionq(x)\displaystyle q( x)q(x) 使得Mq(x)⩾p~(x)\displaystyle Mq( x) \geqslant \tilde{p}( x)Mq(x)⩾p~(x)，然后从q中随机一个样本x∼q(x)\displaystyle x\sim q( x)x∼q(x)，再从均匀分布中随机一个样本u∼U(0,1)\displaystyle u\sim U( 0,1)u∼U(0,1)，如果u⩾p~(x)Mq(x)\displaystyle u\geqslant \frac{\tilde{p}( x)}{Mq( x)}u⩾Mq(x)p~(x)则拒绝，否则接受。为什么这样就能从p中采样？

订正：图中的p\displaystyle pp应该是p~\displaystyle \tilde{p}p~

首先，因为u是0,1区间，所以一定有0⩽uMq(x)⩽Mq(x)\displaystyle 0\leqslant uMq( x) \leqslant Mq( x)0⩽uMq(x)⩽Mq(x)，那么u的作用其实就是选择0到Mq(x)\displaystyle Mq( x)Mq(x)之间的值，如果p~(x)⩽uMq(x)\displaystyle \tilde{p}( x) \leqslant uMq( x)p~(x)⩽uMq(x)，那么就拒绝掉样本，也就是上图白色部分，否则p~(x)⩾uMq(x)\displaystyle \tilde{p}( x) \geqslant uMq( x)p~(x)⩾uMq(x)就是落在黑色区域，则接受。显然，对于x轴上的任意一个点x′\displaystyle x'x′，其对应接受的概率就是p~(x′)Mq(x′)\displaystyle \frac{\tilde{p}( x')}{Mq( x')}Mq(x′)p~(x′)。不过每个点，因为q，抽到的概率都不一样，所以我们要把q抽中的概率也算上去。

那么用这个方法我们可以得到一个序列的样本(x1,accept1,...,xi,accepti)\displaystyle ( x_{1} ,accept_{1} ,...,x_{i} ,accept_{i})(x1,accept1,...,xi,accepti)，其中xi∼q(x)\displaystyle x_{i} \sim q( x)xi∼q(x)是由q\displaystyle qq产生的，而且每个pair(xi,accepti)\displaystyle ( x_{i} ,accept_{i})(xi,accepti)都是相互独立的，那么我们可以证明，p(x∣accept)\displaystyle p( x|accept)p(x∣accept)就等于我们想要的分布p(x)\displaystyle p( x)p(x):

p(x∣accept)=p(accept∣x)q(x)p(accpet)=p~(x)Mq(x)q(x)∫p~(x)Mq(x)q(x)dx=p~(x)∫p~(x)dx=p(x)/Z∫p(x)/Zdx=p(x)\begin{aligned} p( x|accept) & =\frac{p( accept|x) q( x)}{p( accpet)}\\ & =\frac{\frac{\tilde{p}( x)}{Mq( x)} q( x)}{\int \frac{\tilde{p}( x)}{Mq( x)} q( x) dx}\\ & =\frac{\tilde{p}( x)}{\int \tilde{p}( x) dx}\\ & =\frac{p( x) /Z}{\int p( x) /Zdx}\\ & =p( x) \end{aligned} p(x∣accept)=p(accpet)p(accept∣x)q(x)=∫Mq(x)p~(x)q(x)dxMq(x)p~(x)q(x)=∫p~(x)dxp~(x)=∫p(x)/Zdxp(x)/Z=p(x)

通过拒绝采样，我们可以知道，对于任意的分布，我们可以设计一种接受-拒绝的概率，从而使得所有接受的点都是该分布的样本。然而他的问题在于，他对propose distribution q的要求很高，试想下，如果q与真实p的差距过大，那么几乎每个样本都会被拒绝掉，效率很低。那如何设计更好的q\displaystyle qq呢？一个办法是设计一个t(x′∣x)\displaystyle t( x'|x)t(x′∣x)，他从上一次接受的样本作为条件，然后经过转换propose一个新的样本，如果我们能够约束这个转换不会改变目标分布的，那么我们就能快速的找到p(x)\displaystyle p( x)p(x)的样本。这也正是MCMC的思想。

IMCMC

然而MCMC的方法少说也有成百上千种，如何统一起来呢？我们从另一个角度来考虑MCMC。首先MCMC最基本可以从离散马尔科夫链说起，A是一个转移矩阵，π\displaystyle \piπ是一个向量，如果满足下述公式

πA=π\pi A=\pi πA=π

则称A\displaystyle AA的平稳分布是π\displaystyle \piπ。类似的在连续的情况下：

∫t(x′∣x)p(x)dx=p(x′)\int t( x'|x) p( x) dx=p( x') ∫t(x′∣x)p(x)dx=p(x′)

我们希望找到一个转换概率分布t(x′∣x)\displaystyle t( x'|x)t(x′∣x)，且该转换不会改变来自p(x)\displaystyle p( x)p(x)的样本的分布。事实上这个t不一定是随机的，当他是确定的映射时，我们有t(x′∣x)=δ(x′−f(x))\displaystyle t( x'|x) =\delta ( x'-f( x))t(x′∣x)=δ(x′−f(x)),相当于x′=f(x)\displaystyle x'=f( x)x′=f(x)于是
∫δ(x′−f(x))p(x)dx=p(x′)\int \delta ( x'-f( x)) p( x) dx=p( x') ∫δ(x′−f(x))p(x)dx=p(x′)

这意味着我们要找到一个映射函数，使得x′=f(x)\displaystyle x'=f( x)x′=f(x)并且他们的分布一样，又因为根据分布变换公式

px(x)=px′(x)px(x)=px′(f(x))∣df(x)dx∣=px(f(x))∣df(x)dx∣=px′(x)=px(f−1(x))∣df−1(x)dx∣p_{x}( x) =p_{x'}( x)\\ p_{x}( x) =p_{x'}( f( x))\left| \frac{df( x)}{dx}\right| =p_{x}( f( x))\left| \frac{df( x)}{dx}\right| \\ =p_{x'}( x) =p_{x}\left( f^{-1}( x)\right)\left| \frac{df^{-1}( x)}{dx}\right| px(x)=px′(x)px(x)=px′(f(x))∣∣∣∣dxdf(x)∣∣∣∣=px(f(x))∣∣∣∣dxdf(x)∣∣∣∣=px′(x)=px(f−1(x))∣∣∣∣dxdf−1(x)∣∣∣∣
（其实从上述公式中我们就能大致的能看出，f(x)=f−1(x)f(x)=f^{-1}(x)f(x)=f−1(x)这个条件的成立，这就是所谓的Involutive function，IMCMC就是围绕如何设计这个f函数来做的）
不过这样的f不好找，所以一般我们可以用以下这个：

t(x′∣x)=δ(x′−f(x))min⁡{1,p(f(x))p(x)∣∂f∂x∣}⏟Paccept++δ(x′−x)(1−min⁡{1,p(f(x))p(x)∣∂f∂x∣})⏟Preject,\begin{aligned} t(x'|x)=\delta (x'-f(x))\underbrace{\min\biggl\{1,\frac{p(f(x))}{p(x)}\biggl|\frac{\partial f}{\partial x}\biggl|\biggr\}}_{P_{\text{accept}}} +\\ +\delta (x'-x)\underbrace{\biggl( 1-\min\biggl\{1,\frac{p(f(x))}{p(x)}\biggl|\frac{\partial f}{\partial x}\biggl|\biggr\}\biggr)}_{P_{\text{reject}}} , \end{aligned} t(x′∣x)=δ(x′−f(x))Pacceptmin{1,p(x)p(f(x))∣∣∣∣∂x∂f∣∣∣∣}++δ(x′−x)Preject(1−min{1,p(x)p(f(x))∣∣∣∣∂x∂f∣∣∣∣}),

我们一般可以通过拒绝采样的方式来近似这个转换概率(这部分我不是很理解他跟拒绝采样的关系，没看懂)。这意味着，当x′=f(x)\displaystyle x'=f( x)x′=f(x)时，t(x′∣x)=min⁡{1,p(f(x))p(x)∣∂f∂x∣}\displaystyle t(x'|x)=\min\biggl\{1,\frac{p(f(x))}{p(x)}\biggl|\frac{\partial f}{\partial x}\biggl|\biggr\}t(x′∣x)=min{1,p(x)p(f(x))∣∣∣∣∂x∂f∣∣∣∣}，而当x′=x\displaystyle x'=xx′=x时，t(x′∣x)=(1−min⁡{1,p(f(x))p(x)∣∂f∂x∣})\displaystyle t(x'|x)=\biggl( 1-\min\biggl\{1,\frac{p(f(x))}{p(x)}\biggl|\frac{\partial f}{\partial x}\biggl|\biggr\}\biggr)t(x′∣x)=(1−min{1,p(x)p(f(x))∣∣∣∣∂x∂f∣∣∣∣})，显然，如果x∼p(x)\displaystyle x\sim p( x)x∼p(x)，且x′=f(x)\displaystyle x'=f( x)x′=f(x)，我们写的稍微简单一点，那么px′(x′)=t(x′∣x)p(x)=p(f(x))∣∂f∂x∣=px(f−1(x))∣df−1(x)dx∣\displaystyle p_{x'}( x') =t( x'|x) p( x) =p(f(x))\biggl|\frac{\partial f}{\partial x}\biggl| =p_{x}\left( f^{-1}( x)\right)\left| \frac{df^{-1}( x)}{dx}\right|px′(x′)=t(x′∣x)p(x)=p(f(x))∣∣∣∣∂x∂f∣∣∣∣=px(f−1(x))∣∣∣∣dxdf−1(x)∣∣∣∣。事实上，这等价于认为，f\displaystyle ff满足f(x)=f−1(x)\displaystyle f( x) =f^{-1}( x)f(x)=f−1(x)，这样的函数称为involutive function。然而，这个分布只是在两个点之间反复横跳，我们可以引入一个辅助变量v\displaystyle vv，我们随便抽取v\displaystyle vv的样本，从而得到不同的x′\displaystyle x'x′，这样就不会反复横跳了。可以证明，引入一个辅助变量后，只要满足f(x,v)=f−1(x,u)\displaystyle f( x,v) =f^{-1}( x,u)f(x,v)=f−1(x,u)，那么x的平稳分布仍然是目标分布。

t(x′,v′∣x,v)p(x,v)=t(x,v∣x′,v′)p(x′,v′).t(x',v'|x,v)p(x,v)=t(x,v|x',v')p(x',v'). t(x′,v′∣x,v)p(x,v)=t(x,v∣x′,v′)p(x′,v′).

可以证明，t的边缘分布仍然服从detailed balance
t^(x′∣x)=∫dvdv′t(x′,v′∣x,v)p(v∣x)t^(x′∣x)p(x)=t^(x∣x′)p(x′).\begin{aligned} \hat{t} (x'|x)=\int dvdv't(x',v'|x,v)p(v|x) \end{aligned} \begin{aligned} \hat{t} (x'|x)p(x)=\hat{t} (x|x')p(x'). \end{aligned} t^(x′∣x)=∫dvdv′t(x′,v′∣x,v)p(v∣x)t^(x′∣x)p(x)=t^(x∣x′)p(x′).

通过引入辅助变量，我们可以设计出iMCMC的算法，只需保证f是involutive function。

MH

一个特例是f(x,v)=v,x\displaystyle f( x,v) =v,xf(x,v)=v,x，他将x,v\displaystyle x,vx,v的位置交换了一下，这时候iMCMC等价于MH算法。显然他的接收概率为

P=min⁡{1,p(f(x,v))p(x)q(v∣x)}=min⁡{1,p(v)q(x∣v)p(x)q(v∣x)}P=\operatorname{min} \{1,\frac{p(f(x,v))}{p(x)q(v|x)} \}=\operatorname{min} \{1,\frac{p(v)q(x|v)}{p(x)q(v|x)} \} P=min{1,p(x)q(v∣x)p(f(x,v))}=min{1,p(x)q(v∣x)p(v)q(x∣v)}

HMC

另外一个经典的例子是HMC算法，他通过构造一个复杂的函数f\displaystyle ff，从而使得接受率非常高，几乎不会拒绝（实际做的时候就直接接受，不判断拒不拒绝），但是q\displaystyle qq又很简单，实际上启发了很多基于神经网络的MCMC就是在搞这个f\displaystyle ff。

HMC依赖于用leap-frog算子L\displaystyle LL来求解Hamiltonian dynamics的数值积分。对于L:[x(t),v(t)]→[x(t+ε),v(t+ε)]L:[x(t),v(t)]\rightarrow [x(t+\varepsilon ),v(t+\varepsilon )]L:[x(t),v(t)]→[x(t+ε),v(t+ε)]，

v(t+ε/2)=v(t)−ε2∇x(−log⁡p(x(t)))x(t+ε)=x(t)+ε∇v(−log⁡p(v(t+ε/2)))v(t+ε)=v(t+ε/2)−ε2∇x(−log⁡p(x(t+ε)))\left. \begin{array}{ l } v(t+\varepsilon /2)=v(t)-\frac{\varepsilon }{2} \nabla _{x} (-\operatorname{log} p(x(t)))\\ x(t+\varepsilon )=x(t)+\varepsilon \nabla _{v} (-\operatorname{log} p(v(t+\varepsilon /2)))\\ v(t+\varepsilon )=v(t+\varepsilon /2)-\frac{\varepsilon }{2} \nabla _{x} (-\operatorname{log} p(x(t+\varepsilon ))) \end{array}\right. v(t+ε/2)=v(t)−2ε∇x(−logp(x(t)))x(t+ε)=x(t)+ε∇v(−logp(v(t+ε/2)))v(t+ε)=v(t+ε/2)−2ε∇x(−logp(x(t+ε)))

后F\displaystyle FF则是一个简单的将辅助变量取负号的函数F:[x,v]=[x,−v]\displaystyle F:[ x,v] =[ x,-v]F:[x,v]=[x,−v]，可以论文中给出了证明FLk=(FLk)−1\displaystyle FL^{k} =\left( FL^{k}\right)^{-1}FLk=(FLk)−1，而且他的jacobian矩阵的行列式也是为1。里面还有很多算法，这里就不写了，有兴趣自己去看。

参考资料

Neklyudov, Kirill, et al. “Involutive MCMC: a Unifying Framework.” arXiv preprint arXiv:2006.16653 (2020).

How does the proof of Rejection Sampling make sense?

(ML 17.13) Proof of rejection sampling (part 1)

Murphy, Kevin P. Machine learning: a probabilistic perspective. MIT press, 2012.

Bishop, Christopher M. Pattern recognition and machine learning. springer, 2006.

MCMC算法大统一: Involutive MCMC相关推荐

应用随机过程张波商豪_Markov链的应用一：MCMC算法
本文是张迪同学对马尔链的应用的介绍应用一:Markov链在MCMC算法中的应用 1. MCMC概念 MCMC即马尔科夫链蒙特卡洛方法(Markov Chain Monte Carlo).该方法将Ma ...
机器学习算法——马尔可夫链蒙特卡罗（MCMC）
什么是MCMC,什么时候使用它 MCMC只是一种从分布中抽样的算法.这个术语代表"马尔可夫链蒙特卡罗",因为它是一种使用了"马尔可夫链"的"蒙特卡罗& ...
R语言--MCMC算法介绍以及例子
Markov Chain Monte Carlo (MCMC) 是一种用于随机生成满足给定分布的样本的算法.它通过构建一个马尔可夫链来模拟一个随机过程,从而生成样本. 首先,你需要选择一个初始状态,然 ...
人工智能是否存在「大统一理论」？
智源导读:在物理学当中,物理学家们经过近百年的努力,找到了一种理论模型,能够合理解释电磁相互作用.强相互作用和弱相互作用所导致的所有物理现象,这被称为"大统一理论"(Grand U ...
Android王者荣耀模拟金牌,王者荣耀：金牌算法大改，掌握这几点，让你金牌拿到手软...
原标题:王者荣耀:金牌算法大改,掌握这几点,让你金牌拿到手软在王者荣耀中,上分难的问题困扰一直都困扰着很多玩家,即使那些个人能力足够强大,面对最近版本的英雄频繁调整.兵线速度的大幅度提升,也会在没有 ...
【CV】10种轻量级人脸检测算法大PK ｜代码集合开源
点击上方"小白学视觉",选择加"星标"或"置顶" 重磅干货,第一时间送达最近在微信公众号 AIZOO 里看到轻量级人脸检测算法大盘点的文章 ...
爱因斯坦梦断“大统一理论”
来源:数学职业家爱因斯坦发表了他最为得意之作:广义相对论之后,便开始了他的"统一之梦".大有"躲进小楼成一统,管他冬夏与春秋"之势,这一"统&quo ...
GitHub的AI程序员“抄袭”算法大神代码，连原版注释都抄上了
晓查发自凹非寺量子位报道 | 公众号 QbitAI 本周GitHub官方和OpenAI联合发布了一款代码神器AI--GitHub Copilot,只需输入注释,即可自动生成代码,堪称一位&q ...
10种轻量级人脸检测算法大PK
几个月前,AIZOO曾经盘点过最强六大开源轻量级人脸检测项目分析 | 附打包下载,nihate同学将它丰富到10种算法,并用Python.对他们进行了汇总整理,以及效果的对比. Github链接:h ...
手机、平板、PC与智能电视实现数据大统一
进来,由 Ubuntu 手机原型设计引发的风波(数据大贯通),越演越烈,给人的感觉是"大雨欲来,风满楼".这是什么事情呢? 根据4月12日透露出的一份Ubuntu手机的功能设计示意 ...

MCMC算法大统一: Involutive MCMC

标准采样方法

拒绝采样

IMCMC

MH

HMC

参考资料

MCMC算法大统一: Involutive MCMC相关推荐

最新文章

热门文章