MCMC与Gibbs采样

一、均匀随机数生成算法

平方取中法
将2s位十进制数字 $x_n$ 平方后取中间的2s位数字：
$\begin{split} x_{n+1}=\frac{x_n^2}{10^s}\mod{10^{2s}} \\ r_n=\frac{x_n}{10^{2s}} \sim Uniform(0,1) \end{split}$
但其长度较短且均匀性较差，针对性的改进方法是乘积取中法：
$\begin{split} x_{n+2}=\frac{x_{n+1}x_n}{10^s}\mod{10^{2s}} \\ r_n=\frac{x_n}{10^{2s}} \sim Uniform(0,1) \end{split}$
但是随机数长度仍然不够，且对初值依赖比较大。
移位法
$x_{n+1}=x_n2^7+x_n2^{-7} \mod{2^{32}}$
混合线性同余法
$\begin{split} x_{n+1}=\lambda x_n+c \mod{M} \\ r_n=\frac{x_n}{M} \sim Uniform(0,1) \end{split}$
周期达到M的充要条件为：
$\begin{split} &1) c与M互素 \\ &2) \lambda-1为M的每一个素因子p的倍数 \\ &3) 若M为4的倍数，则\lambda-1为4的倍数 \end{split}$
Box-Muller变换
如果相互独立的两个随机变量 $U_1,U_2\sim Uniform(0,1)$ ，则满足 $Z_0=\sqrt{-2\ln{U_1}}\cos(2 \pi U_2) \\ Z_1=\sqrt{-2\ln{U_1}}\sin(2 \pi U_2)$ 的 $Z_0,Z_1\sim N(0,1)$
其他一些分布的随机数生成方法亦可以通过类似的变换得到。

二、markov链

马尔科夫性
下一个状态只依赖于当前状态，而与之前所有状态无关
$P(X_{t+1}=x|X_t, X_{t-1}, \cdots) =P(X_{t+1}=x|X_t)$
马氏定理
对于非周期且任意两点连通的马氏链，其概率转移矩阵 $P=P[P_{ij}]$ 的极限式子 $\lim_{n\to\infty}P_{ij}^n=\pi(j)$ 存在，其中 $\pi$ 称为马氏链的平稳分布，并有
$\begin{split} &1. \qquad \lim_{n \rightarrow \infty} P^n =\begin{bmatrix} \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \end{bmatrix} \\ &2. \qquad \pi(j) = \sum_{i=0}^{\infty}\pi(i)P_{ij} \\ &3. \qquad \pi是\pi P=\pi的唯一非负解 \end{split}$
细致平稳条件
如果非周期马氏链的转移矩阵 $P$ 和分布 $\pi(x)$ 满足
$\pi(i)P_{ij} = \pi(j)P_{ji} \quad\quad \text{for all} \quad i,j$ 则 $\pi(x)$ 是马氏链的平稳分布，上式被称为细致平稳条件 (detailed balance condition)。

三、MCMC(Markov Chain Monte Carlo)

设随机变量 $X \sim p(x)$ ，将markov链的状态空间看做随机变量 $X$ 的所有可能取值。
因为，当马氏链转移步数足够大时，状态分布将收敛到平稳分布。
所以，如果能构造一个转移矩阵为 $P$ 的马氏链，使其平稳分布恰好等于随机变量的概率分布，则可依次生成符合该分布的样本（即每一步的状态值）。
假设已经有一个转移矩阵 $P$ ，通常 $p(i)P_{ij} \neq p(j)P_{ji}$ ，我们构造 $\alpha_{ij}=p(j)P_{ji}$ ，则有
$p(i)P_{ij}\alpha_{ij} = p(j)P_{ji}\alpha_{ji}$ 从而，令 $Q_{ij}=P_{ij}p(j)P_{ji}$ 为状态转移矩阵即可。

MCMC采样算法

因为 $\alpha$ 总是小于1，算法速度比较慢，我们可以令 $\alpha_{ij}= \min\left\{\frac{p(j)P_{ji}}{p(i)P_{ij}},1\right\}$ 来改善算法速度，这个算法称为Metropolis-Hastings算法。

Metropolis-Hastings算法

四、Gibbs Sampling

初始转移矩阵 $P$ 可以人为指定，故而存在改进空间。
对于高纬空间， $\alpha$ 的存在也使得采样算法效率不高。
考察 $x$ 坐标相同的两个点 $A(x_1,y_1),B(x_1,y_2)$ ，我们发现:
$\begin{split} p(x_1,y_1)p(y_2|x_1)=p(x_1)p(y_1|x_1)p(y_2|x_1) \\ p(x_1,y_2)p(y_1|x_1)=p(x_1)p(y_2|x_1)p(y_1|x_1) \end{split}$ 从而有 $p(x_1,y_1)p(y_2|x_1)=p(x_1,y_2)p(y_1|x_1)$ 即 $p(A)p(y_2|x_1)=p(B)p(y_1|x_1)$

于是我们可以如下构造平面上任意两点之间的转移概率矩阵Q：
$\begin{align*} Q(A\rightarrow B) & = p(y_B|x_1) & \text{如果} \quad x_A=x_B=x_1 & \\ Q(A\rightarrow C) & = p(x_C|y_1) & \text{如果} \quad y_A=y_C=y_1 & \\ Q(A\rightarrow D) & = 0 & \text{其它} & \end{align*}$
从而有如下采样算法：