目录

• 蒙特卡罗法是通过从概率模型的随机抽样进行近似数值计算的方法。马尔可夫链蒙特卡罗法则是以马尔可夫链为概率模型的蒙特卡罗法。马尔可夫链蒙特卡罗法构建一个马尔可夫链，使其平稳分布就是要进行抽样的分布，首先基于该马尔可夫链进行随机游走，产生样本的序列，之后使用该平稳分布的样本进行近似数值计算
• 马尔可夫链蒙特卡罗法被应用于概率分布的估计、定积分的近似计算、最优化问题的近似求解等问题，特别是被应用于统计学习中概率模型的学习与推理，是重要的统计学习计算方法

蒙特卡罗法 (Monte Carlo Method)

随机抽样 (random sampling)

• 蒙特卡罗法要解决的问题是，假设概率分布的定义己知，通过抽样获得概率分布的随机样本，并通过得到的随机样本对概率分布的特征进行分析。所以蒙特卡罗法的核心是随机抽样
  
  • e.g. 从样本计算出样本均值，从而估计总体期望
• 一般的蒙特卡罗法有直接抽样法、接受-拒绝抽样法、重要性抽样法等。接受-拒绝抽样法、重要性抽样法适合于概率密度函数复杂 (如密度函数含有多个变量，各变量相互不独立，密度函数形式复杂)，不能直接抽样的情况。下面介绍接受-拒绝抽样法

马尔科夫链蒙特卡罗法也适合于概率密度函数复杂，不能直接抽样的情况，旨在解决一般的蒙特卡罗法，如接受-拒绝抽样法、重要性抽样法，抽样效率不高的问题。一般的蒙特卡罗法中的抽样样本是独立的，而马尔可夫链蒙特卡罗法中的抽样样本不是独立的，样本序列形成马尔科夫链

---

接受-拒绝抽样法 (accept-reject sampling method)

• 假设有随机变量 $x$，取值 $x\in \mathcal{X}$，其概率密度函数为 $p(x)$。目标是得到该概率分布的随机样本，以对这个概率分布进行分析
• 接受-拒绝法的基本想法如下。假设 $p(x)$ 不可以直接抽样。找一个可以直接抽样的分布，称为建议分布 (proposal distribution)。假设 $q(x)$ 是建议分布的概率密度函数，并且有 $cq(x)\ge p(x)$，其中 $c>0$。按照 $q(x)$ 进行抽样，假设得到结果是 $x^{\ast }$，再按 $\frac{p(x^{\ast })}{cq(x^{\ast })}$ 的比例随机决定是否接受 $x^{\ast }$。直观上，落到 $p(x^{\ast })$ 范围内的就接受，落到 $p(x^{\ast })$ 范围外的就拒绝(。接受-拒绝法实际是按照 $p(x)$ 的涵盖面积 (或涵盖体积) 占 $cq(x)$ 的涵盖面积 (或涵盖体积) 的比例进行抽样
  
  
  • 接受-拒绝法的优点是容易实现，缺点是效率可能不高。如果 $p(x)$ 的涵盖体积占 $cq(x)$ 的涵盖体积的比例很低，就会导致拒绝的比例很高，抽样效率很低。注意，一般是在高维空间进行抽样，即使 $p(x)$ 与 $cq(x)$ 很接近，两者涵盖体积的差异也可能很大

数学期望估计 (estimation of mathematical expectation)

• 蒙特卡罗法可用于数学期望估计 (直接抽样法、接受-拒绝抽样法、重要性抽样法…)。假设有随机变量 $x$，目标是求函数 $f(x)$ 关于密度函数 $p(x)$ 的数学期望 $E_{p(x)}[f(x)]$。针对这个问题，蒙特卡罗法按照概率分布 $p(x)$ 独立地抽取 $n$ 个样本后计算函数 $f(x)$ 的样本均值 ${\hat{f}}_n$ 作为数学期望 $E_{p(x)}[f(x)]$ 的近似值
  

定积分的近似计算 (Monte Carlo integration)

• 假设有一个函数 $h(x)$， 目标是计算该函数的积分
  如果能够将函数 $h(x)$ 分解成一个函数 $f(x)$ 和一个概率密度函数 $p(x)$ 的乘积的形式，那么就有
  于是，就可以利用样本均值来近似计算积分。实际上，给定一个概率密度函数 $p(x)$，只要取 $f(x)=\frac{h(x)}{p(x)}$，就可得上式

---

马尔可夫链的性质

• 基础知识见 马尔可夫链 (Markov Chains)

以下介绍离散状态马尔可夫链的性质。可以自然推广到连续状态马尔可夫链

不可约

• 直观上，一个不可约的马尔可夫链，从任意状态出发，当经过充分长时间后，可以到达任意状态

---

例

• 转移概率矩阵
  平稳分布 $\pi =(0,0,1{)}^T$
• 此马尔可夫链，转移到状态 3 后，就在该状态上循环跳转，不能到达状态 1 和状态 2，因此该马尔可夫链是可约的

非周期

• 直观上， 一个非周期性的马尔可夫链，不存在一个状态， 从这一个状态出发， 再返回到这个状态时所经历的时间长呈一定的周期性

---

例

• 转移概率矩阵
  其平稳分布是 $\pi =(1/3,1/3,1/3{)}^T$
• 此马尔可夫链从每个状态出发，返回该状态的时刻都是 3 的倍数，具有周期性，因此该马尔可夫链是周期的

---

定理 19.2

• 不可约且非周期的有限状态马尔可夫链，有唯一平稳分布存在

正常返

• 直观上， 一个正常返的马尔可夫链，其中任意一个状态，从其他任意一个状态出发， 当时间趋于无穷时，首次转移到这个状态的概率不为 0

---

例

• 转移概率矩阵
  
• 当 $p>q$ 时，平稳分布是
  当时间趋于无穷时，转移到任何一个状态的概率不为 0，马尔可夫链是正常返的
• 当 $p\le q$ 时，不存在平稳分布，马尔可夫链不是正常返的

---

定理 19.3

• 不可约、非周期且正常返的马尔可夫链，有唯一平稳分布存在

遍历定理

• 遍历定理的直观解释: 满足相应条件的马尔可夫链，当时间趋于无穷时，马尔可夫链的状态分布趋近于平稳分布，随机变量的函数的样本均值以概率 1 收敛于该函数的数学期望
• 样本均值可以认为是时间均值，而数学期望是空间均值。遍历定理实际表述了遍历性的含义: 当时间趋于无穷时，时间均值等于空间均值。遍历定理的三个条件: 不可约、非周期、正常返，保证了当时间趋于无穷时达到任意一个状态的概率不为 0

---

• 理论上并不知道经过多少次迭代，马尔可夫链的状态分布才能接近于平稳分布，在实际应用遍历定理时，取一个足够大的整数 $m$，经过 $m$ 次迭代之后认为状态分布就是平稳分布，这时计算从第 $m+1$ 次迭代到第 $n$ 次迭代的均值，即
  称为遍历均值; 到时刻 $m$ 为止的时间段称为燃烧期 (burn-in time)

可逆马尔可夫链

• 直观上，如果有可逆的马尔可夫链，那么以该马尔可夫链的平稳分布作为初始分布，进行随机状态转移，无论是面向未来还是面向过去，任何一个时刻的状态分布都是该平稳分布

---

• 定理 19.5 说明，可逆马尔可夫链一定有唯一平稳分布，给出了一个马尔可夫链有平稳分布的充分条件 (不是必要条件)。也就是说，可逆马尔可夫链满足遍历定理 19.4 的条件

证明

马尔可夫链蒙特卡罗法

• 马尔可夫链蒙特卡罗法也是一种随机采样方法。相比传统的蒙特卡罗法，如接受-拒绝法、重要性抽样法，马尔可夫链蒙特卡罗法更适合于随机变量是多元的、密度函数是非标准形式的、随机变量各分量不独立等情况，且抽样效率更高

要解决的问题

• 假设多元随机变量 $x$，满足 $x\in \mathcal{X}$，其概率密度函数为 $p(x)$，$f(x)$ 为定义在 $x\in \mathcal{X}$ 上的函数，目标是获得概率分布 $p(x)$ 的样本集合，以及求函数 $f(x)$ 的数学期望 $E_{p(x)}[f(x)]$

基本想法

• 马尔可夫链蒙特卡罗法的基本想法是: 在随机变量 $x$ 的状态空间 $\mathcal{S}$ 上定义一个满足遍历定理的马尔可夫链
  
  
  
  
  X
  
  
  =
  
  
  {
  
  
  
  X
  
  
  0
  
  
  
  ,
  
  
  .
  
  
  .
  
  
  .
  
  
  ,
  
  
  
  X
  
  
  t
  
  
  
  ,
  
  
  .
  
  
  .
  
  
  .
  
  
  }
  
  
  
  X ={X_0,…, X_t,…}
  
  
  X={X0​,...,Xt​,...}，使其平稳分布就是抽样的目标分布 $p(x)$。然后在这个马尔可夫链上进行随机游走，每个时刻得到一个样本。根据遍历定理，当时间趋于无穷时，样本的分布趋近平稳分布，样本的函数均值趋近函数的数学期望。所以，当时间足够长时 (时刻大于某个正整数 $m$)，在之后的时间 (时刻小于等于某个正整数 $n$，$n>m$) 里随机游走得到的样本集合
  
  
  
  
  {
  
  
  
  x
  
  
  
  m
  
  
  +
  
  
  1
  
  
  
  
  ,
  
  
  
  x
  
  
  
  m
  
  
  +
  
  
  2
  
  
  
  
  ,
  
  
  .
  
  
  .
  
  
  .
  
  
  ,
  
  
  
  x
  
  
  n
  
  
  
  }
  
  
  
  {x_{m +1} , x_{m+2} ,…,x_n}
  
  
  {xm+1​,xm+2​,...,xn​} 就是目标概率分布的抽样结果，得到的函数均值 (遍历均值) 就是要计算的数学期望值:
  

可见，若随机变量 $x$ 是离散的，则应定义相应的离散状态马尔可夫链，状态数即为 $|\mathcal{X}|$；若随机变量 $x$ 是连续的，则应定义相应的连续状态马尔可夫链

• 马尔可夫链蒙特卡罗法中得到的样本序列，相邻的样本点是相关的，而不是独立的。因此，在需要独立样本时，可以在该样本序列中再次进行随机抽样，比如每隔一段时间取一次样本，将这样得到的子样本集合作为独立样本集合

---

难点

• 如何构建具体的马尔可夫链成为这个方法的关键。连续变量的时候，需要定义转移核函数。离散变量的时候，需要定义转移矩阵。一个方法是定义特殊的转移核函数或者转移矩阵，构建可逆马尔可夫链，这样可以保证遍历定理成立。常用的马尔可夫链蒙特卡罗法有 Metropolis-Hastings 算法、吉布斯抽样
• 如何确定收敛步数 $m$，保证样本抽样的无偏性。马尔可夫链蒙特卡罗法的收敛性的判断通常是经验性的，比如，在马尔可夫链上进行随机游走，检验遍历均值是否收敛。具体地，每隔一段时间取一次样本，得到多个样本以后，计算遍历均值，当计算的均值稳定后，认为马尔可夫链已经收敛。再比如，在马尔可夫链上并行进行多个随机游走，比较各个随机游走的遍历均值是否接近一致
• 如何确定迭代步数 $n$，保证遍历均值计算的精度

基本步骤

• (1) 首先，在随机变量 $x$ 的状态空间 $\mathcal{S}$ 上构造一个满足遍历定理的马尔可夫链，使其平稳分布为目标分布 $p(x)$
• (2) 从状态空间的某一点 $x_0$ 出发，用构造的马尔可夫链进行随机游走，产生样本序列 $x_0,...,x_t,..$
• (3) 应用马尔可夫链的遍历定理，确定正整数 $m$ 和 $n$，($m<n$)，得到样本集合
  
  
  
  
  {
  
  
  
  x
  
  
  
  m
  
  
  +
  
  
  1
  
  
  
  
  ,
  
  
  
  x
  
  
  
  m
  
  
  +
  
  
  2
  
  
  
  
  ,
  
  
  …
  
  
  ,
  
  
  
  x
  
  
  n
  
  
  
  }
  
  
  
  {x_{m +1}, x_{m+2}, … ,x_n }
  
  
  {xm+1​,xm+2​,…,xn​}，求得函数 $f(x)$ 的均值 (遍历均值〉
  

马尔可夫链蒙特卡罗法与统计学习 (贝叶斯学习)

• 马尔可夫链蒙特卡罗法在统计学习，特别是贝叶斯学习中，起着重要的作用。主要是因为马尔可夫链蒙特卡罗法可以用在概率模型的学习和推理上

---

• 假设观测数据由随机变量 $y\in \mathcal{Y}$ 表示，模型由随机变量 $x\in \mathcal{X}$ 表示，贝叶斯学习通过贝叶斯定理计算给定数据条件下模型的后验概率，并选择后验概率最大的模型。后验概率
  
• 贝叶斯学习中经常需要进行三种积分运算: 归范化 (normalization) 、边缘化 (marginalization) 、数学期望 (expectation)。当观测数据和模型都很复杂的时候，以上的积分计算变得困难。马尔可夫链蒙特卡罗法为这些计算提供了一个通用的有效解决方案
  
  • 后验概率计算中需要归范化计算:
    
  • 如果有隐变量 $z\in \mathcal{Z}$，后验概率的计算需要边缘化计算:
    
  • 如果有一个函数 $f(x)$，可以计算该函数的关于后验概率分布的数学期望:
    

Metropolis-Hastings 算法

• Metropolis-Hastings 算法是最基本的马尔可夫链蒙特卡罗法， Metropolis 等人在 1953 年提出原始的算法，Hastings 在1970 年对之加以推广，形成了现在的形式

马尔可夫链的构造

• 假设要抽样的概率分布为 $p(x)$。Metropolis-Hastings 算法采用转移核为 $p(x,x^{\prime })$ 的马尔可夫链，该马尔可夫链构成了一个可逆马尔可夫链，它的唯一平稳分布是 $p(x)$:
  其中 $q(x,x^{\prime })$ 和 $\alpha (x,x^{\prime })$ 分别称为建议分布 (proposal distribution) 和 接受分布 (acceptance distribution)

---

建议分布 $q(x,x^{\prime })$

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

---

接受分布 $\alpha (x,x^{\prime })$

---

转移核 $p(x,x^{\prime })$

• 转移核 $p(x,x^{\prime })$ 可以写成
  

在计算过程中，$q(x,x^{\prime })$ 和 $p(x)$ 并不需要归一化，只需要计算出相对大小即可，这对那些难以归一化的待采样概率分布而言十分友好，例如在贝叶斯学习中，后验概率的归一化因子就是一个难以计算的积分项

• 转移核为 $p(x,x^{\prime })$ 的马尔可夫链上的随机游走以以下方式进行。如果在时刻 $(t-1)$ 处于状态 $x$，即 $x_{t-1}=x$，则先按建议分布 $q(x,x^{\prime })$ 抽样产生一个候选状态 $x^{\prime }$. 然后按照接受分布 $\alpha (x,x^{\prime })$ 抽样决定是否接受状态 $x^{\prime }$。以概率 $\alpha (x,x^{\prime })$ 接受 $x^{\prime }$，决定时刻 $t$ 转移到状态 $x^{\prime }$. 而以概率 $1-\alpha (x,x^{\prime })$ 拒绝 $x^{\prime }$. 决定时刻 $t$ 仍停留在状态 $x$。具体地，从区间 $(0,1)$ 上的均匀分布中抽取一个随机数 $u$，决定时刻 $t$ 的状态：
  

---

下面证明转移核为 $p(x,x^{\prime })$ 的马尔可夫链是可逆马尔可夫链 (满足遍历定理)，其平稳分布就是 $p(x)$

证明

• 若 $x=x^{\prime }$. 则式 (19.41) 显然成立
• 若 $x\ne x^{\prime }$. 则
  式 (19.41) 成立
• 由式 (19.41) 知，
  根据平稳分布的定义，$p(x)$ 是马尔可夫链的平稳分布

建议分布的选择

• 建议分布 $q(x,x^{\prime })$ 有多种可能的形式，这里介绍两种常用形式

Metropolis 选择

• Metropolis 选择: 假设建议分布是对称的，即对任意的 $x$ 和 $x^{\prime }$ 有
  这时，接受分布 $\alpha (x,x^{\prime })$ 简化为
  
  
  • 随机游走 Metropolis 算法: Metropolis 选择的一个特例是令 $q(x,x^{\prime })=q(|x-x^{\prime }|)$。例如，
    
  • Metropolis 选择的另一个特例是 $q(x,x^{\prime })$ 取条件概率分布 $p(x^{\prime }|x)$，定义为多元正态分布，其均值是 $x$，其协方差矩阵是常数矩阵 (假设在给定 $x$ 的情况下，$x^{\prime }$ 的分布是均值为 $x$ 的多元正态分布。由于协方差矩阵为常数矩阵，因此 $q(x,x^{\prime })$ 是对称的)
• Metropolis 选择的特点是当 $x^{\prime }$ 与 $x$ 接近时，$q(x,x^{\prime })$ 的概率值高，否则 $q(x,x^{\prime })$ 的概率值低。状态转移在附近点的可能性更大

独立抽样

• 独立抽样: 假设 $q(x,x^{\prime })$ 与当前状态 $x$ 无关，即 $q(x,x^{\prime })=q(x^{\prime })$。建议分布的计算按照 $q(x^{\prime })$ 独立抽样进行。此时，接受分布 $\alpha (x,x^{\prime })$ 可以写成
  其中 $w(x^{\prime })=p(x^{\prime })/q(x^{\prime })$，$w(x)=p(x)/q(x)$
• 独立抽样实现简单，但可能收敛速度慢，通常选择接近目标分布 $p(x)$ 的分布作为建议分布 $q(x)$

Metropolis-Hastings 算法

单分量 Metropolis-Hastings 算法

• 基本思想: 在 Metropolis-Hastings 算法中， 通常需要对多元变量分布进行抽样，有时对多元变量分布的抽样是困难的。可以对多元变量的每一变量的条件分布依次分别进行抽样，从而实现对整个多元变量的一次抽样

---

• 假设马尔可夫链的状态由 $k$ 维随机变量表示
  其中 $x^{(i)}$ 表示马尔可夫链在时刻 $i$ 的状态
• 为了生成容量为 $n$ 的样本集合
  
  
  
  
  {
  
  
  
  x
  
  
  
  (
  
  
  1
  
  
  )
  
  
  
  
  ,
  
  
  .
  
  
  .
  
  
  .
  
  
  ,
  
  
  
  x
  
  
  
  (
  
  
  n
  
  
  )
  
  
  
  
  }
  
  
  
  {x^{(1)},…,x^{(n)} }
  
  
  {x(1),...,x(n)}，单分量 Metropolis-Hastings 算法由下面的 $k$ 步迭代实现 Metropolis-Hastings 算法的一次选代
  
  • 在第 $i$ 次迭代的第 $j$ 步 (也就是采样第 $i$ 个样本的第 $j$ 个属性)，对分量 $x_j$ 根据 Metropolis-Hastings 算法更新，得到其新的取值 $x_j^{(i)}$
  • 首先，由建议分布 $q(x_j^{(i-1)},x_j|x_{-j}^{(i)})$ 抽样产生分量 $x_j$ 的候选值 $x_j^{\prime (i)}$，这里 $x_{-j}^{(i)}$ 表示在第 $i$ 次迭代的第 $(j-1)$ 步后的 $x^{(i)}$ 除去 $x_j^{(i-1)}$ 的所有值，即
    其中分量 $1,2,...,j-1$ 己经更新
  • 然后，按照接受概率
    抽样决定是否接受候选值 $x_j^{\prime (i)}$。如果 $x_j^{\prime (i)}$ 被接受，则令 $x_j^{(i)}=x_j^{\prime (i)}$；否则令 $x_j^{(i)}=x_j^{(i-1)}$。其余分量在第 $j$ 步不改变。马尔可夫链的转移概率为
    

---

• 图 19.10 示意单分量 Metropolis-Hastings 算法的迭代过程。目标是对含有两个变量的随机变量 $x$ 进行抽样。如果变量 $x_1$ 或 $x_2$ 更新，那么在水平或垂直方向产生一个移动，连续水平和垂直移动产生一个新的样本点。注意由于建议分布可能不被接受，Metropolis-Hastings 算法可能在一些相邻的时刻不产生移动
  

吉布斯抽样 (Gibbs sampling)

• 吉布斯抽样 (Gibbs sampling) 可以认为是单分量 Metropolis-Hastings 算法的特殊情况，但是更容易实现，因而被广泛使用

满条件分布 (full conditional distribution)

• 满条件分布: 对于多元联合概率分布 $p(x)=p(x_1,...,x_k)$，其中 $x=(x_1,...,x_k{)}^T$ 为 $k$ 维随机变量。如果条件概率分布 $p(x_I|x_{-I})$ 中所有 $k$ 个变量全部出现，其中
  
  
  
  
  
  x
  
  
  I
  
  
  
  =
  
  
  {
  
  
  
  x
  
  
  i
  
  
  
  ,
  
  
  i
  
  
  ∈
  
  
  I
  
  
  }
  
  
  
  x_I={x_i,i\in I}
  
  
  xI​={xi​,i∈I}，
  
  
  
  
  
  x
  
  
  
  −
  
  
  I
  
  
  
  
  =
  
  
  {
  
  
  
  x
  
  
  i
  
  
  
  ,
  
  
  i
  
  
  ∉
  
  
  I
  
  
  }
  
  
  
  x_{-I} = {x_i ,i\notin I}
  
  
  x−I​={xi​,i∈/​I},
  
  
  
  
  I
  
  
  ⊂
  
  
  K
  
  
  =
  
  
  {
  
  
  1
  
  
  ,
  
  
  2
  
  
  ,
  
  
  .
  
  
  .
  
  
  .
  
  
  ,
  
  
  k
  
  
  }
  
  
  
  I\subset K = {1 ,2,…,k}
  
  
  I⊂K={1,2,...,k} 时，那么称这种条件概率分布为满条件分布
• 满条件分布有以下性质: 对任意的 $x,x^{\prime }\in \mathcal{X}$ 和任意的 $I\subset K$，有
  $$p\left(x_I\mid x_{-I}\right)=\frac{p(x)}{p(x_{-I})}=\frac{p(x)}{\int p(x)\mathrm{d}{x}_I}\propto p(x)$$而且有
  Metropolis-Hastings 算法中，可以利用上述性质，通过满条件分布概率的比计算联合概率的比，从而简化计算

---

例

$$N(1,(x_2-1{)}^{-2})=\frac{1}{\sqrt{2\pi }\sigma }\exp \left\{-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right\}=\frac{|x_2-1|}{\sqrt{2\pi }}\exp \left\{-\frac{(x-1{)}^2(x_2-1{)}^2}{2}\right\}$$

基本原理

• 其基本做法是，从联合概率分布定义满条件概率分布，依次对满条件概率分布进行抽样，得到样本的序列。可以证明这样的抽样过程是在一个马尔可夫链上的随机游走，每一个样本对应着马尔可夫链的状态，平稳分布就是目标的联合分布
  
  • 假设多元变量的联合概率分布为 $p(x)=P(x_1,...,x_k)$。吉布斯抽样从一个初始样本 $x^{(0)}=(x_1^{(0)},x_2^{(0)},...,x_k^{(0)}{)}^T$ 出发，不断进行迭代，每一次迭代得到联合分布的一个样本 $x^{(i)}=(x_1^{(i)},x_2^{(i)},...,x_k^{(i)}{)}^T$。最终得到样本序列
    
    
    
    
    {
    
    
    
    x
    
    
    
    (
    
    
    0
    
    
    )
    
    
    
    
    ,
    
    
    .
    
    
    .
    
    
    .
    
    
    ,
    
    
    
    x
    
    
    
    (
    
    
    n
    
    
    )
    
    
    
    
    }
    
    
    
    {x^{(0)},…,x^{(n)}}
    
    
    {x(0),...,x(n)}
  • 在每次迭代中，依次对 $k$ 个随机变量中的一个变量进行随机抽样。如果在第 $i$ 次迭代中，对第 $j$ 个变量进行随机抽样，那么抽样的分布是满条件概率分布 $p(x_j|x_{-j}^{(i)})$，这里 $x_{-j}^{(i)}$ 表示第 $i$ 次选代中，变量 $j$ 以外的其他变量

---

抽样第 $i$ 个样本 $x^{(i)}$ 的具体做法

• 在第 $i$ 步，首先对第一个变量按照以下满条件概率分布随机抽样得到 $x_1^{(i)}$
  $$p\left(x_1\mid x_2^{(i-1)},\cdots ,x_k^{(i-1)}\right)$$之后依次对第 $j$ 个变量按照以下满条件概率分布随机抽样得到 $x_j^{(i)}$
  最后对第 $k$ 个变量按照以下满条件概率分布随机抽样得到 $x_k^{(i)}$
  
• 于是得到整体样本
  $$x^{(i)}=(x_1^{(i)},x_2^{(i)},...,x_k^{(i)}{)}^T$$

---

吉布斯抽样是单分量 Metropolis-Hastings 算法的特殊情况

• 定义建议分布是当前变量 $x_j$，$j=1,2,...,k$ 的满条件概率分布
  这时，接受概率 $\alpha =1$,
  这里用到 $p(x_{-j})=p(x_{-j}^{\prime })$ 和 $p(•|x_{-j})=p(•|x_{-j}^{\prime })$ (吉布斯抽样中，每一次迭代的第 $j$ 步只改变样本的第 $j$ 个变量，因此 $x_{-j}=x_{-j}^{\prime }$)
• 转移核就是满条件概率分布
  也就是说依次按照单变量的满条件概率分布 $p(x_j^{\prime }|x_{-j})$ 进行随机抽样，就能实现单分量 Metropolis-Hastings 算法

---

吉布斯抽样 v.s. 单分量 Metropolis-Hastings 算法

• 吉布斯抽样对每次抽样的结果都接受，没有拒绝，这一点和一般的 Metropolis-Hastings 算法不同
• 吉布斯抽样适合于满条件概率分布容易抽样的情况，而单分量 Metropolis-Hastings 算法适合于满条件概率分布不容易抽样的情况，这时使用容易抽样的条件分布作建议分布

吉布斯抽样算法

---

抽样计算

• 吉布斯抽样中需要对满条件概率分布进行重复多次抽样。可以利用概率分布的性质提高满条件概率分布抽样的效率。下面以贝叶斯学习为例介绍这个技巧

---

• 设 $y$ 表示观测数据，$\alpha ,\theta ,z$ 分别表示超参数、模型参数、未观测数据，$x=(\alpha ,\theta ,z)$
  
• 贝叶斯学习的目的是估计后验概率分布 $p(x|y)$，求后验概率最大的模型
  
• 现在用吉布斯抽样估计后验概率分布 $p(x|y)$。吉布斯抽样中各个变量 $\alpha ,\theta ,z$ 的满条件分布有以下关系 (下面三式各自省略了上面式子中的无关项):
  其中 ${\alpha }_{-i}$ 表示变量 ${\alpha }_i$ 以外的所有变量，${\theta }_{{}_j}$ 和 $z_{-k}$ 类似。满条件概率分布与若干条件概率分布的乘积成正比，各个条件概率分布只由少量的相关变量组成 (概率图模型中相邻结点表示的变量)。所以，依满条件概率分布的抽样可以通过依这些条件概率分布的乘积的抽样进行。这样可以大幅减少抽样的计算复杂度，因为计算只涉及部分变量

参考文献

• 《统计学习方法》
• 马尔可夫链 (Markov Chains)