任上;唐亚勇
【摘 要】应用M-H综合框架统一描述哈密尔顿蒙特卡洛算法(HMC)、黎曼流形哈密尔顿蒙特卡洛算法(RMHMC)以及拉格朗日蒙特卡洛算法(LMC),分别估计三个算法的最优接受概率并对比.通过模拟有明显几何性质的banana-shaped分布,对比三个算法的接受概率,运用有效样本规模指标(ESS)和无效因子(IF)同时度量各算法生成样本的有效性.分析了三个算法的优势和不足,提出了一些改进设想. 【期刊名称】《滨州学院学报》
【年(卷),期】2017(033)002
沃岭生【总页数】7页(P59-65)
【关键词】HMC;RMHMC;LMC;M-H综合框架;banana-shaped分布
【作 者】任上;唐亚勇
【作者单位】马尔斯四川大学 数学学院,四川 成都 610065;四川大学 数学学院,四川 成都 610065眼动仪
【正文语种】中 文
【中图分类】泰诺星球O211
马尔科夫链蒙特卡罗方法(MCMC)有共同的思想，实现的方式则多种多样。经典MCMC算法在构造马尔科夫链时就有不同的形式，如Metropolis-Hastings算法使用舍选法构造马尔科夫链，Gibbs算法通过条件分布交替更新，尽管同属于马尔科夫链蒙特卡罗方法，却很难用一个综合性的框架来描述这些MCMC算法。Tran&Kohn[1]给出了通用的MCMC算法框架，主要基于M-H算法。
基础的MCMC算法在模拟有明显几何性质的分布时通常效率较低，收敛速度慢，需要的迭代次数多，基于哈密尔顿动力系统的一类马尔科夫链蒙特卡罗方法则能在较少的迭代次数中高效地生成样本。Neal[2]总结了结合分子动力系统的MCMC算法，命名为哈密尔顿蒙特卡洛算法(HMC)，发现在高维分布的模拟中，HMC的抽样效率相比随机游走Metropolis算法显著提高。由于参数空间具有黎曼几何性质，欧几里得空间上的动力系统不适合用来指
导变量更新，Girolami&Calderhead[3]提出了RMHMC算法，使得算法可以局部自适应，提高了HMC算法的效率。在RMHMC算法的基础上，Lan[4]把G(θ)-1p整体替换成速度v。从形式上让表达式变得简洁，同时减少了计算量。
本文在HMC、RMHMC和LMC这三个算法中，分别引入M-H综合框架，重新阐述了三个算法，使得这三种基于哈密尔顿动力系统下的MCMC算法有更加统一的理论描述。从而有利于进行算法对比，并有助于进一步研究这类算法和其他种类MCMC算法的结合。
1.1 基本思想
MCMC算法中的任何一个计算步骤，必然是确定性计算和随机数生成两者之一，而在确定性计算时，可以用映射来描述。因此在传统M-H算法的基础上引入可逆映射，可以解释很大部分的MCMC算法。而根据定理1[1]，可知在传统M-H算法中加入自逆映射构成的M-H综合框架依旧满足不变性条件：
其中Κ(dθk+1|θk)表示转移核，π(θ)为目标密度。
定理1[1] 当映射T(·)是连续可微的自逆映射
时，有θk～π(·)⟹θk+1～π(·)。
定义1 M-H综合框架[1]
(1)生成服从建议分布的随机数υk～q(υk|θk);
(2)把随机数υk映射到建议值ξk,并计算广义接受概率α(θk,υk)，
其中|ϑΓ(θk,υk)|是Γ取值在(θk,υk)上的雅克比行列式；
(3)以广义接受概率判断是否接受更新值，即先生成uk～μ(0,1),然后令
1.2 基于M-H综合框架的HMC算法
HMC算法的核心思想是在MCMC算法中引入哈密尔顿动力系统，把生成建议值的过程看作服从哈密尔顿动力系统的分子运动，将目标分布对应的密度函数f(θ)转化为势能方程U(q)=-log[f(θ)],变量θ作为位置变量q,最后再引入辅助变量p作为动量变量，对应的动能函数为K(p)为各成分的方差)。哈密尔顿动力系统本身是一个关于联合状态(q,p)的能量方程，总能量H(q,p)=U(q)+K(p)具有不变性，其轨迹在一个具有常数概率密度的超曲面中运动。因为满
足时间可逆性、哈密尔顿函数的不变性、保持状态空间体积不变性和辛性质[2]，哈密尔顿动力系统可以用来构造MCMC更新机制。实际应用中，HMC算法的更新分为两步，先随机生成辅助变量p，然后进行Metropolis更新，其中包括L步广义leapfrog迭代。根据M-H综合框架，假设在t时刻，θ(t)=θ,p(t)=p,HMC算法生成t+1时刻的(θ,p)可以描述为以下三步：
(1)生成随机向量p～N(0,I);
(2)依据映射T：[θ,p,r][ξ,w,r′]生成建议值；
(3)计算接受概率
其中|JT(θ,p)|是映射T的雅可比行列式在点(θ,p)处的取值。
米尼兹由于每一步Metropolis更新包含L步广义leapfrog迭代，设映射τ表示一步广义leapfrog迭代，
则映射T=τL,|JT(θ,p)|=|Jτ(θ,p)|L。由于广义leapfrog迭代是对描述哈密尔顿动力系统的微分方程的离散化，由哈密尔顿动力系统的时间可逆性[2]可知，映射τ是自逆的，即τ(τ(θ,υ))=(θ,υ)。映射T是τ的L重复合，
因此，证明了映射T也是自逆的，从而|JT(θ,p)|=1,接受概率α(θ,p)=min{1,exp[H(θ,p)-H(ξ,w)]}。根据哈密尔顿函数的不变性[2]，=0,知道当ε足够小的时候，H(θ,p)-H(ξ,w)无限接近0，因此，接受概率无限接近1。
1.3 基于M-H综合框架的RMHMC算法
统计模型的参数空间是一个黎曼流形，而密度函数空间的几何结构也由这个黎曼流形和相关的度量张量决定[3]。使用MCMC算法时，若已知参数空间对应的几何结构，再适当的引入一个位置特异矩阵G(θ)，自然能够提高转移效率，生成的样本既服从流形的几何性质，又有很好的遍历性。Girolami&Calderhead[3]在HMC算法的基础上，引入向量(速度乘以质量)，对应的动量K(p)等于在矩阵M的度量下的平方范数
矩阵M可以用参数空间对应的黎漫流形的度量张量G(θ)来定义，则的平方范数为‖通常取G(θ)=-Ey|θ[log{p(y,θ)}],即经验Fisher信息矩阵。因为G(θ)是基于观测数据y得到，而不是通过渐进抽样机制得到，更能保证抽样过程中样本不会因为自适应而产生过大的相关性。向量p满足p|θ～N(0,G(θ)),对应的哈密尔顿的函数为
对应的正规则分布的密度函数为exp{-H(θ,p)}=p(θ,p)=f(θ)p(p|θ),其边缘分布的密度函数与目标分布密度函数有同样的核，即
RMHMC算法参数更新的步骤分为生成随机数和Metropolis更新两步，每步Metropolis更新包含L步广义leapfrog迭代。根据定义1的M-H综合框架，假设在t时刻，θ(t)=θ,p(t)=p,RMHMC算法生成t+1时刻的(θ,p)主要分为以下三步：
(1)生成随机数p～N(0,G(θ));
(2)依据映射T：[θ,p,r][ξ,w,r′]生成建议值；
(3)计算接受概率：
其中|JT(θ,p)|是映射T的雅可比行列式在点(θ,p)处的取值。
由于每一步Metropolis更新包含L步广义leapfrog迭代，设映射τ表示一步广义leapfrog迭代，
p(t+1/2)为中间变量，则映射T=τL。由于广义leapfrog迭代是对描述哈密尔顿动力系统的微
分方程的离散化，由哈密尔顿动力系统的时间可逆性可知，映射τ是自逆的即τ(τ(θ,υ))=(θ,υ)。映射T是τ的L重复合，由T=τL知，映射T也是自逆的，从而|JT(θ,p)|=1,接受概率
根据哈密尔顿函数的不变性[2]，=0,知当ε足够小时，H(θ,p)-H(ξ,w)无限接近0，因此，接受概率无限接近1。但实际计算时，RMHMC的每步广义leapfrog迭代都有两个隐式方程(1)和(2)，需要不动点迭代求解，因此会造成一定误差，导致RMHMC接受概率比HMC小一些。
1.4 基于M-H综合框架的LMC算法
在黎曼流形下的哈密尔顿动力系统中，G(θ)-1p多次出现，于是Lan[4]把整体替换成速度v。从形式上让表达式变得简洁，同时满足动力系统的物理定义，动量p等于质量乘以速度。为了用速度v代替动量p，先改写哈密尔顿函数，则
令v=G(θ)-1p,(θ,v)对应的动力系统不再满足哈密尔顿动力系统，而是满足第二类欧拉-拉格朗日方程 。
西方唯美主义
假设在t时刻，θ(t)=θ,v(t)=v,在M-H综合框架下，LMC算法生成t+1时刻的(θ,v)主要分为以下三步：
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议