Metropolis-Hastings 采样原理
我们的⽬标采样分布是$p(z)$，同时我们⼿⾥有⼀个便于随时间进⾏遍历的⻢尔科夫链，它的状态转移矩阵为$Q$。
为了⽅便地在⻢尔科夫链上随时间进⾏状态转移，这⾥的矩阵 Q 设计相当讲究：
\[Q_{ij} = P(x^* \mid x) = N(\mu=x, \sigma=1).pmf(x^*)\]它指的是从$t$时刻的某个指定采样点$x$转移到$t+1$时刻的指定采样点$x^$的概率等于在均值为$x$，⽅差为 1 的整体分布中取得$x^$的概率值，感觉其实应该写成$pdf$的，因为正态分布是⼀个连续型的分布，但是这里因为是计算机采样，实际上我们已经把⽬标分布以及这⾥的正态分布都给离散化了，所以写成$pmf$更符合实际操作的需要。
但是就这么⼀个⽬标概率$p(z)$和⼀个⼋竿⼦打不着的我们⼼中合适的提议矩阵$Q$不能满⻢尔科夫链的细致平稳条件， 即：
\[p(z)Q(z, z^*) \neq p(z^*)Q(z^*, z)\]为了后面整个表达的格式统一， 我们将该状态转移矩阵对应的马尔科夫链的状态$z$转移到状态$z^$的概率记作$Q(z, z^)$。
为了让这个细致平稳条件的等式满足， 我们引入一个接受率因子的表达式：
\[\alpha(z, z^*) = \min(1,\frac{p(z^*)Q(z^*, z)}{p(z)Q(z, z^*)})\]原来的不等式两边各自乘上对应的接受率因子$\alpha$， 细致平稳条件就可以神奇地得到满足：
\[p(z)Q(z, z^*)\alpha(z, z^*) = p(z^*)Q(z^*, z)\alpha(z^*, z)\]为什么此处细致平稳条件可以被满足：
\[\begin{aligned} p(z)Q(z, z^*)\alpha(z, z^*) & = p(z)Q(z, z^*)\min(1, \frac{p(z^*)Q(z^*, z)}{p(z)Q(z, z^*)}) \\ & = \min(p(z)Q(z, z^*), p(z^*)Q(z^*, z)) \\ & = p(z^*)Q(z^*, z) \min(1, \frac{p(z)Q(z, z^*)}{p(z^*)Q(z^*, z)}) \\ & = p(z^*)Q(z^*, z)\alpha(z^*, z) \end{aligned}\]此时我们证明了细致平稳条件的成立。
但是此时$p(z)$并不是$Q$对应的马尔科夫链的平稳分布， 而是跌价了接受率$\alpha$之新的转移过程。
M-H 采样步骤
如果我们要采样$N$个样本， 那么第$i$轮的采样过程如下：
- 第一步： 我们从$[0, 1]$均匀分布中随机抽取一个数$u$， 即$u \thicksim U(0, 1)$。
- 第二步： 我们利用状态转移概率矩阵，将上一轮的样本$z^{(i - 1)}$依概率抽取出这一步应该转移到的状态$z^$， 即： $z^ \thicksim Q(z \mid z^{(i-1)})$。
- 第三步：计算出接受率$\alpha = \min(1, \frac{p(z^)Q(z^, z)}{p(z)Q(z, z^*)})$。
- 第四步：判断， 如果此时$u \leq \alpha$， 那么我们把$z^$作为第$i$轮的采样值， 即$z^{(i)} = z^$， 如果$u \geq \alpha$，抛弃$z^*$， 我们沿用上一轮的采样值$z^{i-1}$作为本轮的采样值， 即$z^{(i)} = z^{(i-1)}$。
经过燃烧期之后，我们取连续$N$个时间点的状态值， 就可以作为我们目标分布的样本近似了。
⼀般而⾔，我们会选择⼀个⽐较适合进⾏状态转移操作的⻢尔科夫链，这⾥描述为：作为状态空间$S$⾥的任意两个状态，第$i$个状态$x$，第$j$个状态为$x^$，$Q_{ij} = P(x^ \mid x) = N(\mu = x, \sigma=1).pmf(x^*)$
满足这个条件的$Q$矩阵是一个非常好的选择：
- 第一， 它是合理的， 矩阵的任意一行之和显然为1， 因此满足状态转移矩阵的条件。
就是在一个以$x$为均值， 1为方差的正态分布中， 状态空间中所有可选状态的概率之和为1， 这是概率分布的定义决定的。
-
第二， 它是方便的， 当前状态为$x$的情况下， 我们决定下一个状态$x^*$的过程， 就是从正态分布$N(\mu=x, \sigma=1)$中采得一个样本的过程。
-
第三， 它是对称的， 即$Q_{ij} = Q_{ji}$， 也就是说$N(\mu=x, \sigma=1).pmf(x^) = N(\mu=x^, \sigma=1).pmf(x)$， 这样使得$Q(z^, z) = Q(z, z^)$， 接受率得以简化：
这就是 Metropolis-Hastings 算法的最核心的部分。
Gibbs ⽅法核⼼流程
吉布斯采样是⼀种针对⾼维分布的采样⽅法，假设我们待采样的⾼维⽬标分布为$p(z_1, z_2, z_3, …, z_m)$， 吉布斯采样的⽬标就是从这个⾼维分布中采样出⼀组样本$z^{(1)}, z^{(2)}, z^{(3)}, …, z^{(N)}$， 使得这一组样本能够服从$m$维分布$p(z)$。
我们先约定⼀个记法：按照上述假设，从分布中采得的样本$z$是⼀个拥有$m$维属性值的样本值， $z_i$表⽰某个样本$z$的第$i$维属性， $z_{-i}$表⽰除开第$i$维样本外，剩余的$m-1$维属性： $z_1, z_2, …, z_{i-1}, z_{i+1}, …, z_m$。
具体如何采样，核⼼在七个字：⼀维⼀维地采样。
我们以三维随机变量$z$为例：
\[p(z) = p(z_1, z_2, z_3)\]⾸先在 0 时刻，我们给定⼀个初值，也就是$z^{(0)}$的三个属性的值：
\[z_1^{(0)}, z_2^{(0)}, z_3^{(0)}\]此时我们来看下⼀个时刻，也就是时刻 1 的采样值如何⽣成：
\[z_1^{(1)} \thicksim p(z_1 \mid z_2^{(0)}, z_3^{(0)})\] \[z_2^{(1)} \thicksim p(z_2 \mid z_1^{(1)}, z_3^{(0)})\] \[z_3^{(1)} \thicksim p(z_3 \mid z_1^{(1)}, z_2^{(1)})\]这样就获得了时刻1的包含三维属性值的完整采样值。
那么接着往下递推，概况起来，通过$t$时刻的采样值递推到$t+1$时刻，也是如上述⽅法那样⼀维⼀维的⽣成采样值的各个属性：固定其他维的属性值，通过条件概率，依概率⽣成$t+1$时刻的某⼀维度的属性值。
\[z_1^{(t+1)} \thicksim p(z_1 \mid z_2^{(t)}, z_3^{(t)})\] \[z_2^{(t+1)} \thicksim p(z_2 \mid z_1^{(t+1)}, z_3^{(t)})\] \[z_3^{(t+1)} \thicksim p(z_3 \mid z_1^{(t+1)}, z_2^{(t+1)})\]这样就由$t$时刻的采样值$z^{(t)}$得到了下一时刻$t+1$的采样值$z^{(t+1)}$。
按照该方法循环下去， 就可以采样出$N$个样本值： $z^{(1)}, z^{(2)}, z^{(3)}, …, z^{(N)}$， 丢弃掉前面一部分燃烧期的样本， 剩下的样本就服从我们的目标高维分布。
Gibbs 采样的合理性
为什么通过这种采样⽅法采得的样本能够服从⽬标分布$p(z)$？
我们拿 Gibbs 采样和 M-H ⽅法进⾏类⽐，我们会发现其实 Gibbs 采样就是 M-H 采样的⼀种特殊情况，这个结论就可以证 明 Gibbs 采样的可⾏性，注意逻辑是这样的，我们已经证明了 M-H 是符合细致平稳条件的，而 Gibbs 是 M-H 的⼀种特殊情况，因 此 Gibbs 也符合细致平稳条件，因此可⾏。
⾸先，我们采样的⽬标分布是$p(z)$，而我们实际使⽤的提议转移过程$Q$就是概率$p$本⾝： $Q(z^, z) = p(z_i \mid z^_{-i})$，在这个采样过程，好像没有接受率 的出现。
我们重新看接受率$\alpha(z^*, z)$的定义式：
\[\alpha(z^*, z) = \min(1, \frac{p(z^*)Q(z^*, z)}{p(z)Q(z, z^*)})\]重点就在$\frac{p(z^)Q(z^, z)}{p(z)Q(z, z^*)}$这个式子：
\[p(z) = p(z, z_{-i}) = p(z_i \mid z_{-i})p(z_{-i})\] \[\frac{p(z^*)Q(z^*, z)}{p(z)Q(z, z^*)} = \frac{[p(z_i^* \mid z_{-i}^*)p(z_{-i}^*)]p(z_i \mid z_{-i}^*)}{[p(z_i \mid z_{-i})p(z_{-i})]p(z_i^* \mid z_{-i})}\]在上面$z \rightarrow z^$当中， 我们通过固定采样值$z$的$-i$维属性， 单采$z^$的第$i$维属性值， 那么在一次属性采样前后$z_{-i}$和$z_{-i}^*$是一码事儿。
什么意思？
那么我们把上面式子中所有的$z_{-i}^*$都替换成$z_{-i}$， 就有：
\[\frac{p(z^*)Q(z^*, z)}{p(z)Q(z, z^*)} = \frac{[p(z_i^* \mid z_{-i}^*)p(z_{-i}^*)]p(z_i \mid z_{-i}^*)}{[p(z_i \mid z_{-i})p(z_{-i})]p(z_i^* \mid z_{-i})} = \frac{[p(z_i^* \mid z_{-i})p(z_{-i})]p(z_i \mid z_{-i})}{[p(z_i \mid z_{-i})p(z_{-i})]p(z_i^* \mid z_{-i})}=1\]也就是说Gibbs采样中的接受率$\alpha$为：
\[\alpha(z^*, z) = \min(1, \frac{p(z^*)Q(z^*, z)}{p(z)Q(z, z^*)}) = \min(1, 1) = 1\]因此满足细致平稳条件的原始等式：
\[p(z)Q(z, z^*)\alpha(z, z^*) = p(z^*)Q(z^*, z)\alpha(z^*, z)\]就等效为：
\[p(z)Q(z, z^*)*1 = p(z^*)Q(z^*, z) * 1 \Rightarrow p(z)Q(z, z^*) = p(z^*)Q(z^*, z)\]其中：$Q(z^, z) = p(z_i \mid z_{-i}^)$
因此$p(z)$和$Q(z , z^*)$满足细致平稳条件， 按照$Q$描述的状态转移过程所得到的采样样本集， 在数量足够多的情况，就近似为目标分布$p(z)$。
提议转移概率矩阵$Q$就定位了$p$的满条件概率， 那么很显然就有一个限制条件， 那就是$p$的所有满条件概率都是易于采样的。
Gibbs 采样代码试验
⽐如⽬标分布是⼀个⼆维的⾼斯分布，它的均值向量：$\mu = \begin{bmatrix}3 \ 3\end{bmatrix}$， 协方差矩阵$\Sigma = \begin{bmatrix}1 & 0.6 \ 0.6 & 1\end{bmatrix}$， 即$z$包含两维属性$(z_1, z_2), z \thicksim N(\mu, \Sigma)$。
我们通过吉布斯采样的⽅法，来获取服从上述分布的⼀组样本值：
$z$包含两维属性$(z_1, z_2), z \thicksim N(\mu, \Sigma)$， 其中$\mu = \begin{bmatrix}\mu_1 \ \mu_2\end{bmatrix}， \Sigma = \begin{bmatrix}\delta_{11} & \delta_{12} \ \delta_{21} & \delta_{22}\end{bmatrix}$。
那么条件分布$z_1 \mid z_2$和$z_2 \mid z_1$分别服从下面的形式：
\[z_1 \mid z_2 \thicksim N(\mu_1 + \frac{\delta_{12}}{\delta_{22}}(z_2 - \mu_2), \delta_{11} - \frac{\delta_{12}^2}{\delta_{22}})\] \[z_2 \mid z_1 \thicksim N(\mu_2 + \frac{\delta_{12}}{\delta_{11}}(z_1 - \mu_2), \delta_{22} - \frac{\delta_{12}^2}{\delta_{11}})\]这就是二维高斯分布采样过程中，逐维采样所依从的满条件分布， 它是高斯分布， 因此操作尤其简便。
