前置知识

1. PDF：概率密度函数

2. CDF：累积概率密度函数

3. 示性函数(Indicator function)

   $$\begin{array}{}\text{(1)}& \begin{array}{c}{\mathbb{I}}_X(x)={\mathbb{I}}_{\{x\in X\}}=\{\begin{array}{c}1x\in X\\ 0x\notin X\end{array}\end{array}\end{array}$$

   示性函数有如下性质：示性函数的期望等于下标事件发生的概率

   $$\begin{array}{}\text{(2)}& \mathbf{E}{\mathbb{I}}_{\{x\in X\}}=1\cdot P\{x\in X\}+0\cdot P\{x\notin X\}=P\{x\in X\}\end{array}$$

4. 指数分布

   $$\begin{array}{}\text{(3)}& \begin{array}{rl}f(x)& =\theta \cdot e^{-\theta x}(x>0)\\ F(x)& =1-e^{-\theta x}(x>0)\end{array}\end{array}$$

5. 高斯分布

   $$\begin{array}{}\text{(4)}& P(x)=\varphi (x)=\frac{1}{\sqrt{2\pi }\cdot \sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}\end{array}$$

MCMC简介

蒙特卡罗法（Monte Carlo method）是通过从概率模型随机抽样然后进行近似数值计算的方法。马尔可夫链蒙特卡罗法（Markov Chain Monte Carlo, MCMC），则是以马尔可夫链为概率模型的蒙特卡罗法。MCMC构建一个马尔可夫链，使其平稳分布就是要进行抽样的分布，首先基于该马尔可夫链进行随机游走，产生样本的序列，之后使用该平稳分布的样本进行近似数值计算。

MCMC包含MH（Metropolis-Hastings）算法和吉布斯抽样算法。

MCMC可用于概率分布的估计、定积分的近似计算、最优化问题求解等问题。

蒙特卡罗法

随机抽样

蒙特卡罗法要解决的问题是：假设概率分布已知，通过抽样获得概率分布的随机样本，然后使用随机样本对概率分布的特征进行分析。例如，通过计算样本均值，从而估计总体期望。

因此，蒙特卡洛法的核心是随机抽样。

一般的蒙特卡罗法包括：直接抽样、接受-拒绝抽样、重要性抽样。后两者适合于概率密度函数复杂，不能直接抽样的情况。

直接抽样（逆采样）

对于均匀分布，一般可以使用线性同余发生器生成（0,1）之间的伪随机数。那么，对于其它分布，将其采样转换为（0,1）均匀分布进行采样。转换的方法为：

1. 确定待采样分布的概率密度函数：$f(x)$$f(x)$，待采样的随机变量为:$x$$x$

2. 假设（0,1）均匀分布的随机变量为：$z$$z$

3. 求该概率密度函数的累积概率密度函数：$F(x)$$F(x)$

4. 由于$0\le z\le 1$$0 \le z \le 1$，$0\le F(x)\le 1$$0 \le F(x) \le 1$，因此可以定义$z=F(x)$$z = F(x)$，解得：$x=G(z)$$x = G(z)$。

   所以，通过均匀分布采样得到$z$$z$，带入$G(z)$$G(z)$即可得到$x$$x$。

   下图以指数分布进行映射示例：

   

接受-拒绝采样

假设$p(x)$$p(x)$不可以直接抽样。可以找一个可以直接抽样的分布，称为建议分布。假设$q(x)$$q(x)$是建议分布的概率密度函数，并且$k\cdot q(x)>p(x),k>0$$k\cdot q(x) \gt p(x) \; ,k \gt 0$。定义接受-拒绝算法如下：

输入：待抽样概率密度函数$p(x)$$p(x)$，满足条件的$k\cdot q(x)$$k \cdot q(x)$，

输出：随机样本$x_1,x_2,\cdots ,x_n$$x_1, x_2, \cdots, x_n$，其中，$n$$n$为参数

步骤：

1. 从$q(x)$$q(x)$中抽样，得到样本$x^{\ast }$$x^*$

2. 从(0,1)均匀分布中抽样，得到$u$$u$

3. 如果$u\le \frac{p(x^{\ast })}{k\cdot q(x^{\ast })}$$\Large u \le \frac{p(x^*)}{k \cdot q(x^*) }$，则将$x^{\ast }$$x^*$作为抽样样本，否则，回到步骤1

4. 直到得到$n$$n$个随机样本，结束。

重要性采样

函数$f(x)$$f(x)$在分布$p(x)$$p(x)$下的期望也可以写为：

重要性采样（Importance Sampling）是通过引入重要性权重，将分布$p(x)$$p(x)$下$f(x)$$f(x)$的期望变为在分布$q(x)$$q(x)$下$f(x)w(x)$$f(x)w(x)$的期望，从而可以近似为：

案例：$f(x)\sim \mathbf{N}(0,1)$$f(x) \sim \mathbf{N}(0,1)$，使用采样方法求$P(X>8)$$P(X > 8)$的概率。

解答：如果直接从$f(x)$$f(x)$中采样，只有极少数的点会大于8，因此考虑使用重要性采样法求解。

定义$g(x)\sim \mathbf{N}(8,1)$$g(x) \sim \mathbf{N}(8,1)$，将目标问题进行如下转换：

定义$w(x)={\mathbb{I}}_{x>8}{\mathbf{e}}^{32-8x}$$w(x) = \mathbb{I}_{x > 8} \mathbf{e} ^{32-8x}$。此时，可以将原问题转为函数$w(x)$$w(x)$关于概率分布$g(x)$$g(x)$期望的问题。接下来，可以从$g(x)$$g(x)$中进行采样，得到样本序列$x_1,x_2,\cdots ,x_N$$x_1,x_2,\cdots,x_N$，将采样点带入$w(x)$$w(x)$，可以得到对概率的估计值：

应用

期望计算

假设有随机变量$x\in \mathcal{X}$$x \in \mathcal{X}$，其概率密度函数为$p(x)$$p(x)$，$f(x)$$f(x)$为定义在$\mathcal{X}$$\mathcal{X}$上的函数，目标是求函数$f(x)$$f(x)$关于密度函数$p(x)$$p(x)$的数学期望$E_{p(x)}[f(x)]$$E_{p(x)}[f(x)]$。针对该问题，蒙特卡罗法通过如下步骤进行解决：

1. 按照概率分布$p(x)$$p(x)$独立地抽取$n$$n$个样本$x_1,x_2,\cdots ,x_n$$x_1, x_2, \cdots, x_n$

2. 计算函数$f(x)$$f(x)$的样本均值${\hat{f}}_n$$\hat{f}_n$

   $$\begin{array}{}\text{(9)}& {\hat{f}}_n=\frac{1}{n}\sum _{i=1}^{n}f\left(x_i\right)\end{array}$$

3. 根据大数定律，当样本容量足够大时，样本均值以概率1收敛于数学期望：

   $$\begin{array}{}\text{(10)}& {\hat{f}}_n\to E_{p(x)}[f(x)],n\to \mathrm{\infty }\end{array}$$

4. 综上述，数学期望的近似值为：

积分计算

假设有一个函数$h(x)$$h(x)$，目标是计算该函数的积分：${\int }_{\mathcal{X}}h(x)\mathrm{d}x$$\int_{\mathcal{X}} h(x) \mathrm{d} x$

如果将$h(x)$$h(x)$分解为一个函数$f(x)$$f(x)$和一个概率密度函数$p(x)$$p(x)$乘积的形式，那么有：

因此，函数$h(x)$$h(x)$的积分可以表示为函数$f(x)$$f(x)$关于概率密度函数$p(x)$$p(x)$的期望。而数学期望又可以通过样本均值进行估算。于是，可以通过样本均值近似计算积分，即：

马尔可夫链

基本定义

给定一个随机变量序列$X=\{X_0,X_1,\cdots ,X_t,\cdots \}$$X = \{X_0, X_1, \cdots, X_t, \cdots \}$，$X_t$$X_t$表示时刻$t$$t$的随机变量，该随机变量在不同时刻的取值集合相同，即状态空间相同，表示为$\mathcal{S}$$\mathcal{S}$。随机变量可以是离散的，也可以是连续的。

如果随机变量的状态仅依赖于前一个状态，即：

该性质称为马尔可夫性。

具有马尔可夫性的随机序列$X=\{X_0,X_1,\cdots ,X_t,\cdots \}$$X = \{X_0, X_1, \cdots, X_t, \cdots \}$称为马尔可夫链，或马尔科夫过程。条件概率分布$P(X_t\mid X_{t-1})$$P(X_t \mid X_{t-1})$称为马尔可夫链的转移概率分布。

若转移概率分布$P(X_t\mid X_{t-1})$$P(X_t \mid X_{t-1})$与$t$$t$无关，即：

则称该马尔可夫链为时间齐次的马尔可夫链。本文的马尔可夫链都是时间齐次的。

离散状态马尔可夫链

转移概率矩阵和状态分布

若马尔可夫链在时刻$t-1$$t-1$处于状态$j$$j$，在时刻$t$$t$移动到状态$i$$i$，将转移概率记作：

马尔可夫链的转移概率$p_{ij}$$p_{ij}$可由矩阵表示，即：

考虑马尔可夫链$X=\{X_0,X_1,\cdots ,X_t,\cdots \}$$X = \{X_0, X_1, \cdots, X_t, \cdots \}$在时刻$t(t=0,1,2,\cdots )$$t \;(t=0,1,2,\cdots)$的概率分布，称为时刻$t$$t$的状态分布，记作：

对于初始分布$\pi (0)$$\pi(0)$，其向量通常只有一个分量为1，其余为0，表示马尔可夫链从一个具体状态开始。

马尔可夫链$X$$X$在时刻$t$$t$的状态分布，可以由时刻$t-1$$t-1$的状态分布以及转移概率分布决定：

这是因为（对时刻$t$$t$状态为$i$$i$进行推导，其余状态类似）：

假设有三个状态，那么转移公式为：

马尔可夫链在时刻$t$$t$的状态分布，可以通过递推得到。由上述公式：

通过递推得到：

其中，$P^t$$P^t$称为$t$$t$步转移概率矩阵，表示时刻$0$$0$从状态$j$$j$出发，时刻$t$$t$到达状态$i$$i$的$t$$t$步转移概率：

平稳分布

对马尔可夫链$X=\{X_0,X_1,\cdots ,X_t,\cdots \}$$X = \{X_0, X_1, \cdots, X_t, \cdots \}$，其状态空间为$\mathcal{S}$$\mathcal{S}$，转移概率矩阵为$P=(p_{ij})$$P=(p_{ij})$，如果存在状态空间$\mathcal{S}$$\mathcal{S}$上的一个分布：

使得：

则称$\pi$$\pi$为马尔可夫链$X=\{X_0,X_1,\cdots ,X_t,\cdots \}$$X = \{X_0, X_1, \cdots, X_t, \cdots \}$的平稳分布。

直观上，如果马尔可夫链的平稳分布存在，那么以该分布为初始分布，然后进行随机状态转移，之后任何一个时刻的状态分布都是平稳分布。

连续状态马尔可夫链

将离散状态马尔可夫链中的转移概率矩阵替换为转移核，其它概率类似，不再赘述。

马尔可夫链性质

1. 不可约

   从任意状态出发，当经过充分长时间后，可以到达任意状态。

2. 非周期

   不存在一个状态，从该状态出发，再返回该状态，经历的时间呈周期性。

3. 正常返

   对其中任意一个状态，从其它任意一个状态出发，当时间趋于无穷时，首次转移到这个状态的概率不为0。

4. 遍历定理

   如果马尔可夫链是不可约、非周期且正常返的，则该马尔可夫链有唯一平稳分布。

   满足相应条件的马尔可夫链，当时间趋于无穷时，其状态分布接近于平稳分布。（注：可以理解为首先初始化一个状态分布，然后通过无数次迭代转移，最后状态分布趋于平稳分布）。

5. 可逆马尔可夫链

   如果对马尔可夫链，满足：

   $$\begin{array}{}\text{(27)}& {\pi }_j\cdot P(X_t=i\mid X_{t-1}=j)={\pi }_i\cdot P(X_{t-1}=j\mid X_t=i),i,j=1,2,\cdots \end{array}$$

   或者简写为（此式称为细致平衡方程）：

   $$\begin{array}{}\text{(28)}& {\pi }_j\cdot p_{ij}={\pi }_i\cdot p_{ji}\end{array}$$

   则称该马尔可夫链为可逆的马尔可夫链。

6. （定理）满足细致平衡方程的状态分布$\pi$$\pi$，就是该马尔可夫链的平稳分布，即$P\pi =\pi$$P \pi = \pi$，证明如下：

   $$\begin{array}{}\text{(29)}& P{\pi }_i=\sum _j{p}_{ij}{\pi }_j=\sum _j{p}_{ji}{\pi }_i={\pi }_i\sum _j{p}_{ji}={\pi }_i,i=1,2,\cdots \end{array}$$

马尔可夫链蒙特卡罗法

基本步骤

马尔可夫链蒙特卡罗法的基本步骤如下：

1. 在随机变量$x$$x$的状态空间$\mathcal{S}$$\mathcal{S}$上构造一个满足遍历定理的马尔可夫链，使其平稳分布为目标分布$p(x)$$p(x)$

2. 从状态空间的某一点$x_0$$x_0$出发，用构造的马尔可夫链进行随机游走，产生样本序列$x_0,x_1,x_2,\cdots ,x_t,\cdots$$x_0, x_1, x_2,\cdots,x_t,\cdots$

3. 应用马尔可夫链的遍历定理，确定正整数$m$$m$和$n$$n$，$m<n$$m \lt n$，前$m$$m$个样本处于燃烧期，被丢弃掉。得到样本集合$x_{m+1},x_{m+2},\cdots ,x_n$$x_{m+1}, x_{m+2}, \cdots, x_n$，求得函数$f(x)$$f(x)$的均值（遍历均值）

   $$\begin{array}{}\text{(30)}& \hat{E}f=\frac{1}{n-m}\sum _{i=m+1}^{n}f\left(x_i\right)\end{array}$$

需要考虑的几个重要问题：

1. 如何定义马尔可夫链，保证马尔可夫链蒙特卡罗法的条件成立；
2. 如何确定收敛步数$m$$m$，保证样本抽取的无偏性；
3. 如何确定迭代步数$n$$n$，保证遍历均值计算的精度。

Metropolis-Hastings 算法

马尔可夫链

假设要抽样的概率分布为$p(x)$$p(x)$，M-H算法采用转移核为$p(x,x^{\prime })$$p(x, x^{\prime})$的马尔可夫链：

建议分布$q(x,x^{\prime })$$q\left(x, x^{\prime}\right)$为另一个马尔可夫链的转移核，并且是不可约的，同时是一个容易抽样的分布。

接受分布$\alpha (x,x^{\prime })$$\alpha\left(x, x^{\prime}\right)$定义如下：

那么，转移核$p(x,x^{\prime })$$p\left(x, x^{\prime}\right)$可以写为：

转移核为$p(x,x^{\prime })$$p\left(x, x^{\prime}\right)$的马尔可夫链上的随机游走以如下方式进行：

1. 假设时刻$t-1$$t-1$状态为$x$$x$，即$x_{t-1}=x$$x_{t-1}=x$

2. 按建议分布$q(x,x^{\prime })$$q\left(x, x^{\prime}\right)$抽样产生一个候选状态$x^{\prime }$$x^{\prime}$

3. 按照接受分布$\alpha (x,x^{\prime })$$\alpha\left(x, x^{\prime}\right)$抽样决定是否接受$x^{\prime }$$x^{\prime}$：接受概率为$\alpha (x,x^{\prime })$$\alpha\left(x, x^{\prime}\right)$，拒绝概率为$1-\alpha (x,x^{\prime })$$1-\alpha\left(x, x^{\prime}\right)$。具体地，从区间(0,1)上均匀分布抽取一个随机数$u$$u$，决定时刻$t$$t$的状态：

   $$\begin{array}{}\text{(34)}& \begin{array}{c}{x}_t=\{\begin{array}{cc}{x}^{\prime },& u⩽\alpha (x,x^{\prime })\\ x,& u>\alpha (x,x^{\prime })\end{array}\end{array}\end{array}$$

   如果接受，则时刻$t$$t$转移到$x^{\prime }$$x^{\prime}$，否则，仍停留在状态$x$$x$

可以证明，转移核为$p(x,x^{\prime })$$p\left(x, x^{\prime}\right)$的马尔可夫链是可逆马尔可夫链（满足遍历定理），其平稳分布就是$p(x)$$p(x)$，即要抽样的目标分布。

证明过程：

1. 如果$x=x^{\prime }$$x=x^{\prime}$，显然$p(x)p(x,x^{\prime })=p\left(x^{\prime }\right)p(x^{\prime },x)$$p(x) p\left(x, x^{\prime}\right)=p\left(x^{\prime}\right) p\left(x^{\prime}, x\right)$成立

2. 如果$x\ne x^{\prime }$$x \neq x^{\prime}$，则

   $$\begin{array}{}\text{(35)}& \begin{array}{rl}p(x)p(x,x^{\prime })& =p(x)q(x,x^{\prime })min\{1,\frac{p\left(x^{\prime }\right)q(x^{\prime },x)}{p(x)q(x,x^{\prime })}\}\\ & =min\{p(x)q(x,x^{\prime }),p\left(x^{\prime }\right)q(x^{\prime },x)\}\\ & =p\left(x^{\prime }\right)q(x^{\prime },x)min\{\frac{p(x)q(x,x^{\prime })}{p\left(x^{\prime }\right)q(x^{\prime },x)},1\}\\ & =p\left(x^{\prime }\right)q(x^{\prime },x)\alpha (x^{\prime },x)\\ & =p\left(x^{\prime }\right)p(x^{\prime },x)\end{array}\end{array}$$

   因此，$p(x)p(x,x^{\prime })=p\left(x^{\prime }\right)p(x^{\prime },x)$$p(x) p\left(x, x^{\prime}\right)=p\left(x^{\prime}\right) p\left(x^{\prime}, x\right)$也成立，该马尔可夫链为可逆马尔可夫链。

3. 由$p(x)p(x,x^{\prime })=p\left(x^{\prime }\right)p(x^{\prime },x)$$p(x) p\left(x, x^{\prime}\right)=p\left(x^{\prime}\right) p\left(x^{\prime}, x\right)$可知：

   $$\begin{array}{}\text{(36)}& \begin{array}{rl}\int p(x)p(x,x^{\prime })\mathrm{d}x& =\int p\left(x^{\prime }\right)p(x^{\prime },x)\mathrm{d}x\\ & =p\left(x^{\prime }\right)\int p(x^{\prime },x)\mathrm{d}x\\ & =p\left(x^{\prime }\right)\end{array}\end{array}$$

   根据平稳分布的定义，$p(x)$$p(x)$是马尔可夫链的平稳分布。

建议分布

建议分布$q(x,x^{\prime })$$q\left(x, x^{\prime}\right)$有多种可能的形式，以下介绍常用的两种。

1. 对称式

   假设建议分布是对称的，即对任意的$x$$x$和$x^{\prime }$$x^{\prime}$，有：

   $$\begin{array}{}\text{(37)}& q(x,x^{\prime })=q(x^{\prime },x)\end{array}$$

   如果选择该分布，也可称为Metropolis算法。此时，接受分布$\alpha (x,x^{\prime })$$\alpha\left(x, x^{\prime}\right)$可以简化为：

   $$\begin{array}{}\text{(38)}& \alpha (x,x^{\prime })=min\{1,\frac{p\left(x^{\prime }\right)}{p(x)}\}\end{array}$$

   Metropolis选择的一个特例是$q(x,x^{\prime })$$q(x,x^{\prime})$取条件概率分布$p(x^{\prime }\mid x)$$p(x^{\prime} \mid x)$，定义为多元正态分布，其均值为$x$$x$，其协方差矩阵是常数矩阵。

   Metropolis选择的特点是当$x^{\prime }$$x^{\prime}$与$x$$x$接近时，$q(x,x^{\prime })$$q(x,x^{\prime})$的概率值高，否则$q(x,x^{\prime })$$q(x,x^{\prime})$的概率值低。状态转移在附近点的可能性更大。

    

2. 独立抽样

   假设$q(x,x^{\prime })$$q(x,x^{\prime})$与当前状态$x$$x$无关，即$q(x,x^{\prime })=q(x^{\prime })$$q(x,x^{\prime}) = q(x^{\prime})$，建议分布的计算按照$q(x^{\prime })$$q(x^{\prime})$独立抽样进行。

   独立抽样实现简单，但是收敛速度可能较慢，通常选择接近目标分布$p(x)$$p(x)$的分布作为建议分布$q(x)$$q(x)$。

满条件分布

马尔可夫链蒙特卡罗法的目标分布通常是多元联合概率分布$p(x)=p(x_1,x_2,\cdots ,x_k)，\text{其}\text{中}x=(x_1,x_2,\cdots ,x_k{)}^T$$p(x)=p(x_1, x_2, \cdots,x_k)，其中x=(x_1, x_2, \cdots,x_k)^T$为$k$$k$维随机变量。如果条件概率分布$p(x_I\mid x_{-I})$$p(x_I \mid x_{-I})$中所有$k$$k$个变量都出现，其中$x_I=\{x_i,i\in I\},x_{-I}=\{x_i,i\notin I\},I\subset K=\{1,2,\cdots ,k\}$$x_I=\left\{x_i, i \in I\right\}, x_{-I}=\left\{x_i, i \notin I\right\}, I \subset K=\{1,2, \cdots, k\}$，那么这种条件概率分布为满条件概率分布。

满条件概率分布有以下性质，对任意的$x,x^{\prime }\in \mathcal{X}$$x,x^{\prime} \in \mathcal{X}$和任意的$I\subset K$$I \subset K$，有：

在MH算法中，可以利用该性质简化计算，提高计算效率。具体地，$\frac{p(x_I^{\prime }\mid x_{-I}^{\prime })}{p(x_I\mid x_{-I})}$$\frac{p\left(x_I^{\prime} \mid x_{-I}^{\prime}\right)}{p\left(x_I \mid x_{-I}\right)}$比$\frac{p\left(x^{\prime }\right)}{p(x)}$$\frac{p\left(x^{\prime}\right)}{p(x)}$更容易计算。

MH算法步骤

输入：待抽样的目标分布$p(x)$$p(x)$，函数$f(x)$$f(x)$，建议分布$q(x)$$q(x)$

输出：$p(x)$$p(x)$的随机样本$x_{m+1},x_{m+2},\cdots ,x_n$$x_{m+1}, x_{m+2}, \cdots, x_n$，函数样本均值$f_{mn}$$f_{mn}$

参数：收敛步数$m$$m$，迭代步数$n$$n$

步骤：

1. 任意选择一个初始值$x_0$$x_0$

2. 对$i=1,2,\cdots ,n$$i=1,2,\cdots,n$循环执行

   a) 设状态$x_{i-1}=x$$x_{i-1} = x$，按照建议分布$q(x,x^{\prime })$$q(x, x^{\prime})$随机抽取一个候选状态$x^{\prime }$$x^{\prime}$

   b) 计算接受概率

   $$\begin{array}{}\text{(40)}& \alpha (x,x^{\prime })=min\{1,\frac{p\left(x^{\prime }\right)q(x^{\prime },x)}{p(x)q(x,x^{\prime })}\}\end{array}$$

   c) 从区间(0,1)中按均匀分布随机抽取一个数$u$$u$

   如$u\le \alpha (x,x^{\prime })$$u \le \alpha\left(x, x^{\prime}\right)$，则状态$x_i=x^{\prime }$$x_i = x^{\prime}$，否则，状态$x_i=x$$x_i = x$

3. 得到样本集合$x_{m+1},x_{m+2},\cdots ,x_n$${x_{m+1}, x_{m+2}, \cdots, x_n}$，计算

   $$\begin{array}{}\text{(41)}& f_{mn}=\frac{1}{n-m}\sum _{i=m+1}^{n}f\left(x_i\right)\end{array}$$

吉布斯抽样

基本原理

吉布斯抽样是MH算法的特殊情况，多用于多元变量联合分布的抽样和估计。其基本做法是，从联合概率分布定义满条件概率分布，依次对满条件概率分布进行抽样，得到样本的序列。

假设多元变量的联合概率分布为$p(x)=p(x_1,x_2,\cdots ,x_k)$$p(x) = p(x_1, x_2, \cdots, x_k)$。吉布斯抽样从一个初始样本$x^{(0)}=(x_1^{(0)},x_2^{(0)},\cdots ,x_k^{(0)}{)}^{\mathbf{T}}$$x^{(0)} = (x_1^{(0)},x_2^{(0)},\cdots,x_k^{(0)})^\mathbf{T}$出发，不断迭代，每次迭代得到联合分布的一个样本$x^{(i)}=(x_1^{(i)},x_2^{(i)},\cdots ,x_k^{(i)}{)}^{\mathbf{T}}$$x^{(i)} = (x_1^{(i)},x_2^{(i)},\cdots,x_k^{(i)})^\mathbf{T}$。最终，得到抽样样本序列$x^{(0)},x^{(1)},\cdots ,x^{(n)}$$x^{(0)},x^{(1)},\cdots,x^{(n)}$。

在每次迭代中，依次对$k$$k$个随机变量中的一个变量进行抽样。例如，在第$i$$i$次迭代中，对第$j$$j$个变量进行随机抽样，那么抽样的分布是满条件概率分布$p(x_j\mid x_{-j}^{(i)})$$p(x_j \mid x_{-j}^{(i)})$。

假设在第$i-1$$i-1$步得到样本$(x_1^{(i-1)},x_2^{(i-1)},\cdots ,x_k^{(i-1)}{)}^{\mathbf{T}}$$(x_1^{(i-1)},x_2^{(i-1)},\cdots,x_k^{(i-1)})^\mathbf{T}$，在第$i$$i$步，首先对第一个变量按照以下满条件概率分布随机抽样：

于是得到第$i$$i$步的整体样本$x^{(i)}=(x_1^{(i)},x_2^{(i)},\cdots ,x_k^{(i)}{)}^{\mathbf{T}}$$x^{(i)} = (x_1^{(i)},x_2^{(i)},\cdots,x_k^{(i)})^\mathbf{T}$。

吉布斯抽样是单分量MH算法的特殊情况。定义建议分布是当前变量$x_j(j=1,2,\cdots ,k)$$x_j \; (j=1,2,\cdots, k)$的满条件概率分布：

这时，接受概率$\alpha =1$$\alpha = 1$：

这里用到$p\left(x_{-j}\right)=p\left(x^{\prime }{}_{-j}\right)$$p\left(x_{-j}\right)=p\left(x^{\prime}{ }_{-j}\right)$和$p(\cdot \mid x_{-j})=p(\cdot \mid x^{\prime }{}_{-j})$$p\left(\cdot \mid x_{-j}\right)=p\left(\cdot \mid x^{\prime}{ }_{-j}\right)$

转移核就是满条件概率分布：

吉布斯抽样对每次抽样的结果都接受，没有拒绝，这一点和一般的MH算法不同。

这里，假设满条件概率分布$p(x_j^{\prime }\mid x_{-j})$$p\left(x_j^{\prime} \mid x_{-j}\right)$不为0，即马尔可夫链是不可约的。

吉布斯抽样算法步骤

输入：目标概率分布的密度函数$p(x)$$p(x)$，函数$f(x)$$f(x)$

输出：$p(x)$$p(x)$的随机样本$x_{m+1},x_{m+2},\cdots ,x_n$$x_{m+1}, x_{m+2}, \cdots, x_n$，函数样本均值$f_{mn}$$f_{mn}$

参数：收敛步数m，迭代步数n

步骤：

(1) 初始化样本$x^{(0)}=(x_1^{(0)},x_2^{(0)},\cdots ,x_k^{(0)}{)}^{\mathbf{T}}$$x^{(0)} = (x_1^{(0)},x_2^{(0)},\cdots,x_k^{(0)})^\mathbf{T}$

(2) 对步数$i$$i$进行循环

设第$i-1$$i-1$次迭代的样本为$x^{(i-1)}=(x_1^{(i-1)},x_2^{(i-1)},\cdots ,x_k^{(i-1)}{)}^{\mathbf{T}}$$x^{(i-1)} = (x_1^{(i-1)},x_2^{(i-1)},\cdots,x_k^{(i-1)})^\mathbf{T}$，则第$i$$i$次迭代进行如下几步操作：

(1) 由满条件概率分布$p(x_1\mid x_2^{(i-1)},\cdots ,x_k^{(i-1)})$$p\left(x_1 \mid x_2^{(i-1)}, \cdots, x_k^{(i-1)}\right)$抽取$x_1^{(i)}$$x_1^{(i)}$

$⋮$$\vdots$

(j) 由满条件概率分布$p(x_j\mid x_1^{(i)},\cdots ,x_{j-1}^{(i)},x_{j+1}^{(i-1)},\cdots ,x_k^{(i-1)})$$p\left(x_j \mid x_1^{(i)}, \cdots, x_{j-1}^{(i)}, x_{j+1}^{(i-1)}, \cdots, x_k^{(i-1)}\right)$抽取$x_j^{(i)}$$x_j^{(i)}$

$⋮$$\vdots$

(k) 由满条件概率分布$p(x_k\mid x_1^{(i)},\cdots ,x_{k-1}^{(i)})$$p\left(x_k \mid x_1^{(i)}, \cdots, x_{k-1}^{(i)}\right)$抽取$x_k^{(i)}$$x_k^{(i)}$

得到第$i$$i$次迭代值$x^{(i)}={(x_1^{(i)},x_2^{(i)},\cdots ,x_k^{(i)})}^{\mathrm{T}}$$x^{(i)}=\left(x_1^{(i)}, x_2^{(i)}, \cdots, x_k^{(i)}\right)^{\mathrm{T}}$

(3) 得到样本集合

(4) 计算

吉布斯抽样示例

假设要采样的是二维正态分布$N(\mu ,\mathrm{\Sigma })$$N(\mu, \Sigma)$，其中：

首先，需要求状态转移满条件分布：

证明：已知

因此：

确定了目标分布的满条件分布，接下来可以使用吉布斯抽样得到样本：

1. 采样一个初始状态：$(x_1^{(0)},x_2^{(0)})$$(x_1^{(0)}, x_2^{(0)})$，设燃烧期样本数为$m$$m$，需要得到的样本数为$n$$n$

2. for i in range(1, m+n):

   $x_1^{(i)}\sim P(x_1\mid x_2^{(i-1)})$$x_1^{(i)} \sim P\left(x_1 \mid x_2^{(i-1)} \right)$

   $x_2^{(i)}\sim P(x_2\mid x_1^{(i-1)})$$x_2^{(i)} \sim P\left(x_2 \mid x_1^{(i-1)} \right)$

3. 得到样本集合$\{(x_1^{(m+1)},x_2^{(m+1)}),(x_1^{(m+2)},x_2^{(m+2)}),\dots ,(x_1^{(m+n)},x_2^{(m+n)})\}$$\left\{\left(x_1^{(m+1)}, x_2^{(m+1)}\right),\left(x_1^{(m+2)}, x_2^{(m+2)}\right), \ldots,\left(x_1^{(m+n)}, x_2^{(m+n)}\right)\right\}$

参考文档

1. 李航 《统计学习方法》
2. 第十九课 马尔可夫链蒙特卡罗法 (julyedu.com)
3. 蒙特卡罗算法 - 知乎 (zhihu.com)
4. (1 封私信 / 62 条消息) 如何理解蒙特卡洛方法的接受拒绝采样？ - 知乎 (zhihu.com)
5. 使用蒙特卡洛计算定积分（附Python代码） - Ji-Huan Guan (guanjihuan.com)
6. (145条消息) MCMC方法与示例禾雨辰子的博客-CSDN博客_mcmc方法
7. 你一定从未看过如此通俗易懂的马尔科夫链蒙特卡罗方法(MCMC)解读(上) - 知乎 (zhihu.com)
8. 你一定从未看过如此通俗易懂的马尔科夫链蒙特卡罗方法(MCMC)解读(下) - 知乎 (zhihu.com)
9. (144条消息) python实现Metropolis采样算法实例_David-Chow的博客-CSDN博客
10. Monte Carlo笔记-2：MCMC，附简易python代码 - simplex - 博客园 (cnblogs.com)
11. 用Python实现马尔可夫链蒙特卡罗 - 知乎 (zhihu.com)
12. 如何快速理解马尔科夫链蒙特卡洛法？ - 知乎 (zhihu.com)
13. 机器学习中的数学-07讲-常见的采样方法4 (重要性采样1)哔哩哔哩bilibili
14. Markov Chain Monte Carlo (MCMC) — Computational Statistics in Python 0.1 documentation (duke.edu)