统计计算第六、七节课，Mante Calor方法（三）——MCMC

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

MCMC方法是将抽样与马氏链结合的抽样方法，这次的内容主要是MCMC方法下的两个子方法，第一个是Metropolis-Hastings（MH）算法，第二个是Gibb Sampling（GS），最后还会谈到一些MCMC的评估方法

背景：之前拒绝抽样方法和IS方法有两个缺点：
（1）q(x)很难找到合适的，而且M较难计算
（2）高维时抽样将变得困难

MH方法尝试解决缺点（1），GS方法尝试解决缺点（2）

MH方法

直观：q(x)很难找到合适的，虽然之前我们讲过自适应IS方法，但是q(x)只能限制在某个分布族，而且计算起来较为复杂，MH方法则是将MC和马氏过程结合，希望构造的马氏过程最终能够收敛到目标分布

马氏过程简介

我们在这里提一些接下来将用到的离散时间马氏过程的性质，由于我们不追求严格定义，所以下文在推导一些式子时过程也不会很严格。

离散空间版本（略）

连续空间版本

注：以下提到的集合均为正测度集合
转移密度函数：记从x转移到y的转移密度函数为$p(x,y)$
转移概率：记从x转移到集合A的转移概率为$P(x,A)$
不可约性：$\exists n>0,s.t.P_n(x,A)>0$
非周期性（略，严格定义有点麻烦，大概是说只要N足够大，对任意n>N和集合A，都有$P_n(x,A)>0$）
不变分布：$\pi (A)={\int }_{\Omega }P(x,A)\pi (dx)$，如果$\pi (x)$有密度函数$\dot{\pi }(x)$，则可写成$\dot{\pi }(y)={\int }_{\chi }p(x,y)\dot{\pi }(x)dx$
细致平衡条件：$\pi (dx)P(x,dy)=\pi (dy)P(y,dx)$，如果$\pi (x)$有密度函数$\dot{\pi }(x)$，则可写成$\dot{\pi }(x)p(x,y)=\dot{\pi }(y)p(y,x)$，由细致平衡条件可以推出不变分布为$\pi (x)$
遍历定理：若马氏链不可约，非周期，有不变分布$\pi$时，则不变分布唯一且马氏链是遍历的，即$$forallA\subseteq \Omega ，\underset{n\to \infty }{\lim }P(x,A)=\pi (A)$$从遍历定理中可以得到我们在最开始想要的结果，即马氏链最终能够收敛到目标分布，只要我们能构造一个马氏链满足不可约，非周期以及目标分布是不变分布$\pi (x)$
遍历定理加强版：条件同上，则对任意满足$E_{\pi }f(x)<\infty$的$f(x)$，都有$$\underset{T\to \infty }{\lim }\frac{1}{T}\sum _{t=1}^{T}f(X_t)=\int f(x)\pi (x)dx$$这个结果可以直接得到我们Mante Calor方法一开始想要计算的积分。

可以看出问题的关键是如何构造转移密度函数（以下称为proposal）能够让马氏链收敛到目标分布。

MH方法

直观：我们希望构造一个和最初的MC方法类似的算法，只不过这次的proposal是适应的，即随着当前所在的位置而变化

算法
step1：构造转移密度函数$Q(x,y)$，通常为以当前位置为中心的均匀分布或正态分布，主要是好算，设目标分布为$P(x)$
step2：从Q中抽样得到y
step3：计算接受概率 a ( x , y ) = min ⁡ { 1 , P ( y ) Q ( y , x ) P ( x ) Q ( x , y ) } a(x,y)=\min\{1,\frac{P(y)Q(y,x)}{P(x)Q(x,y)}\} a(x,y)=min{1,P(x)Q(x,y)P(y)Q(y,x)​}

性质
可以看出proposal与当前所在的位置有关，特别的如果proposal，MH方法就和拒绝抽样方法相似

收敛性
验证细致平衡条件$\pi (dx)P(x,dy)=\pi (dy)P(y,dx)$即可

优点
$Q(x,y)$需要满足的条件少且容易达到

缺点
想找到一个能让马氏链收敛速度较快的$Q(x,y)$很困难，下面我们将会谈到这种方法构造出的抽样过程的评估方法

GS方法

GS方法是一种特殊的MH方法，对MH方法在抽样高维变量上作出一些改进

算法
step1：选择$i$，从分布$p(x_i|x_1,...,x_{i-1},x_{i+1},...,x_d)$中抽样
step2：更新$x_i$
从算法步骤中我们需要考虑三件事：
（1）如何选择$i$：通常是在{1,…,d}中轮流一圈，或者从其中随机抽样，这里需要说的是，分量更新顺序的不同会导致总体抽样效果的不同
（2）条件分布需要容易从中抽样，不然GB方法就失去其处理高维变量的优势，一般来说，条件分布需要有显式表达式
（3）与MH方法的关系：上面我们提到GB方法是特殊的MH方法，原因在于在GB方法中，$Q(x,y)$是$Q((x_i,x_{-i}),(y_i,x_{-i}))=p(y_i|x_{-i})$，此时$a(x,y)$恒为1

如何更好地应用GS方法
主要还是如何设计抽样顺序，当某个分量上的抽样效果不好时（下面会介绍诊断中的评估准则），可以在该分量上多进行几次抽样并不断更新该分量，只取最后一次作为输出，即在该分量上走了更多步

MCMC的诊断

Mixing Rate（混合率）

首先画出MCMC的变量的轨迹图，如果是高维变量则画出分量的轨迹图（轨迹图是将分量第$i$次取值的折线图），如果轨迹图上下震动较为剧烈，则说明混合得较好，抽样效率较高

Autocorrelation（自相关）

滞后s步的自相关系数的计算公式为$$\rho (s)=\frac{cov(X_i,X_{i+s})}{Var(X_i)}$$其中协方差和方差都用抽样得到的样本估计，然后画出$\rho (s)$的图像，如果随着s的增大，$\rho (s)$很快降为0，就说明混合得较好

Effective Sample Size（有效样本）

计算公式为$$ESS=\frac{n}{1+2{\sum }_{i=1}^{\infty }\rho (i)}$$这个指标是衡量平均自相关性的，直观上还可以解释为相对于多少个相互独立样本

是否漏找mode（密度较高的点）

MCMC可能陷入某个密度高的点附近出不来，这是极其致命的，但上述指标检测不出这件事，目前也没有普适的指标可以用，一般的做法是同时运行多组MCMC，看是否漏找mode