​一、蒙特卡洛方法

1. 蒙特卡洛方法Monte Carlo 可以通过采用随机投点法来求解不规则图形的面积。 求解结果并不是一个精确值，而是一个近似值。当投点的数量越来越大时，该近似值也越接近真实值。
2. 蒙特卡洛方法也可以用于根据概率分布来随机采样的任务。

1.1 布丰投针问题

1. 布丰投针问题是1777年法国科学家布丰提出的一种计算圆周率的方法：随机投针法。其步骤为：
2. 首先取一张白纸，在上面绘制许多条间距为 d 的平行线。
3. 取一根长度为 l,l\lt d 的针，随机地向纸上投掷 n 次，观测针与直线相交的次数，记做 m 。
4. 计算针与直线相交的概率 p=\frac{m}{n} 。可以证明这个概率 p=\frac{2l}{\pi\times d} 。因此有： \pi = 2 \frac{n\times l}{m\times d}

• 由于向纸上投针是完全随机的，因此用二维随机变量 (X,Y) 来确定针在纸上的具体位置。其中：
• X 表示针的中点到平行线的距离，它是 [0,d/2] 之间的均匀分布。
• Y 表示针与平行线的夹角，它是 [0,\frac{\pi}{2}] 之间的均匀分布。

​

​

当 X\lt \frac{l}{2} \sin Y 时，针与直线相交。

由于 X,Y 相互独立，因此有概率密度函数：

p(X=x,Y=y) = \begin{cases} \frac{4}{\pi d},& 0\le x\le d/2,0\le y\le \pi/2\\ 0,& \text{else} \end{cases}

因此，针与直线相交的概率为：

P\{X\lt \frac l2\sin Y\}=\int\int\_{x\lt \frac l2 \sin y} p(x,y) dxdy=\int\_{x=0}^{x=\frac l2\sin y}\int\_{y=0}^{y= \pi/ 2} \frac {4}{\pi d} dxdy=\frac{2l}{\pi d}

根据 \frac{2l}{\pi d}=\frac{m}{n} 即可得证。

• 布丰投针问题中，蒙特卡洛方法是利用随机投点法来求解面积 \int\int\_{x\lt \frac l2 \sin y} p(x,y) dxdy 。因为曲线的积分就是面积，这里的曲线就是概率密度函数 p(X,Y) 。

1.2 蒙特卡洛积分

1. 对于函数 f(x)，其在区间 [a,b] 上的积分 \int\_a^b f(x) dx 可以采用两种方法来求解：投点法、期望法。
2. 投点法求积分：对函数 f(x)，对其求积分等价于求它的曲线下方的面积。 此时定义一个常数 M，使得 M\gt \max\_{a \le x\le b} f(x) ，该常数在区间 [a,b] 上的面积就是矩形面积 M(b-a) 。 随机向矩形框中随机的、均匀的投点，设落在函数 f(x) 下方的点为绿色，落在 f(x) 和 M 之间的点为红色。则有：落在 f(x) 下方的点的概率等于 f(x) 的面积比上矩形框的面积 。 具体做法是：从 [a,b] 之间的均匀分布中采样 x\_0 ，从 [0,M] 之见的均匀分布中采样 y\_0 ， (x\_0,y\_0) 构成一个随机点。
3. 若 y\_0 \le f(x\_0)，则说明该随机点在函数 f(x) 下方，染成绿色。
4. 若 f(x\_0) \lt y\_0 \le M，则说明该随机点在函数 f(x) 上方，染成红色。

假设绿色点有 n\_1 个，红色点有 n\_2 个，总的点数为 n\_1+n\_2 ，因此有：\int\_a^b f(x)dx = \frac{n\_1}{n\_1+n\_2}\times M(b-a) 。

​

​

1. 期望法求积分：假设需要求解积分 I=\int\_a^b f(x)dx ，则任意选择一个概率密度函数 p(x)，其中 p(x) 满足条件： \int\_a^b p(x) dx = 1\\ \text{if} \;f(x)\ne 0\;\text{then}\; p(x)\ne 0 ,\quad a\le x\le b 令： f^\*(x)=\begin{cases} \frac{f(x)}{p(x)},& p(x)\ne 0\\ 0,& p(x)=0 \end{cases} 则有：I=\int\_a^b f(x)dx=\int\_a^b f^\*(x)g(x)dx ，它刚好是一个期望：设随机变量 X 服从分布 X\sim p(x)，则 I=\mathbb E\_{X\sim p}[f^\*(X)] 。 则期望法求积分的步骤是：
2. 任选一个满足条件的概率分布 p(x) 。
3. 根据 p(x)，生成一组服从分布 p(x) 的随机数 x\_1,x\_2,\cdots,x\_N 。
4. 计算均值 \bar I=\frac 1N\sum\_{i=1}^N f^\*(x\_i) ，并将 \bar I 作为 I 的近似。
5. 在期望法求积分中，如果 a,b 均为有限值，则 p(x) 可以取均匀分布的概率密度函数： p(x) = \begin{cases} \frac {1}{b-a},& a\le x\le b\\ 0,& \text{else} \end{cases} 此时 f^\*(x) = (b-a) f(x)， \bar I = \frac{b-a}{N}\sum\_{i=1}^Nf(x\_i) 。 其物理意义为：\frac{\sum\_{i=1}^Nf(x\_i)}{N} 为在区间 [a,b] 上函数的平均高度，乘以区间宽度 b-a 就是平均面积。

​

1. 对于期望 \mathbb E\_p[f(X)]，如果 p(x) 或者 f(x) 的表达式比较复杂，则也可以转化为另一个期望的计算。 选择一个比较简单的概率密度函数 q(x)，根据： \mathbb E\_p[f(X)]=\int f(x) p(x)dx=\int f(x)\frac {p(x)}{q(x)}q(x)dx 令 f^\*(x) = f(x)\frac {p(x)}{q(x)}，则原始期望转换为求另一个期望 \mathbb E\_q[f^\*(X)] 。此时可以使用期望法求积分的策略计算。

1.3 蒙特卡洛采样

1. 采样问题的主要任务是：根据概率分布 p(x) ，生成一组服从分布 p(x) 的随机数 x\_1,x\_2,\cdots 。
2. 如果 p(x) 就是均匀分布，则均匀分布的采样非常简单。
3. 如果 p(x) 是非均匀分布，则可以通过均匀分布的采样来实现。其步骤是：
4. 首先根据均匀分布 U(0,1) 随机生成一个样本 z\_i 。
5. 设 \tilde P(x) 为概率分布 p(x) 的累计分布函数：\tilde P(x)=\int\_{-\infty}^x p(z) dz 。 令 z\_i=\tilde P(x\_i) ，计算得到 x\_i=\tilde P^{-1}(z\_i) ，其中 \tilde P^{-1} 为反函数，则 x\_i 为对 p(x) 的采样。

​

​

1. 通过均匀分布的采样的原理：假设随机变量 Z,X 满足 Z=\tilde P(X)，则 X 的概率分布为： p\_Z(z) \frac{d}{d x} \tilde P(x) 因为 Z 是 [0,1] 上面的均匀分布，因此 p\_Z(z)=1 ； \tilde P(x) 为概率分布 p(x) 的累计分布函数，因此 \frac{d}{d x} \tilde P(x)=p\_X(x) 。因此上式刚好等于 p(x) ，即：x\_i 服从概率分布 p(x) 。 这其中有两个关键计算：
2. 根据 p(x)，计算累计分布函数 \tilde P(x)=\int\_{-\infty}^x p(z) dz 。
3. 根据 z=\tilde P(x) 计算反函数 x=\tilde P^{-1}(z) 。

如果累计分布函数无法计算，或者反函数难以求解，则该方法无法进行。

1. 对于复杂的概率分布 p(x) ，难以通过均匀分布来实现采样。此时可以使用接受-拒绝采样 策略。
2. 首先选定一个容易采样的概率分布 q(x) ，选择一个常数 k ，使得在定义域的所有位置都满足 p(x) \le k\times q(x) 。
3. 然后根据概率分布 q(x) 随机生成一个样本 x\_i 。
4. 计算 \alpha\_i =\frac{p(x\_i)}{kq(x\_i)} ，以概率 \alpha\_i 接受该样本。 具体做法是：根据均匀分布 U(0,1) 随机生成一个点 u\_i 。如果 u\_i \le \alpha\_i ，则接受该样本；否则拒绝该样本。 或者换一个做法：根据均匀分布 U(0,kq(x\_i)) 生成一个随机点，如果该点落在灰色区间（(p(x\_i),kq(x\_i)])）则拒绝该样本；如果该点落在白色区间（[0,p(x\_i)]）则接受该样本。

​

​

1. 接受-拒绝采样 在高维的情况下会出现两个问题：
2. 合适的 q 分布比较难以找到。
3. 难以确定一个合理的 k 值。

这两个问题会导致拒绝率很高，无效计算太多。

二、马尔可夫链

1. 马尔可夫链是满足马尔可夫性质的随机过程。 马尔可夫链 X\_1,X\_2,\cdots 描述了一个状态序列，其中每个状态值取决于前一个状态。X\_t 为随机变量，称为时刻 t 的状态，其取值范围称作状态空间。 马尔可夫链的数学定义为： P(X\_{t+1}\mid X\_t,X\_{t-1},\cdots,X\_1)=P(X\_{t+1}\mid X\_t) 。

2.1 马尔可夫链示例

1. 社会学家把人按照经济状况分成三类：下层、中层、上层。用状态 1,2,3 代表着三个阶层。社会学家发现：决定一个人的收入阶层的最重要因素就是其父母的收入阶层。
2. 如果一个人的收入属于下层，则他的孩子属于下层的概率是 0.65，属于中层的概率是 0.28，属于上层的概率是 0.07 。
3. 如果一个人的收入属于中层，则他的孩子属于下层的概率是 0.15，属于中层的概率是 0.67，属于上层的概率是 0.18 。
4. 如果一个人的收入属于上层，则他的孩子属于下层的概率是 0.12，属于中层的概率是 0.36，属于上层的概率是 0.52 。

从父代到子代，收入阶层的变化的转移概率如下：

1. 使用矩阵的表示方式，转移概率矩阵记作： \mathbf P= \begin{bmatrix} 0.65&0.28&0.07\\ 0.15&0.67&0.18\\ 0.12&0.36&0.52 \end{bmatrix} 假设当前这一代人在下层、中层、上层的人的比例是概率分布 \vec\pi\_0=(\pi\_0(1),\pi\_0(2),\pi\_0(3)) ，则：
2. 他们的子女在下层、中层、上层的人的概率分布是 \vec\pi\_1=(\pi\_1(1),\pi\_1(2),\pi\_1(3))=\vec\pi\_0 \mathbf P
3. 他们的孙子代的分布比例将是 \vec\pi\_2=(\pi\_2(1),\pi\_2(2),\pi\_2(3))=\vec\pi\_1 \mathbf P=\vec\pi\_0 \mathbf P^{2}
4. ....
5. 第 n 代子孙在下层、中层、上层的人的比例是 \vec\pi\_n=(\pi\_n(1),\pi\_n(2),\pi\_n(3))=\vec\pi\_{n-1} \mathbf P=\vec\pi\_0 \mathbf P^{n}
6. 假设初始概率分布为 \pi\_0=(0.72,0.19,0.09) ，给出前 14 代人的分布状况： xxxxxxxxxx  0 0.72 0.19 0.09  1 0.5073 0.3613 0.1314  2 0.399708 0.431419 0.168873  3 0.34478781 0.46176325 0.19344894  4 0.3165904368 0.4755635827 0.2078459805  5 0.302059838985 0.482097475693 0.215842685322  6 0.294554638933 0.485285430346 0.220159930721  7 0.290672521545 0.486874112293 0.222453366163  8 0.288662659788 0.487677173087 0.223660167125  9 0.28762152488 0.488086910874 0.224291564246  10 0.287082015513 0.488297220381 0.224620764107  11 0.286802384833 0.488405577077 0.22479203809  12 0.286657431274 0.488461538107 0.224881030619  13 0.286582284718 0.488490482311 0.22492723297  14 0.28654332537 0.488505466739 0.224951207891 可以看到从第 9 代开始，阶层分布就趋向于稳定不变。
7. 如果换一个初始概率分布为 \vec\pi\_0=(0.51,0.34,0.15) ，给出前 14 代人的分布状况： xxxxxxxxxx  0 0.51 0.34 0.15  1 0.4005 0.4246 0.1749  2 0.345003 0.459586 0.195411  3 0.31663917 0.47487142 0.20848941  4 0.3020649027 0.4818790066 0.2160560907  5 0.294550768629 0.48521729983 0.220231931541  6 0.290668426368 0.486853301457 0.222478272175  7 0.288659865019 0.487671049342 0.223669085639  8 0.28761985994 0.488085236095 0.224294903965  9 0.287081082851 0.488296834394 0.224622082755  10 0.286801878943 0.488405532034 0.224792589023  11 0.286657161801 0.488461564615 0.224881273584  12 0.286582142693 0.488490512087 0.224927345221  13 0.28654325099 0.488505487331 0.224951261679  14 0.286523087645 0.488513240994 0.224963671362 可以发现到第 8 代又收敛了。
8. 两次不同的初始概率分布，最终都收敛到概率分布 \vec\pi=(0.286 , 0.489, 0.225) 。 这说明收敛的行为和初始概率分布 \vec\pi\_0 无关，而是由概率转移矩阵 \mathbf P 决定的。 计算 \mathbf P^{n}： xxxxxxxxxx  0 [ 0.65  0.28  0.07 0.12  0.36  0.52]  1 [ 0.4729  0.3948  0.1323 0.1944  0.462   0.3436] ...  18 [ 0.28650397  0.48852059  0.22497545 0.28649994  0.48852213  0.22497793]  19 [ 0.28650272  0.48852106  0.22497622 0.28650063  0.48852187  0.2249775 ]  20 [ 0.28650207  0.48852131  0.22497661 0.28650099  0.48852173  0.22497728]  21 [ 0.28650174  0.48852144  0.22497682 0.28650118  0.48852166  0.22497717]  ... 可以看到 ： \mathbf P^{18}=\mathbf P^{19}=\cdots=\begin{bmatrix} 0.286& 0.489 & 0.225 \\ 0.286 & 0.489 & 0.225 \\ 0.286& 0.489 & 0.225 \\ \end{bmatrix} 发现当 n 足够大的时候， 矩阵 \mathbf P^{n} 收敛且每一行都稳定收敛到 \vec\pi=(0.286 , 0.489, 0.225) 。

2.2 平稳分布

1. 马尔可夫链定理：如果一个非周期马尔可夫链具有转移概率矩阵 \mathbf P，且它的任何两个状态是联通的，则有： \lim\_{n\rightarrow \infty} \mathbf P^{n}=\begin{bmatrix} \pi {(1)}&\pi {(2)}&\cdots&\pi {(j)}&\cdots\\ \pi {(1)}&\pi {(2)}&\cdots&\pi {(j)}&\cdots\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ \pi {(1)}&\pi {(2)}&\cdots&\pi {(j)}&\cdots\\ \vdots&\vdots&\ddots&\vdots&\vdots \end{bmatrix}\\ \pi {(j)} =\sum\_{i=0}^{\infty} \pi {(i)} P\_{i,j} 其中：
2. 1,2,\cdots,j,\cdots 为所有可能的状态。
3. P\_{i,j} 是转移概率矩阵 \mathbf P 的第 i 行第 j 列的元素，表示状态 i 转移到状态 j 的概率。
4. 概率分布 \vec\pi 是方程 \vec\pi\mathbf P= \vec\pi 的唯一解，其中 \vec\pi=(\pi {(1)},\pi {(2)},\cdots,\pi {(j)},\cdots),\sum\_{j=0}^{\infty}\pi {(j)} =1 。 称概率分布 \vec \pi 为马尔可夫链的平稳分布。
5. 注意，在马尔可夫链定理中：
6. 马尔可夫链的状态不要求有限，可以是无穷多个。
7. 非周期性在实际任务中都是满足的。
8. 两个状态的连通指的是：状态 i 可以通过有限的 n 步转移到达 j (并不要求从状态 i 可以直接一步转移到状态 j ）。 马尔可夫链的任何两个状态是联通的含义是：存在一个 n ，使得矩阵 \mathbf P^{n} 中的任何一个元素的数值都大于零。
9. 从初始概率分布 \vec\pi\_0 出发，在马尔可夫链上做状态转移，记时刻 i 的状态 X\_i 服从的概率分布为 \vec\pi\_i， 记作 X\_i \sim \vec \pi\_i ，则有： X\_0 \sim \vec\pi\_0\\ X\_1 \sim \vec\pi\_1\\ \cdots\\ X\_n \sim \vec\pi\_n\\ X\_{n+1} \sim \vec \pi\_{n+1}\\ \cdots 假设到达第 n 步时，马尔可夫链收敛，则有： X\_n \sim \vec\pi\\ X\_{n+1} \sim \vec\pi\\ \cdots 所以 X\_n,X\_{n+1},X\_{n+2},\cdots 都是同分布的随机变量（当然它们并不独立）。 如果从一个具体的初始状态 x\_0 开始，然后沿着马尔可夫链按照概率转移矩阵做调整，则得到一个转移序列 x\_0,x\_1,\cdots,x\_{n},x\_{n+1},\cdots 。 根据马尔可夫链的收敛行为，当 n 较大时， x\_{n},x\_{n+1},\cdots 将是平稳分布 \vec \pi 的样本。
10. 定理：如果非周期马尔可夫链的转移矩阵 \mathbf P 和某个分布 \vec \pi 满足：\pi(i)P\_{i,j}=\pi(j)P\_{j,i} ，则 \vec \pi 是马尔可夫链的平稳分布。 这被称作马尔可夫链的细致平稳条件detailed balance condition ，其证明如下： \pi(i)P\_{i,j}=\pi(j)P\_{j,i} \rightarrow \sum\_{i=1}^{\infty}\pi(i)P\_{i,j}=\sum\_{i=1}^{\infty}\pi(j)P\_{j,i}= \pi(j)\sum\_{i=1}^{\infty}P\_{j,i}=\pi(j) \rightarrow \vec\pi\mathbf P=\vec\pi .

三、MCMC 采样

1. 概率图模型中最常用的采样技术是马尔可夫链蒙特卡罗方法Markov Chain Monte Carlo:MCMC。 给定连续随机变量 X \in \mathcal X 的概率密度函数 \tilde p(x) ，则 X 在区间 \mathbb A 中的概率可以计算为： P(\mathbb A)=\int\_\mathbb A\tilde p(x)dx 如果函数 f: \mathcal X \longmapsto \mathbb R， 则可以计算 f(X) 的期望：\mathbb E\_{X \sim \tilde p(x)}[f(X)]=\int f(x)\tilde p(x)dx 。
2. 如果 X 不是一个单变量，而是一个高维的多元变量 \vec X ，且服从一个非常复杂的分布，则对于上式的求积分非常困难。为此，MCMC先构造出服从 \tilde p 分布的独立同分布随机变量 \mathbf {\vec x}\_1,\mathbf {\vec x}\_2,\cdots,\mathbf {\vec x}\_N ， 再得到 \mathbb E\_{\vec X\sim \tilde p(\mathbf {\vec x})}[f(\vec X)] 的无偏估计： \mathbb E\_{\vec X\sim \tilde p(\mathbf {\vec x})}[f(\vec X)]=\frac 1N \sum\_{i=1}^{N}f(\mathbf {\vec x}\_i)
3. 如果概率密度函数 \tilde p(\mathbf {\vec x}) 也很复杂，则构造服从 \tilde p 分布的独立同分布随机变量也很困难。MCMC 方法就是通过构造平稳分布为 \tilde p(\mathbf {\vec x}) 的马尔可夫链来产生样本。

3.1 MCMC 算法

1. MCMC 算法的基本思想是：先设法构造一条马尔可夫链，使其收敛到平稳分布恰好为 \tilde p 。然后通过这条马尔可夫链来产生符合 \tilde p 分布的样本。最后通过这些样本来进行估计。 这里马尔可夫链的构造至关重要，不同的构造方法将产生不同的MCMC算法。Metropolis-Hastings:MH算法是MCMC的重要代表。
2. 假设已经提供了一条马尔可夫链，其转移矩阵为 \mathbf Q 。目标是另一个马尔科夫链，使转移矩阵为 \mathbf P 、平稳分布是 \tilde p 。 通常 \tilde p(i)Q\_{i,j} \ne \tilde p(j)Q\_{j,i} ，即 \tilde p 并不满足细致平稳条件不成立。但是可以改造已有的马尔可夫链，使得细致平稳条件成立。 引入一个函数 \alpha(i,j) ，使其满足：\tilde p(i)Q\_{i,j}\alpha(i,j)=\tilde p(j)Q\_{j,i}\alpha(j,i) 。如：取 \alpha(i,j)=\tilde p(j)Q\_{j,i} ，则有： \tilde p(i)Q\_{i,j}\alpha(i,j)=\tilde p(i)Q\_{i,j}\tilde p(j)Q\_{j,i}=\tilde p(j)Q\_{j,i}\tilde p(i)Q\_{i,j}=\tilde p(j)Q\_{j,i}\alpha(j,i) 令： Q^{\prime}\_{i,j}=Q\_{i,j}\alpha(i,j),Q^{\prime}\_{j,i}=Q\_{j,i}\alpha(j,i) ，则有 \tilde p(i)Q^{\prime}\_{i,j}=\tilde p(j)Q^{\prime}\_{j,i} 。其中 Q^{\prime}\_{i,j} 构成了转移矩阵 \mathbf Q^{\prime} 。而 \mathbf Q^{\prime} 恰好满足细致平稳条件，因此它对应的马尔可夫链的平稳分布就是 \tilde p 。
3. 在改造 \mathbf Q 的过程中引入的 \alpha(i,j) 称作接受率。其物理意义为：在原来的马尔可夫链上，从状态 i 以 Q\_{i,j} 的概率跳转到状态 j 的时候，以 \alpha(i,j) 的概率接受这个转移。
4. 如果接受率 \alpha(i,j) 太小，则改造马尔可夫链过程中非常容易原地踏步，拒绝大量的跳转。这样使得马尔可夫链遍历所有的状态空间需要花费太长的时间，收敛到平稳分布 \tilde p 的速度太慢。
5. 根据推导 \alpha(i,j)=\tilde p(j)Q\_{j,i} ，如果将系数从 1 提高到 K ，则有： \alpha^{\*}(i,j)=K\tilde p(j)Q\_{j,i}=K\alpha (i,j)\\ \alpha^{\*}(j,i)=K\tilde p(i)Q\_{i,j}=K\alpha (j,i) 于是： \tilde p(i)Q\_{i,j}\alpha^{\*}(i,j)=K\tilde p(i)Q\_{i,j}\alpha(i,j)=K\tilde p(j)Q\_{j,i}\alpha(j,i)=\tilde p(j)Q\_{j,i}\alpha^{\*}(j,i) 。因此，即使提高了接受率，细致平稳条件仍然成立。
6. 将 \alpha(i,j),\alpha(j,i) 同比例放大，取：\alpha(i,j)=\min\{\frac{\tilde p(j)Q\_{j,i}}{\tilde p(i)Q\_{i,j}},1\} 。
7. 当 \tilde p(j)Q\_{j,i}= \tilde p(i)Q\_{i,j} 时： \alpha(i,j)=\alpha(j,i)=1，此时满足细致平稳条件。
8. 当 \tilde p(j)Q\_{j,i}\gt \tilde p(i)Q\_{i,j} 时： \alpha(i,j)=1,\alpha(j,i)=\frac{\tilde p(i)Q\_{i,j}}{\tilde p(j)Q\_{j,i}}，此时满足细致平稳条件。
9. 当 \tilde p(j)Q\_{j,i} \lt \tilde p(i)Q\_{i,j} 时： \alpha(i,j)=\frac{\tilde p(j)Q\_{j,i}}{\tilde p(i)Q\_{i,j}},\alpha(j,i)=1，此时满足细致平稳条件。
10. MH 算法：
11. 输入：
12. 先验转移概率矩阵 \mathbf Q
13. 目标分布 \tilde p
14. 输出： 采样出的一个状态序列 \{x\_0,x\_1,\cdots,x\_{n},x\_{n+1},\cdots\}
15. 算法：
16. 初始化 x\_0
17. 对 t=1,2,\cdots 执行迭代。迭代步骤如下：
18. 根据 Q(x^{\*} \mid x\_{t-1}) 采样出候选样本 x^{\*} ，其中 Q 为转移概率函数。
19. 计算 \alpha(x^{\*} \mid x\_{t-1}) ： \alpha( x^{\*} \mid x\_{t-1})=\min \left(1,\frac{\tilde p( x^{\*})Q( x \_{t-1} \mid x^{\*})}{\tilde p( x \_{t-1})Q( x^{\*} \mid x \_{t-1})} \right)
20. 根据均匀分布从 (0,1) 内采样出阈值 u，如果 u \le \alpha( x^{\*} \mid x\_{t-1}) ，则接受 x^{\*} ， 即 x\_{t}= x^{\*} 。否则拒绝 x^{\*} ， 即 x\_{t}= x\_{t-1} 。
21. 返回采样得到的序列 \{x\_0,x\_1,\cdots,x\_{n},x\_{n+1},\cdots\}

注意：返回的序列中，只有充分大的 n 之后的序列 \{ x\_{n}, x\_{n+1},\cdots\} 才是服从 \tilde p 分布的采样序列。

3.2 Gibbs 算法

1. MH算法不仅可以应用于一维空间的采样，也适合高维空间的采样。 对于高维的情况，由于接受率 \alpha 的存在（通常 \alpha \lt 1），MH算法的效率通常不够高，此时一般使用 Gibbs 采样算法。
2. 考虑二维的情形：假设有概率分布 \tilde p(x,y) ，考察状态空间上 x 坐标相同的两个点 A(x\_1,y\_1),B(x\_1,y\_2) ，可以证明有： \tilde p(x\_1,y\_1)\tilde p(y\_2\mid x\_1)=\tilde p(x\_1)\tilde p(y\_1\mid x\_1)\tilde p(y\_2\mid x\_1)\\ \tilde p(x\_1,y\_2)\tilde p(y\_1\mid x\_1)=\tilde p(x\_1)\tilde p(y\_2\mid x\_1)\tilde p(y\_1\mid x\_1) 于是 \tilde p(x\_1,y\_1)\tilde p(y\_2\mid x\_1)=\tilde p(x\_1,y\_2)\tilde p(y\_1\mid x\_1) 。则在 x=x\_1 这条平行于 y 轴的直线上，如果使用条件分布 \tilde p(y\mid x\_1) 作为直线上任意两点之间的转移概率，则这两点之间的转移满足细致平稳条件。 同理：考察 y 坐标相同的两个点 A(x\_1,y\_1),C(x\_2,y\_1) ， 有 \tilde p(x\_1,y\_1)\tilde p(x\_2\mid y\_1)=\tilde p(x\_2,y\_1)\tilde p(x\_1\mid y\_1) 。在 y=y\_1 这条平行于 x 轴的直线上，如果使用条件分布 \tilde p(x\mid y\_1) 作为直线上任意两点之间的转移概率，则这两点之间的转移满足细致平稳条件。 可以构造状态空间上任意两点之间的转移概率矩阵 \mathbf Q： 对于任意两点 A=(x\_A,y\_A),B=(x\_B,y\_B)， 令从 A 转移到 B 的概率为 Q(A\rightarrow B)：
3. 如果 x\_A=x\_B=x^{\*} ，则 Q(A\rightarrow B)=\tilde p(y\_B\mid x^{\*}) 。
4. 如果 y\_A=y\_B=y^{\*} ，则 Q(A\rightarrow B)=\tilde p(x\_B\mid y^{\*}) 。
5. 否则 Q(A\rightarrow B)=0 。

采用该转移矩阵 \mathbf Q ，可以证明：对于状态空间中任意两点 A,B，都满足细致平稳条件：

\tilde p(A)Q(A\rightarrow B)=\tilde p(B)Q(B\rightarrow A)

于是这个二维状态空间上的马尔可夫链将收敛到平稳分布 \tilde p(x,y) ，这就是吉布斯采样的原理。

1. Gibbs 算法：
2. 输入：目标分布 \tilde p(\mathbf{\vec x}) ，其中 \mathbf{\vec x}\in \mathbb R^n
3. 输出： 采样出的一个状态序列 \{\mathbf{\vec x}\_0,\mathbf{\vec x}\_1,\cdots\}
4. 算法步骤：
5. 初始化：随机初始化 \mathbf{\vec x}\_0,t=0 。
6. 执行迭代，迭代步骤如下：
7. 随机或者以一定次序遍历索引 1,2,\cdots,n 。遍历过程为（设当前索引为 i ）：
8. 将 x\_{1},\cdots,x\_{i-1},x\_{i+1},\cdots,x\_{n} 保持不变，计算条件概率： \tilde p(x\_{i}\mid x\_{1}=x\_{t,1},\cdots,x\_{i-1}=x\_{t,i-1},x\_{i+1}=x\_{t,i+1},\cdots,x\_{n}=x\_{t,n}) 。 该条件概率就是状态转移概率 Q(A\rightarrow B)
9. 根据条件概率 \tilde p(x\_{i}\mid x\_{1}=x\_{t,1},\cdots,x\_{i-1}=x\_{t,i-1},x\_{i+1}=x\_{t,i+1},\cdots,x\_{n}=x\_{t,n}) 对分量 x\_{i} 进行采样，假设采样结果为 x\_{t,i}^\* ，则获得新样本 \mathbf{\vec x}\_{t+1}=(x\_{t,1},\cdots,x\_{t,i-1},x\_{t,i}^\*,x\_{t,i+1},\cdots,x\_{t,n})^T 。
10. 令 t\leftarrow t+1，继续遍历。
11. 最终返回一个状态序列 \{\mathbf{\vec x}\_0,\mathbf{\vec x}\_1,\cdots\} 。
12. 吉布斯采样Gibbs sampling 有时被视作MH算法的特例，它也使用马尔可夫链获取样本。