统计计算第七、八节课，Mante Calor方法（四）——辅助变量法

这是我上的统计计算课讲的主要内容，写在这可以互相交流，有些地方我不是很理解的会标出来（用加粗斜体*标出），求大佬在留言处表达自己的看法，另外如果有啥问题也可以在留言处留言，如果我看到了会回复

这次的主要内容有Slice方法和HMC方法，其中Slice方法是辅助变量法和GS方法的结合，HMC方法可以算是Slice方法的改进，当然也会谈一些其他零碎的方法

解决的问题：MCMC方法可能会混合不好或遗漏mode，实际上遗漏mode也算是一种混合不好

辅助变量法简介

直观想法：我们可以想象一个物体在三维空间上的运动会比在二维空间上的运动更自由，可以选择的方向更多，类似的，当我们在原来的分布上进行抽样混合得不好时，我们可以增加变量，构造新的高维空间上的分布，只要新分布在原来变量上的边际分布和原分布相同，那么我们在新分布上进行抽样后取原来变量上的值就是想要的

下面我们先举个辅助变量法的例子说明这一点

平行退火法

算法
设目标分布$f(x)=e^{-U(x)}$
step1：取$1=T_0<T_1<\cdots <T_K$，构造密度函数$f_k(x)$满足$$f_k(x)\propto e^{-\frac{U(x)}{T_k}}$$目标是从联合密度函数$p(x_0,...,x_K)\propto f_0(x)\cdots f(x_K)$中抽样
step2：同时跑K+1个MCMC，设当前样本点为$x_{t-1}^0,...,x_{t-1}^K$，则更新后的样本点为$x_t^0,...,x_t^K$
step3：选择两个温度$i,j$，交换$x_t^i,x_t^j$，其中交换的接受概率为 min ⁡ { 1 , f i ( x t j ) f j ( x t i ) f i ( x t i ) f j ( x t j ) } \min\{1,\frac{f_i(x_t^j)f_j(x_t^i)}{f_i(x_t^i)f_j(x_t^j)}\} min{1,fi​(xti​)fj​(xtj​)fi​(xtj​)fj​(xti​)​}，返回step2

可以看出step3是此算法的核心步骤，通过计算我们可以知道当温度较高时，密度函数的图像趋于水平线，在不同mode之间的跳跃也就比较容易，因为mode之间的“山谷”较为宽敞，接受概率较高，所以把高温度的点和低温度的点交换就能让低温度的点在不同的mode之间跳跃，恰好目标分布是温度最低的点，这也就增强的目标分布的样本点的混合性

收敛性
根据之前的直观想法，只需验证$x_0,...,x_K$的联合分布在$x_0$上的边际分布是目标分布即可，首先可以验证step2是进行的MCMC过程的不变分布为$p(x_0,...,x_K)$，由各分量间的独立性得到，其次可以验证step3也保持不变分布为$p(x_0,...,x_K)$，可以直接验证细致平衡条件，因此最终的不变分布即为$p(x_0,...,x_K)$，而$p(x_0,...,x_K)$在$x_0$上的边际分布就是目标分布

如何能让平行退火算法表现得更好
调整超参数：温度的设计，如何设计合适的温度是最重要的，温度的数目太多影响收敛效率，温度太低难以在不同mode间跳跃，温度太高会忽视目标分布本身的mode特征，所以调整温度的超参数是一件很困难的事
交换的选择：这个和（三）中的GS方法一样，选择的方法有很多，可以设计成固定的，也可以随机选择，这件事比调超参数更没头绪

以上是对辅助变量法的一个简介，下面进入本次的重点Slice和HMC

Slice方法

直观想法
把从目标分布抽样看作从目标分布与X轴围成的区域上的均匀分布抽样，从这片区域上的均匀分布抽样又可以用GS方法

算法
设目标分布为$f(x)$，构造新分布$f(x,u)=f(x)f(u|x)$，其中$U$是服从均匀分布的随机变量，$f(u|x)$是$(0,f(x))$上的均匀分布
step1：给定当前状态x，从$U(0,,f(x))$中抽样，更新u
step2：给定当前状态u，令 S = { x : f ( x ) ≥ u } S=\{x:f(x)\geq u\} S={x:f(x)≥u}，从$U(S)$中抽样，更新x，返回step1

改进
可以看出step2是较为困难的，因为$S$的范围不好确定，并且从$U(S)$中抽样也比较困难，所以进行如下改进
step2’(a)step-out：找一个x的邻域使得此邻域两端的函数值小于x的函数值
（b）shrinkage：从上述邻域的均匀分布上抽样，满足新样本的函数值大于x的函数值，更新x

可以看出step-out步更加困难，因为如果邻域选小没达到要求的话，首先要一步一步地向外扩展，影响效率，其次容易跳不出x所处的mode，这就没达到最初的目的了，如果邻域选太大了，shrinkage步骤的拒绝概率就会变高，影响效率

HMC方法

目前抽样的主流方法几乎都是在这个方法上的改进

背景
之前的MCMC方法的proposal虽然与当前所在的点有关，但形式上都是均匀分布，正态分布之类的普适的东西，而且分布中的参数也是固定的，没有用到目标分布在当前点的形状的信息，即导数的信息，这样的抽样在高维时很没有效率，会出现（三）中MCMC诊断中的各种问题。
在介绍HMC方法之前，先介绍一个简单的方法来理解HMC的直观想法

Metropolis Adjusted Langevin Algorithm方法(MALA)

这是在MCMC上的一个简单的改进，利用了目标分布的一阶导数信息，可以类比于梯度下降法

改进的地方
设目标分布为$f(x)$
$$Q(x,y)=N(x+\delta ▽logf(x),{\sigma }^{2}I)$$其中$\delta$取太大和太小都不好， 实践中经常取$\frac{{\sigma }^2}{2}$

HMC方法

HMC的直观想法就是来自上面的MALA，不过HMC的类比更像是更加高效的矩方法，使用的是更加精细的微分方程模型

记号
设目标分布为$p(x)\propto e^{-U(x)}$，$U(x)=-log(P(x))$其中P(x)是没有规范化过的，成U(x)为势能函数；引入动量$r$，$r$携带动能$K(r)=\frac{1}{2}{r}^T{M}^{-1}r$，构造总能量函数$H(x,r)=U(x)+K(r)$，也称作哈密顿函数。
令$x,r$是时间$t$的函数$x(t),r(t)$，则x，r需满足哈密顿方程$$\frac{dx}{dt}=\frac{dH}{dr},\frac{dr}{dt}=-\frac{dH}{dx}$$此方程在数学上的直观解释为$r$的的方向为目标分布的导数方向（可以先把M当作单位矩阵），$x$前进的方向是$r$
在物理上的直观解释为x可以类比为力，r可以类比为速度（在物体质量为1时）
则$(x,r)$的联合分布为$$p(x,r)\propto e^{-H(x,r)}\propto p(x)N(r|0,M)$$x的边际分布就是p(x)

算法
step1：从$N(0,M)$中抽样，更新r
step2：模拟上述微分方程$$(x,r)=(x(0),r(0))\to (x(t),r(t)),(x^{\prime },r^{\prime })=(x(t),-r(t))$$

收敛性：（留个坑，没想明白）
性质：满足哈密顿方程的$x(t),r(t)$有很好的性质
（1）轨道可逆性（验证细致平衡的必要条件）
（2）$H(x(t),r(t))$恒为常数（只需证导数为零），这点保证了step2中接受概率为1
数值模拟微分方程：关于微分方程的模拟，我们不采用欧拉格式，而采用下面更为精细的蛙跳（leap-frog）格式$$\begin{gathered} r(t+\frac{ϵ}{2})=r(t)-\frac{ϵ}{2}\frac{\partial U}{\partial x}(x(t)) \\ x(t+ϵ)=x(t)+ϵ\frac{\partial K}{\partial r}(r(t+\frac{ϵ}{2})) \\ r(t+ϵ)=r(t+\frac{ϵ}{2})-\frac{ϵ}{2}\frac{\partial U}{\partial x}(x(t+ϵ)) \end{gathered}$$可以发现此数值格式保持性质（1），但不保持性质（2），所以应该在算法中增加step3：接受概率为 a ( x ′ , r ′ ∣ x , r ) = min ⁡ { 1 , e − H ( x ′ , r ′ ) + H ( x , r ) } a(x',r'|x,r)=\min\{1,e^{-H(x',r')+H(x,r)}\} a(x′,r′∣x,r)=min{1,e−H(x′,r′)+H(x,r)}

改进（RMHMC方法–基于黎曼流形的方法）

上述方法中M是固定的，但我们希望利用目标分布局部的几何性质来更新M，这种方法可以类比为牛顿法。最开始的HMC方法有个缺点，在模拟时会在等高线为狭窄的椭圆的短轴方向上不断震荡，这样会让效率变低，我们希望轨道能沿着狭窄的椭圆的长轴方向前进，并突破当前mode，提高效率，为此，我们希望利用目标分布的变量所在空间的二阶导信息。

M的选择
首先我们再明确一下问题，我们使用贝叶斯方法求参数的后验分布，并从改后验分布中抽样，设对数似然函数为$L(\theta |y)=log(f(y|\theta ))=\sum log(f(y_i|\theta ))$，对数后验分布为$log(p(\theta |y))$，则我们的目标分布就是后验分布，现在我们想利用的是$\theta$所在空间的二阶导信息，而每一个$\theta$代表一个似然函数，所以我们考虑的是函数空间，选择的距离为之前定义过的KL距离，但由于该距离不具备对称性，所以需要稍微修改一下，重新定义对称KL距离$$D_{KL}^S(p||q)=D_{KL}(p||q)+D_{KL}(q||p)$$通过对$\theta$泰勒展开并进行简单的估计可以算出$$D_{KL}^S(f(y|\theta +△\theta )||f(y|\theta ))\approx △{\theta }^T{E}_y({▽}_{\theta }L(\theta |y{)}^T{▽}_{\theta }L(\theta |y))△\theta =△{\theta }^{T}G(\theta )△\theta$$其中$G(\theta )$是Fisher信息量，则重新定义哈密顿函数为(x相对于$\theta$)$$H(x,r)=U(x)+\frac{1}{2}log((2\pi {)}^d|G(x)|)+\frac{1}{2}{r}^{T}G(x{)}^{-1}r$$此时密度函数为$$p(x,r)\propto p(x)N(r|0,G(x))$$可以看出r的分布与x有关了，因此模拟哈密顿方程的数值格式也需要改变$$\begin{gathered} r(t+\frac{ϵ}{2})=r(t)-\frac{ϵ}{2}{▽}_{x}H(x(t),r(t+\frac{ϵ}{2})) \\ x(t+ϵ)=x(t)+\frac{ϵ}{2}(G(x{)}^{-1}+G(x+ϵ{)}^{-1})r(t+\frac{ϵ}{2})) \\ r(t+ϵ)=r(t+\frac{ϵ}{2})-\frac{ϵ}{2}{▽}_{x}H(x(t+ϵ),r(t+\frac{ϵ}{2})) \end{gathered}$$上式中第一个式子是隐式格式，因此我们将使用欧拉格式得到$r(t+\frac{ϵ}{2})$，至此我们成功地选择了根据x变化而变化的M，而且也得到了（x，r）的数值更新格式，改进后的方法的x(t)的轨道将会平滑很多，效率也有所提高

超参数t的选择
如果t太小，有可能混合不均匀，如果t太大，根据性质2有可能转回初始点，这时候需要设置一个停止准则，直观上可以选择No-U-turn sampler，这里不细谈，总之要尽可能跑得远。