MCMC-Markov Chain Monte Carlo

MCMC-Markov Chain Monte Carlo

什么是采样…

如果说相对一个高斯分布进行采样
$$f(x)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}}$$
其含义是我要随机的去取一堆x的值，这些x的值符合$f(x)$的概率密度分布。但是应该如何采呢？首先采样要保证足够的随机性，也就是说所有的样本值我都有相同的概率取到，那么对x轴所有的数均匀采样，每个x的概率密度为$\frac{1}{N}$,显然不符合我们要的$f(x)$分布。

因此采样的等可能性，并不是针对样本值的等可能性，不是对x轴的等可能性，因为我们要使得$f(x)$大的地方多取到，$f(x)$小的的地方少取到。 难不成我每取一个x就计算一下$f(x)$，$f(x)$大的我就要，小的我就适当抛弃？那这就是个唯心的问题了，不是随机采样了。那怎么办？

再来观察一下正态分布的cdf，概率分布函数
当我们对y轴，也就是概率取等可能采样，$u\sim U(0,1)$,后用cdf的反函数求出对应$x=cd{f}^{-1}(u)$,可得很显然，这时候得到的x是我们想要的采样，那么为什么？

对y轴的均匀采样使得$\Delta cdf$相同，而$\Delta cdf=\sum pdf$,也就是说$cdf$的增加值等于这个小的区间内所有样本值$x$的$pdf$的累加，那意味着$\Delta cdf$相同的情况下，$pdf$越高的地方，需要累加的$x$的数量越少，对应$x$轴上的区间越窄，$x$越密集。此乃 inverse of CDF。

所以说，均匀采样，等概率的采样，是针对样本的概率分布函数均匀采，从而得到的样本值才会符合我们要的$pdf$。

而问题在于不是所有的概率分布都可以有逆函数，大多数情况不能这么简单的去做，于是有了许多牛逼的采样方法。

Rejection Sampling

绿色的曲线使我们要去采样的复杂的分布$p(x)$，红色的曲线使我们生成的辅助的概率分布$q(x)$，乘以M是为了保证红色曲线处处都在绿色曲线上部，其实红色曲线并不是一个正确的概率分布，对绿色曲线积分，一定是1，因为积分就是总的概率，为红色积分一定大于1，哪有概率大于1的，不合理，但也无所谓，反正是为了方便采样辅助用的。

进行采样：

1. 随机给一个$x$，随机生成一个$u\sim U(0,1)$,
2. 如果$u<\frac{p(x)}{Mq(x)}$,则这个样本值$x$我要了，否则扔了，
3. 再重新取x 重新取u

这里的均匀性体现在$u$上，因为$u$是一个0到1的均匀分布的取值，因此即便有时候$\frac{p(x)}{Mq(x)}$很小，也就是说绿色和红色差的很大，这个$x$也还是有可能取到，只不过取到的可能性小了很多。而反之绿色和红色曲线越接近的地方，这个局部取到的$x$我接受的概率就大，取到的x就多就密。要是两个曲线能完全重合，那采样的点完全符合绿色的分布，没一丁点毛病。

缺点是肯定有的么，比如红色曲线怎么选，选不好你采样100个，可能80个都要被扔掉，效率太低。

Adaptive Rejection Sampling

所以有了ARS。

目的是为了生成尽可能好的，接近目标概率分布辅助概率分布。

拿正态分布为例子，搞一下它的概率密度函数
$$log(f(x))=log(\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}})=-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}\cdot \frac{1}{\sigma \sqrt{2\pi }}$$
这是一个凸函数，下左图画反了，右图为正态分布的原概率密度函数。

在二次曲线上取很多点，求其切线，由许多切线段形成的折线将二次曲线包起来，由凹函数的性质可知这个折线段肯定在这个二次函数上方。

接下来将这个折线段取指数，相当于映射到了原概率密度函数的空间中，得到的就是一段一段的指数函数的曲线，因为原折线段在二次曲线上，那取了指数也都在原概率密度函数上，这样生成的曲线肯定在目标函数上方，这样就可以用这个辅助曲线去做Rejection Sampling了，而且当切线段取得很密很密的时候，得到的指数函数结合体就会很接近正态分布函数，就很完美。

但是这里是正态分布，别的分布要保证$logf(x)$是一个凸函数，convex，才能用这个法子。

Importance Sampling

直接上公式
$$\begin{array}{rl}{E}_{p(x)}[f(x)]& =\int f(x)p(x)dx\\ & =\int \underset{new\widetilde{f(x)}}{\underbrace{f(x)\frac{p(x)}{q(x)}}}q(x)dx\\ & =\frac{1}{N}\sum _{n=1}^{N}f(x^i)\frac{p(x^i)}{q(x^i)}\end{array}$$
目的是求函数$f(x)$在$p(x)$分布下的期望，而$p(x)$如果很复杂，这个积分不好算，这时候引入新的概率分布$q(x)$,这时候变成求$f(x)$在$p(x)$分布下每个样本带有一定权重的期望。

这个跟之前的Rejection Sampling很想，都要考虑原概率分布和新引进的概率分布之间的比值，但区别是Rejection Sampling考虑完比值要决定此数之去留，这里来了就留下只不过权重有高低。

MCMC-Markov Chain Monte Carlo

上文提到的每种采样，每次采得的样本之间没有关系，先取哪个后取哪个无所谓总之最后取出来一堆数据就好。但是MCMC与众不同，下一个数据点与上一个数据点有关联。

说到这种下一个状态和上一个状态有关的情况，必然想到卡尔科夫链，MC的条件就是
$$P(X_{t+1}=x_{t+1}|X_t,...X_0)=P(X_{t+1}=x_{t+1}|X_t)$$
下一个状态只和上一个状态有关系。
其状态转移矩阵（马尔科夫矩阵）为
$$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right]$$
其中$a_{ij}$代表上一个状态为$i$时下一个状态为$j$的概率，其性质为每一行的总和为1。这会导致转移矩阵得到特征值一个为1，其余皆为小于1的数。而
$$\begin{array}{r}P(Xn)=A^{n}P(X_0)\end{array}$$
其中已知矩阵A的特征向量线性无关，因此可以用特征向量来表示$\overrightarrow{P(x^0)}=a_1\overrightarrow{v_1}+a_2\overrightarrow{v_2}+...+a_n\overrightarrow{v_n}$
因此上式得
$$\begin{array}{rl}\overrightarrow{P(Xn)}& =a_1{A}^n\overrightarrow{v_1}+a_2{A}^n\overrightarrow{v_2}+...+a_n{A}^n\overrightarrow{v_n}\\ \stackrel{A\overrightarrow{v}=\lambda \overrightarrow{v}}{\to }\overrightarrow{P(Xn)}& =a_1{\lambda }_1^k\overrightarrow{v_1}+a_2{\lambda }_2^k\overrightarrow{v_2}+...+a_n{\lambda }_n^k\overrightarrow{v_n}\\ \stackrel{\forall {\lambda }_i⩽1}{\to }\overrightarrow{P(Xn)}& =0+0+\cdots +a_k{1}^n\overrightarrow{v_k}+\cdots +0\end{array}$$
因此最终随着一步步的计算，事件的概率总会收敛到一个定值，就是稳态。

那么高潮来了，如果在$n$步之后到达稳态，也就是说从n步之后不管怎样都一定符合稳态概率$P$，那如果我们就从第$n$步开始按着符合马尔科夫矩阵的状态转移概率的方式一步一步采样的话，不就能保证，所有采到的样本一定符合$P$分布么，妙啊。

这就是MCMC方法的核心思想。

细致平稳条件

思想归思想，怎么实现是另一个问题， 按照之前的思路，先有矩阵A，求出对应稳态概率，但是现在我们有了稳态概率$\pi (x)$,难道要反求出矩阵A，然后用A对样本进行迭代？够呛，因为往往$\pi (x)$是连续函数，咋反求，反正不现实。

所以现在是办法是构建一个转移概率函数$k(x_t+1|x_t)$,使得样本以此为规律迭代时，迭代出的每一个杨蓓都符合$\pi (x)$分布，还是符合MCMC的思想，只不过跟状态转移矩阵没关系了。

首先根据边缘概率公式
$${\pi }_{t+1}(x_{t+1})={\int }_{x_t}{\pi }_t(x_t)\cdot k(x_{t+1}|x_t)d{x}_t$$

类比
$$p(y)={\int }_{x}p(y|x)p(x)dx$$

因为要上一步和下一步采样的值符合同样的概率分布，因此${\pi }_{t+1}={\pi }_t$，得到
$$\pi (x_{t+1})={\int }_{x_t}\pi (x_t)\cdot k(x_{t+1}|x_t)d{x}_t\text{\hspace{1em}(1)}$$

这样我们就得到了转移概率$k$和目标概率$\pi$之间应该有的关系，也就是说我们要找到合适的$k$分布使得$k$与$\pi$符合上述公式，但是问题是$k$我好找，证明符合上述公式太费劲，因此我们需要一个另外的好证明的条件，只要满足那个条件，就可以说符合马尔科夫稳态，那就是detailed balance

detailed balance 定理
当$k$与$\pi$满足公式
$$\pi (x)k(x^{\ast }|x)=\pi (x^{\ast })k(x|x\ast )\text{\hspace{1em}(2)}$$
时，说明此时$\pi$分布是$k$转移概率的稳态分布

也就是说，对于任意两个样本A，B，A本身概率乘以从A跳到B的概率等于B本身的概率乘以从B跳到A的概率。那么这个定理有什么意义？

意义在意，只要满足公式（2），就能满足公式（1）。

证明
$$\begin{array}{rl}{\int }_{x_t}\pi (x_t)\cdot k(x_{t+1}|x_t)d{x}_t& \stackrel{带入式(2)}{\to }{\int }_{x_t}\pi (x_{t+1})\cdot k(x_t|x_{t+1})d{x}_t\\ & =\pi (x_{t+1})\cdot {\int }_{x_t}k(x_t|x_{t+1})d{x}_t\\ & =\pi (x_{t+1})\cdot 1\end{array}$$

妙啊。。。
但这是一个充分非必要条件，反过来不能证明。

有了以上的定理我们就可以去大胆的寻找合适的$k$分布，不同的$k$，即转移概率的分布就代表了不同的算法。其中常用的为M-H算法和(Gibbs Sampling)吉布斯采样算法

Metropolis-Hasting (M-H) 算法

1. 此刻有样本$x_t$
2. 生成$u\sim U(0,1)$,根据$q(x_{t+1}|x_t)$生成$x_{t+1}$, 此处$q$随便取啥都行
3. 判断，若$u<\min \left(1,\frac{\pi (x_{t+1})k(x_t\mid x_{t+1})}{\pi (x_t)k(x_t\mid x_{t+1})}\right)$则接受$x_{t+1}$，否则$x_{t+1}=x_t$

注意此处的$q$并不是代表符合了detailed balance 的一个转移概率,不是之前说的$k$，只是一个可以随便取的转移概率分布，需要满足detailed balance的不是$q$,而是$\min \left(1,\frac{\pi (x_{t+1})k(x_t\mid x_{t+1})}{\pi (x_t)k(x_t\mid x_{t+1})}\right)$,这才是$k$,那么它满不满足呢？

证明
$$\begin{array}{rl}\pi (x_t)\cdot k(x_{t+1}|x_t)& =\pi (x_t)\cdot q(x_{t+1}|x_t)\cdot \min \left(1,\frac{\pi (x_{t+1})q(x_t\mid x_{t+1})}{\pi (x_t)q(x_{t+1}\mid x_t)}\right)\\ & =\min \left(\pi (x_t)\cdot q(x_{t+1}|x_t),\pi (x_{t+1})q(x_t\mid x_{t+1})\right)\\ & =\pi (x_{t+1})\cdot q(x_t|x_{t+1})\cdot \min \left(\frac{\pi (x_t)q(x_{t+1}\mid x_t)}{\pi (x_{t+1})q(x_t\mid x_{t+1})},1\right)\\ & =\pi (x_{t+1})\cdot k(x_t|x_{t+1})\end{array}$$

因此此时$\min \left(1,\frac{\pi (x_{t+1})k(x_t\mid x_{t+1})}{\pi (x_t)k(x_t\mid x_{t+1})}\right)$就是我们要的转移概率分布，它即与目标分布$\pi$有关，又符合detailed balance分布，妙啊~~~

需要注意的是，在MH算法里面如何选择$q$是很重要的，虽然$q$可以随便选，但是算法中是靠着$\frac{\pi (x_t)q(x_{t+1}\mid x_t)}{\pi (x_{t+1})q(x_t\mid x_{t+1})}$这个比例来进行数据筛选，因此$q$当然还是选的越接近$\pi$分布越好。

Gibbs Sampling-吉布斯采样算法

之前提到的采样办法，样本都是一维的，就是都沿着$x$轴采样就好了，那如果样本是高维的怎么办，这就要用到Gibbs Sampling。

对于一个符合联合概率分布$p(x,y,z)$的分布，我们想去采集样本 { ( x 1 , y 1 , z 1 ) , ( x 2 , y 2 , z 2 ) , ( x 3 , y 3 , z 3 ) , ⋯   , ( x n , y n , z n ) } \{ (x_1,y_1,z_1), (x_2,y_2,z_2), (x_3,y_3,z_3),\cdots ,(x_n,y_n,z_n) \} {(x1​,y1​,z1​),(x2​,y2​,z2​),(x3​,y3​,z3​),⋯,(xn​,yn​,zn​)},那吉布斯采样的办法就是每次只沿着某一条坐标轴采样，在采样这条坐标轴的时候，其他坐标轴上的样本值不变。
$$\begin{array}{rl}{x}_2\sim & p(x|y_1,z_1)\\ y_2\sim & p(y|x_2,z_1)\\ z_2\sim & p(z|x_2,y_2)\\ & ⋮\\ x_n\sim & p(x|y_{n-1},z_{n-1})\\ y_n\sim & p(y|x_n,z_{n-1})\\ z_n\sim & p(z|x_n,y_n)\end{array}$$
那既然吉布斯采样时MCMC的一种，那必然要先证明他符合detailed balance条件，那就直观一点用二维的分布去证明。
按照吉布斯采样的法则，每步只对一个维度采样，那可以先选$x$维度，即$y$维度保持不变，也就是从$A(x_1,y_1)$跳到$C(x_2,y_1)$。
$$\begin{array}{r}p(x_1,y_1)\cdot p(x_2|y_1)=p(x_1|y_1)\cdot p(y_1)\cdot p(x_2|y_1)\\ p(x_2,y_1)\cdot p(x_1|y_1)=p(x_2|y_1)\cdot p(y_1)\cdot p(x_1|y_1)\end{array}$$
上面等式右边完全一样，那么
$$p(x_1,y_1)\cdot p(x_2|y_1)=p(x_2,y_1)\cdot p(x_1|y_1)$$
由此可以证出这种采样符合detailed balance定理，可以保证这样一步步走下去，最终样本会收敛到稳态概率$p(x,y)$

或者从另一个角度去证明，从一个更泛化的角度去证明。
对于联合分布$p(x_1,x_2,x_3,\cdots ,x_n)$,定义 X − i = { x 1 , x 2 , ⋯   , x i − 1 , x i + 1 , ⋯   , x n } X_{-i}=\{ x_1, x_2,\cdots, x_{i-1}, x_{i+1}, \cdots, x_n \} X−i​={x1​,x2​,⋯,xi−1​,xi+1​,⋯,xn​}，就是排除$x_i$那一项的集合。那么先来计算一个东西$\frac{p(x^{\ast })p(x_i|X_{-i}^{\ast })}{p(x)p(x_i^{\ast }|X_{-i})}$
$$\begin{array}{rl}\frac{p(x^{\ast })p(x_i|X_{-i}^{\ast })}{p(x)p(x_i^{\ast }|X_{-i})}& =\frac{p(x^{\ast }|X_{-i}^{\ast })p(X_{-i}^{\ast })p(x_i|X_{-i}^{\ast })}{p(x|X_{-i})p(X_{-i})p(x_i^{\ast }|X_{-i})}\end{array}$$
这一步$p(x^{\ast })$替换为$p(x^{\ast }|X_{-i}^{\ast })p(X_{-i}^{\ast })$怎么来的，$p(x^{\ast }|X_{-i}^{\ast })p(X_{-i}^{\ast })=p(x^{\ast },X_{-i}^{\ast })$才对啊，对这个联合分布取积分才会得到边缘分布$p(x^{\ast })$啊。。。
因为在吉布斯采样中，每次改变的只是一个维度上的量，其他维度上的值是固定的，因此此时的联合分布因为其他维度固定住了，坍缩成了一维的，就算是对其他维度求积分，其他维度的可能取值也只有一个，就是此刻的$X_{-i}^{\ast }$。

接着上面的式子往下算，刚说了其他维度的值固定，那么也就是说$X_{-i}^{\ast }=X_{-i}$,那么接着往下算上面的公式。
$$\frac{p(x^{\ast }|X_{-i}^{\ast })p(X_{-i}^{\ast })p(x_i|X_{-i}^{\ast })}{p(x|X_{-i})p(X_{-i})p(x_i^{\ast }|X_{-i})}=\frac{p(x^{\ast }|X_{-i})p(X_{-i})p(x_i|X_{-i})}{p(x|X_{-i})p(X_{-i})p(x_i^{\ast }|X_{-i})}=1$$
算这个等于1有啥意义？现在我们来对比一下H-M 算法中的那个算子$\min \left(1,\frac{\pi (x_{t+1})k(x_t\mid x_{t+1})}{\pi (x_t)k(x_t\mid x_{t+1})}\right)$和$\frac{p(x^{\ast })p(x_i|X_{-i}^{\ast })}{p(x)p(x_i^{\ast }|X_{-i})}$。$p(x^{\ast })$不就是$\pi (x_{t+1})$,就是在这个选中的维度下的单目标概率分布么。$p(x_i|X_{-i}^{\ast })$不就是$k(x_t\mid x_{t+1})$，就是别的维度都固定只在这个维度下的转移概率么。那最后$\min \left(1,1\right)$不就意味着所有新采的样本我都接受么。

也就是说，吉布斯采样是一种特殊的H-M采样，新样本接受率为100%的H-M采样。
妙啊#%￥&*@

Slice Sampling

超纲了，日后再说。。。