前言

感觉理论学的太多了，该做点实践操作一下了

第十四章

在我们之前学习的基于概率的模型中，进行预测时会涉及三类变量集合：需要推断的目标变量集合 $Y$、可观测到的变量集合 $O$，以及其他相关变量集合 $R$。无论模型采用生成式形式 $P(Y,R,O)$ 还是判别式形式 $P(Y,R|O)$，最终都需要得到我们关心的条件概率分布 $P(Y|O)$。为了得到这个分布，必须消除变量集合 $R$ 的影响。具体来说，这要求我们对 $R$ 中的所有变量执行求和或积分操作。需要注意的是，即使 $R$ 中的每个变量都是二值的（即只能取两种状态），计算过程中也需要处理高达 $2^{|Y|+|R|}$ 种可能的组合情况，这个计算量在实际应用中往往是不可行的。
为了解决这个问题，图论提供了一种直观的的语言来表达变量之间的条件独立性关系。因此一个复杂的联合概率分布，可以被分解为一系列更小的、更易于处理的局部概率关系的乘积。从而解决我们的维度灾难。

一般说来概率模型图有两大类：

• 贝叶斯网络：使用有向无环图来表示变量间的因果关系或依赖关系。箭头从“原因”指向“结果”。
• 马尔可夫随机场：使用无向图来表示变量间的软约束或相关性，但不强调方向性。

14.1 隐马尔可夫模型

隐马尔可夫模型（Hidden Markov Model, HMM）是动态贝叶斯网络的一个特例，专门用于建模时序数据或序列数据。该模型在语音识别、自然语言处理等众多领域中有广泛应用。

隐马尔可夫模型的结构可以通过以下图示直观表示：

       y_1 -------> y_2 -------> y_3 -------> ... -------> y_t
        |            |            |                       |
        |            |            |                       |
        v            v            v                       v
       o_1          o_2          o_3          ...        o_t

在这个图示中：

• $y_t$ 表示系统在时刻 $t$ 的隐藏状态，这些状态通常是无法直接观测到的。
• $o_t$ 表示在时刻 $t$ 基于隐藏状态 $y_t$ 生成的可观测值。

接下来，我们形式化地定义模型的各个组成部分：

1. 状态空间 Y={s1,s2,…,sN}\mathcal{Y} = \{s_1, s_2, \dots, s_N\}Y={s1​,s2​,…,sN​}：包含所有可能的隐藏状态，共有 $N$ 个不同的状态。
2. 观测空间 X={o1,o2,…,oM}\mathcal{X} = \{o_1, o_2, \dots, o_M\}X={o1​,o2​,…,oM​}：包含所有可能的观测值，共有 $M$ 个不同的观测结果。
3. 模型参数 $\lambda =[\mathbf{A},\mathbf{B},\pi ]$：
   • 初始状态概率向量 $\pi$：这是一个 $N$ 维向量，其中第 $i$ 个分量 ${\pi }_i$ 表示初始时刻（$t=1$）状态 $y_1$ 取值为 $s_i$ 的概率，即 ${\pi }_i=P(y_1=s_i)$，其中 $s_i\in \mathcal{Y}$。
   • 状态转移概率矩阵 $\mathbf{A}$：这是一个 $N\times N$ 的方阵，其中元素 $A_{ij}$ 表示在 $y_{t-1}=s_i$ 的条件下，$y_t=s_j$ 的概率，即 $A_{ij}=P(y_t=s_j|y_{t-1}=s_i)$。
   • 观测发射概率矩阵 $\mathbf{B}$：这是一个 $N\times M$ 的矩阵，其中元素 $B_{jk}$ 表示在 $y_t=s_j$ 的条件下，观测到 $o_k$ 的概率，即 $B_{jk}=P(x_t=o_k|y_t=s_j)$，其中 $s_j\in \mathcal{Y}$ 且 $o_k\in \mathcal{X}$。

图中蕴含了隐马尔可夫模型基于的两个关键假设：

1. 齐次马尔可夫性假设：当前状态 $y_t$ 仅依赖于前一个状态 $y_{t-1}$，而与更早的状态无关。数学上表示为：
   $$P(y_t=s_j|y_{t-1}=s_i,y_{t-2}=s_k,\dots ,y_1=s_l)=P(y_t=s_j|y_{t-1}=s_i)$$
   这正是马尔可夫链的核心性质。
2. 观测独立性假设：当前观测 $o_t$ 仅依赖于当前状态 $y_t$，而与历史状态和观测无关。数学表示为：
   $$P(x_t=o_k|y_t=s_j,y_{t-1}=s_i,\dots ,y_1=s_l,x_{t-1},\dots ,x_1)=P(x_t=o_k|y_t=s_j)$$

基于上述假设，我们可以推导出所有变量的联合概率分布：
$$P(x_1,y_1,\dots ,x_n,y_n)=P(y_1)P(x_1|y_1)\prod _{i=2}^{n}P(y_i|y_{i-1})P(x_i|y_i)$$
这个分布的推导是直接应用链式法则和模型假设的结果。

现在，我们详细说明如何使用HMM生成一个长度为 $n$ 的观测序列。

初始时刻（$t=1$）：
我们从初始概率分布 $\pi$ 中随机抽样，得到第一个隐藏状态 $y_1$。假设抽样结果为 $y_1=s_i$。然后，根据当前状态 $y_1=s_i$，我们查看观测概率矩阵 $\mathbf{B}$ 的第 $i$ 行 ${\mathbf{B}}_i$，该行给出了在状态 $s_i$ 条件下各个观测值的概率分布。从这个分布中随机抽样，得到第一个观测值 $o_1$。

迭代过程（对于 $t=2$ 到 $t=n$ 的每个时刻）：

1. 假设在时刻 $t-1$ 的状态为 $y_{t-1}=s_i$。
2. 为了确定时刻 $t$ 的状态：
   • 查看状态转移概率矩阵 $\mathbf{A}$ 的第 $i$ 行 ${\mathbf{A}}_i$，它表示从状态 $s_i$ 转移到其他状态的概率分布。
   • 根据 ${\mathbf{A}}_i$ 进行随机抽样，得到新状态 $y_t=s_j$。
3. 生成观测值：
   • 根据当前状态 $y_t=s_j$，查看 $\mathbf{B}$ 的第 $j$ 行 ${\mathbf{B}}_j$。
   • 从 ${\mathbf{B}}_j$ 分布中抽样得到观测值 $o_t$。

重复上述迭代步骤，直到处理完所有 $n$ 个时刻。最终，我们得到：

• 隐藏状态序列：$Y=(y_1,y_2,...,y_n)$
• 对应的观测序列：$O=(o_1,o_2,...,o_n)$

最后，我们探讨HMM能够解决的三个基本问题：

概率计算问题：给定一个完整的HMM模型参数 $\lambda =[\mathbf{A},\mathbf{B},\pi ]$ 和一个具体的观测序列 {x1,…,xn}\{x_1, \dots, x_n\}{x1​,…,xn​}，需要计算这个观测序列由该模型生成的总概率 $P(x_1,\dots ,x_n|\lambda )$。
解决方法是采用前向算法，该算法通过动态规划高效地计算观测序列的概率。

预测问题：在给定模型 $\lambda$ 和观测序列 {x1,…,xn}\{x_1, \dots, x_n\}{x1​,…,xn​} 的条件下，需要找到最有可能对应的隐藏状态序列 {y1∗,…,yn∗}\{y_1^*, \dots, y_n^*\}{y1∗​,…,yn∗​}。这个问题可以形式化为找到使条件概率 $P(y_1,\dots ,y_n|x_1,\dots ,x_n,\lambda )$ 最大的状态序列，即：
{y1∗,…,yn∗}=arg⁡max⁡y1∈Y,…,yn∈Y P(y1,…,yn∣x1,…,xn,λ) \{y_1^*, \dots, y_n^*\} = \underset{y_1 \in \mathcal{Y}, \dots, y_n \in \mathcal{Y}}{\arg\max} \, P(y_1, \dots, y_n | x_1, \dots, x_n, \lambda) {y1∗​,…,yn∗​}=y1​∈Y,…,yn​∈Yargmax​P(y1​,…,yn​∣x1​,…,xn​,λ)
解决方法是使用维特比算法，这是一种基于动态规划的方法，能够有效地找到最优路径。

学习问题：当只有一个观测序列 {x1,…,xn}\{x_1, \dots, x_n\}{x1​,…,xn​} 时，需要通过调整模型参数 $\mathbf{A}$、$\mathbf{B}$ 和 $\pi$，使得该观测序列的生成概率 $P(x_1,\dots ,x_n|\lambda )$ 最大化。
解决方法是采用鲍姆-韦尔奇算法，这是期望最大化算法在HMM中的具体应用，通过迭代方式优化模型参数。

14.2 马尔可夫随机场

隐马尔可夫模型是一种有向图模型，这自然表明节点之间存在单向的依赖关系。然而在现实问题中，许多关系并非单向，而是对称且相互影响的。对于这种关系，使用有向边来描述可能显得不自然甚至错误，因此我们需要一种无向图来刻画它，这就是马尔可夫随机场（MRF）。

现在，给定一个无向图 $G=(V,E)$，我们如何定义一个与该图结构相符的联合概率分布 $P(\mathbf{X})$？与贝叶斯网络不同，无向图中不存在“父节点”的概念，因此无法直接通过条件概率的乘积来构造联合分布。为了解决这个问题，我们需要采用一种新的建模方式，这引出了马尔可夫随机场（MRF）中的几个核心概念。

在构建联合概率分布时，MRF 不仅关注单个节点，更着眼于图中被称为团（clique）的子结构。团是图中节点的子集，满足该子集中任意两个不同的节点之间都存在一条边。特别地，单个节点本身也被视为一个团。由于团之间可能存在包含关系，我们尤其关注那些无法再扩展的团：如果一个团不能在保持团性质的前提下加入图中任何其他节点，则称其为极大团。

形式上说，对于 $n$ 个变量 x={x1,…,xn}\mathbf{x} = \{x_1, \dots, x_n\}x={x1​,…,xn​}，所有团构成的集合记为 $\mathcal{C}$。对于每个团 $Q\in \mathcal{C}$，其对应的变量集合记为 ${\mathbf{x}}_Q$。

为了对团内的变量依赖关系进行建模，我们引入势函数，也称因子，记为 ${\varphi }_Q({\mathbf{x}}_Q)$。该函数将团 $Q$ 中变量的每一个可能取值组合映射为一个非负实数，用于刻画这些变量之间的局部相互作用强度。这个分数代表了这种组合的兼容性：

• 分数值高：表示团内变量的这种取值组合非常和谐、兼容，更可能出现。
• 分数值低：表示这种组合不太和谐、不兼容，不那么可能出现。

有了这些概念后，我们下一步就是要把这些局部的势函数组合成一个全局的联合概率分布 $P(\mathbf{X})$。为此，我们需要引入Hammersley-Clifford 定理：任何一个严格为正的、并且其条件独立性关系可以被一个无向图 $G$ 描述的概率分布，必然可以被表示为图中所有极大团上的势函数的乘积形式。这个形式被称为吉布斯分布：
$$P(X_1,\dots ,X_n)=\frac{1}{Z}\prod _{C\in \mathcal{C}}{\varphi }_C(X_C)$$
其中 $Z$ 是归一化因子，也称为配分函数。它的计算方式是把所有可能的全局状态 $X$ 的“总分”加起来：
$Z={\sum }_X{\prod }_{C\in \mathcal{C}}{\varphi }_C(X_C)$……可见这其实是在计算所有变量的指数级数量的可能取值组合进行求和，这根本不可精确计算，所以在现实的算法中总会近似或者绕开它的计算，这里我们不展开它，后面会讨论。

接着我们看一下势函数，一般来说，我们常取指数函数作为势函数：
$${\varphi }_C(X_C)=e^{-E_C(X_C)}$$
其中 $E_C(X_C)$ 被称为能量函数。势函数值高，对应能量值低。常见形式为
$$E_C(X_C)=\sum _{u,v\in C,u\ne v}{\alpha }_{uv}{x}_u{x}_v+\sum _{v\in C}{\beta }_v{x}_v$$

其实马尔可夫随机场来源于物理中的统计力学，所以这里能看到很多类似的名词，然而我的统计力学学的有限，这里就不展开了。

最后，我们回顾一下之前提到的条件独立性。在马尔可夫随机场中，条件独立性十分直观。我们首先定义：集合 $C$ 分离集合 $A$ 和集合 $B$，如果从集合 $A$ 中任意一个节点出发，到集合 $B$ 中任意一个节点的所有路径，都必须经过集合 $C$ 中的至少一个节点。

全局马尔可夫性：如果一个概率分布 $P(X)$ 满足与图 $G$ 相关的全局马尔可夫性，那么对于任意三个节点集合 $A,B,C$，只要 $C$ 在图 $G$ 上分离了 $A$ 和 $B$，那么在概率分布 $P$ 中，变量集合 $X_A$ 和 $X_B$ 在给定变量集合 $X_C$ 的条件下是条件独立的。

全局马尔可夫性给出了 MRF 的独立性，其实全局马尔可夫性和以下两个性质完全等价：

• 局部马尔可夫性：任意一个变量，在给定其所有邻居节点的条件下，与图中所有其他非邻居节点都是条件独立的。
• 成对马尔可夫性：图中任意两个没有边直接相连的节点，在给定所有其他节点的条件下，是条件独立的。

14.3 条件随机场

在HMM中，我们依赖于齐次马尔可夫性和观测独立性这两个假设。这些假设在某些情况下可能过于严格，无法满足实际需求。例如，在自然语言处理中，一个词的标签往往依赖于整个句子的上下文，而不仅仅是当前单词；此外，当我们只关心判别任务时，MRF中的配分函数计算困难，并且需要对不必要的观察序列$X$进行建模。

为了解决这些问题，条件随机场（CRF）应运而生。它被设计为一种条件概率的判别式模型$P(Y|X)$，利用马尔可夫随机场（MRF）的框架，但条件是整个观察数据$X$。

CRF依旧是一个无向图$G=(V,E)$，其中节点代表变量，边代表依赖关系，但整个图都依赖于观察序列$X$。$y_v$ 表示与结点 $v$ 对应的标记变量，$n(v)$ 表示结点 $v$ 的邻接结点。如果图 $G$ 的每个变量 $y_v$ 都满足马尔可夫性，即P(yv∣x,yV∖{v})=P(yv∣x,yn(v))P(y_v | \mathbf{x}, \mathbf{y}_{V \setminus \{v\}}) = P(y_v | \mathbf{x}, \mathbf{y}_{n(v)})P(yv​∣x,yV∖{v}​)=P(yv​∣x,yn(v)​)，那么 $(\mathbf{y},\mathbf{x})$ 就构成一个条件随机场。

也就是说，给定观测变量 $\mathbf{x}$ 的条件下，对于图 $G$ 中的任意一个标记变量 $y_v$，其条件概率分布仅依赖于它的直接邻居 ${\mathbf{y}}_{n(v)}$，而与所有非邻居节点是条件独立的。最常见的图$G$是线性链结构，称为线性链CRF。
根据无向图模型的Hammersley-Clifford定理，一个马尔可夫随机场的概率分布可以被表示为其最大团的势函数的乘积。对于条件随机场，其条件概率分布 $P(\mathbf{y}|\mathbf{x})$ 也以类似的方式通过一系列特征函数和对应的权重来定义。
对于给定的观测序列 $\mathbf{x}=(x_1,x_2,\dots ,x_n)$ 和标记序列 $\mathbf{y}=(y_1,y_2,\dots ,y_n)$，线性链CRF的条件概率分布 $P(\mathbf{y}|\mathbf{x})$ 被定义为：
$P(\mathbf{y}|\mathbf{x})=\frac{1}{Z(\mathbf{x})}\exp \left({\sum }_{i=1}^n\left({\sum }_{j=1}^{K_t}{\lambda }_j{t}_j(y_{i-1},y_i,\mathbf{x},i)+{\sum }_{k=1}^{K_s}{\mu }_k{s}_k(y_i,\mathbf{x},i)\right)\right)$
其中：

• $t_j(y_{i-1},y_i,\mathbf{x},i)$是转移特征函数，这些特征函数关注相邻标记之间的依赖关系
• $s_k(y_i,\mathbf{x},i)$是状态特征函数，这些特征函数关注当前标记 $y_i$ 与当前位置 $i$ 的观测 ${\mathbf{x}}_i$ 之间的依赖关系
• $K_t$ 和 $K_s$ 分别是转移特征函数和状态特征函数的数量。
• 特征函数是CRF的核心要素和强大所在。它是一个实值函数，可以任意定义在整个观测序列 $\mathbf{x}$ 和整个标记序列 $\mathbf{y}$ 上，以捕捉它们之间的特定关联，以及 $\mathbf{y}$ 内部的依赖关系
  其他我们都很熟悉了，书上和这里直接就中断了。

14.4 学习和推断

之前三节分别介绍了HMM，MRF和CRF，他们都是概率图模型，接下来自然就是要使用他们去回答我们感兴趣的问题，得到这个回答的过程就是推断。

在概率的语境下，推断的过程实际上是将我们关心的变量保留下来，同时去除其他变量的过程。我们可以将这些变量$X$分为三类：

• 查询变量$X_F$：我们感兴趣的变量
• 证据变量$X_E$：我们观察到的变量
• 无关变量$X_Z$
  此时为了得到我们感兴趣问题的概率分布$P(X_F|X_E=x_e)$，我们必须对我们不关心的变量全部求和或积分，这个过程称为边缘化，以求和为例，即
  $$P(X_F|X_E=x_e)=\frac{P(X_F,X_E=x_e)}{P(X_E=x_e)}=\frac{{\sum }_{x_Z}P(X_F,X_E=x_e,X_Z=x_Z)}{{\sum }_{x_F}{\sum }_{x_Z}P(X_F,X_E=x_e,X_Z=x_Z)}$$

为了构建一个合适的概率图模型，我们自然需要对其参数进行学习。如果将参数变量也视为上述变量中的一部分，那么学习过程就和推断一致了（这也是书上没有单独撰写学习部分的原因）。另外，这里我也将无关变量表示进来。

可以看到，边缘化的困难在于不同变量的状态可能非常多。我们只关心 $X_Z$ 中的变量，假设每个变量的取值个数为 $k$ 个，那么我们就必须对 $k^{|X_Z|}$ 项进行计算，这显然是指数级的复杂度，属于经典的NP难题。
但即使面临这样的困难，我们仍然可以通过一些技巧来简化计算，而不必一开始就求助于近似算法。其中，我们将介绍两种精准算法：变量消除法和信念传播法。

变量消除法
这是最直观的一种方法。首先，我们考虑一个吉布斯分布 $P(X)=\frac{1}{Z}{\prod }_i{\varphi }_i(X_i)$，忽略归一化常数 $Z$，那么我们的边缘化计算可以写为：
$$\tilde{P}(X_F,X_E=x_e)=\sum _{x_Z}\left(\prod _i{\varphi }_i(\dots )\right)$$
这个乘积的组合数量极其庞大，如前所述，是指数级的，因此不能直接硬算。但图结构经常呈现出一些特点，例如叶子节点的存在。我们先考虑一个最简单的结构：一个待消除的变量 $A\in X_Z$ 只出现在一个因子中，比如说 ${\varphi }_A(A,B)$，其中 $B$ 是 $A$ 的邻居节点集合。现在，我们重写计算：
P~(XF,XE=xe)=∑xZ∖{A}∑A(ϕA(A,B)⋅∏j≠Aϕj)=∑xZ∖{A}((∏j≠Aϕj)⋅∑AϕA(A,B)) \tilde{P}(X_F, X_E=x_e) = \sum_{x_{Z} \setminus \{A\}} \sum_A \left( \phi_A(A,B) \cdot \prod_{j \neq A} \phi_j \right)=\sum_{x_{Z} \setminus \{A\}} \left( \left(\prod_{j \neq A} \phi_j \right) \cdot \sum_A \phi_A(A,B) \right) P~(XF​,XE​=xe​)=xZ​∖{A}∑​A∑​ ​ϕA​(A,B)⋅j=A∏​ϕj​ ​=xZ​∖{A}∑​ ​ ​j=A∏​ϕj​ ​⋅A∑​ϕA​(A,B) ​
可以发现，对 $A$ 的求和是相对简单的，并且之后再也没有变量 $A$ 了，它变成了一个依赖于 $B$ 的中间因子 $m_A(B)$。这就是消去法的原理——乘法分配律的应用。
通过多次进行相同的操作，我们就能得到较为简单的计算。但需要注意以下几点：：

• 计算顺序非常重要。例如，考虑一个星型网络，中间节点 $A$ 连接到周围的节点 $B、C、D$ 上。如果先消去 $A$，我们需要将所有与 $A$ 相关的因子相乘，生成一个巨大的、依赖于 {B,C,D}\{B,C,D\}{B,C,D} 的新因子。而如果我们先消除叶子节点 $B、C、D$，每次只会生成一个只依赖于 $A$ 的小因子。
• 对于树宽很大的图，中间因子会很大，算法往往不可行
• 如果换一个查询变量$X_F^{\prime }$，那么要从头算起——这一点可以通过下面的信念传播法处理

信念传播法
信念传播算法将每个变量视为一个独立的计算单元。这些变量通过向邻居节点发送消息来共享信息。一条消息 $m_{i\to j}(X_j)$ 可以直观地理解为：节点 $i$ 基于其所掌握的所有信息，向节点 $j$ 传达它认为 $j$ 可能处于的状态。当消息传递过程达到稳定状态后，每个节点都可以根据接收到的消息计算出自己的边缘概率分布，这被称为信念。通过这种方式，我们能够高效地完成所有变量的边缘概率计算，而无需像变量消除法那样进行重复计算。

在变量消除法中，消除变量的过程在信念传播法中对应着消息的传递。具体来说，消息的计算公式如下：
mi→j(Xj)=∑Xiψi(Xi)ψij(Xi,Xj)(∏k∈N(i)∖{j}mk→i(Xi)) m_{i \to j}(X_j) = \sum_{X_i} \psi_i(X_i) \psi_{ij}(X_i, X_j) \left( \prod_{k \in N(i) \setminus \{j\}} m_{k \to i}(X_i) \right) mi→j​(Xj​)=Xi​∑​ψi​(Xi​)ψij​(Xi​,Xj​) ​k∈N(i)∖{j}∏​mk→i​(Xi​) ​
其中：

• ∏k∈N(i)∖{j}mk→i(Xi)\prod_{k \in N(i) \setminus \{j\}} m_{k \to i}(X_i)∏k∈N(i)∖{j}​mk→i​(Xi​) 表示节点 $i$ 从除了 $j$ 以外的所有其他邻居那里收到的消息，这里 $N(i)$ 是节点 $i$ 的所有邻居集合。这部分汇总了除了来自 $j$ 的信息外，$i$ 所能获取的关于自身状态的所有外部信息。
• ${\psi }_i$ 是节点势函数，它蕴含节点本身的信息。特别地，如果变量 $X_i$ 是证据变量且观测值为 $x_i^{\prime }$，那么 ${\psi }_i(X_i=x_i^{\prime })=1$，而对所有其他状态 ${\psi }_i(X_i\ne x_i^{\prime })=0$。
• ${\psi }_{ij}(X_i,X_j)$ 是边势函数，它代表变量 $X_i$ 和 $X_j$ 之间的依赖关系，即之前的因子。

一旦所有消息都计算完毕，任何一个节点 $i$ 的信念 $b_i(X_i)$，即其未归一化的边缘概率，可以通过将其本地势函数与所有收到的消息相乘得到：
$$b_i(X_i)\propto {\psi }_i(X_i)\prod _{k\in N(i)}{m}_{k\to i}(X_i)$$
正如之前所述，该算法要求图中不存在环结构。那么在有环的情况下，这个算法还能有效吗？事实上我们知道，此时会出现消息循环传递的情况，因此该算法不仅不再精确，也不保证能够收敛——但也可能在某些情况下达到收敛状态。此时算法就转变为一种近似推断方法，此时称为循环信念传播。事实上，显示运用中近似推断才是我们的主流，接下来我们就进入近似推断的部分。

14.5 近似推断

精准推断虽然在理论上很理想，但正如之前所述，其计算量往往非常巨大，以至难以获得实际结果。在这种情况下，我们就需要寻求高质量的近似推断方法。

近似推断通常有两种思路：

• 通过采样来估计我们感兴趣的统计量，这种方法并不直接处理概率分布本身
• 通过简单的概率分布逼近真实的复杂分布

下面我们分别介绍这两种方法中的代表性技术。

马尔可夫链蒙特卡洛
马尔可夫链蒙特卡洛（MCMC）是一种基于采样和大数定律的统计方法。例如，当我们无法直接计算后验概率 $P(X_F|X_E=x_e)$ 时，如果能够从中抽取样本 $x_F^{(1)},x_F^{(2)},\dots ,x_F^{(N)}$，那么：

• 状态 $a$ 的概率可以近似表示为 $P(X_F=a|X_E=x_e)\approx \frac{1}{N}{\sum }_{i=1}^N\mathbb{I}(x_F^{(i)}=a)$
• 函数 $f$ 的期望可以近似表示为 $\mathbb{E}[f(X_F)|X_E=x_e]\approx \frac{1}{N}{\sum }_{i=1}^{N}f(x_F^{(i)})$

这里大数定律发挥了关键作用。然而，如果没有具体的概率分布，我们如何进行采样呢？MCMC 方法指出，可以构建一个马尔可夫链，使其平稳分布恰好等于目标概率分布，从而通过采样获得所需数据。暂且不论这个链的具体构造方法，我们先考虑链达到平稳状态的条件。在马尔可夫链理论中，我们多次接触到状态转移的概念：对于任意两个状态 $x$ 和 $x^{\prime }$，如果用 $T(x^{\prime }|x)$ 表示从状态 $x$ 转移到 $x^{\prime }$ 的概率，用 $p(x)$ 表示平稳分布，那么当分布达到平稳时，必须满足如下条件：
$$p(x)T(x^{\prime }|x)=p(x^{\prime })T(x|x^{\prime })$$

这个条件被称为细致平衡条件，是马尔可夫链达到平稳状态的充分必要条件。现在我们构造链的动机就清晰了：我们需要设计一个好的$T$使得分布能平稳在目标分布上。

MCMC方法中的一种重要算法是Metropolis-Hastings（MH）算法，它基于“拒绝采样”的思想来逼近目标分布 $p$。简单来说，该算法的核心步骤如下：

• 首先引入一个简单的提议分布 $Q(x^{\prime }|x)$，这个分布应该容易从中抽样。给定当前状态 $x^{(t-1)}$，我们从提议分布中抽取一个候选状态 $x^{\ast }$：$x^{\ast }\sim Q(x|x^{(t-1)})$。
• 然后以一定的接受概率 $A(x^{\ast }|x^{(t-1)})$ 来决定是否接受这个提议。如果接受，则链转移到新状态 $x^{(t)}=x^{\ast }$；如果拒绝，则链保持原状态 $x^{(t)}=x^{(t-1)}$。

因此，从 $x^{(t-1)}$ 转移到 $x^{\ast }$ 的总概率是提议概率乘以接受概率：$P(x^{\ast }|x^{(t-1)})=Q(x^{\ast }|x^{(t-1)})A(x^{\ast }|x^{(t-1)})$。代入细致平衡条件，我们有：
$$p(x^{(t-1)})Q(x^{\ast }|x^{(t-1)})A(x^{\ast }|x^{(t-1)})=p(x^{\ast })Q(x^{(t-1)}|x^{\ast })A(x^{(t-1)}|x^{\ast })$$

现在，问题转化为如何设计接受概率 $A$。整理上述方程可得：
$$\frac{A(x^{\ast }|x^{(t-1)})}{A(x^{(t-1)}|x^{\ast })}=\frac{p(x^{\ast })Q(x^{(t-1)}|x^{\ast })}{p(x^{(t-1)})Q(x^{\ast }|x^{(t-1)})}$$

为了使接受概率尽可能大，通常选择：
$$A(x^{\ast }|x^{(t-1)})=\min \left(1,\frac{p(x^{\ast })Q(x^{(t-1)}|x^{\ast })}{p(x^{(t-1)})Q(x^{\ast }|x^{(t-1)})}\right)$$

可以发现，所有概率项都以相同的形式出现在分子和分母中，这意味着配分函数 $Z$ 可以被约掉，这正是MH算法强大之处。现在我们看一下具体的算法：
初始化: 随机选择一个初始状态 $x^{(0)}$。
开始迭代：对于 $t=1,2,\dots ,N$：

• 从提议分布中抽取候选样本：$x^{\ast }\sim Q(x|x^{(t-1)})$。
• 计算接受率：$$r=\frac{p(x^{\ast })Q(x^{(t-1)}|x^{\ast })}{p(x^{(t-1)})Q(x^{\ast }|x^{(t-1)})}$$
• 从均匀分布中抽取随机数：$u\sim U[0,1]$。
  • 如果 $u<\min (1,r)$，则接受提议：$x^{(t)}=x^{\ast }$。
  • 否则，拒绝提议：$x^{(t)}=x^{(t-1)}$。
• 扔掉前B项：收集样本 {x(B+1),x(B+2),…,x(N)}\{x^{(B+1)}, x^{(B+2)}, \dots, x^{(N)}\}{x(B+1),x(B+2),…,x(N)}。

最后扔掉前面 $B$ 项的原因在于链的初始阶段可能尚未达到平稳分布，因此这些样本不够稳定。

可以证明，吉布斯抽样是MH算法的一个特例，其提议分布恰好是真实的条件分布 $Q(x^{\prime }|x)=p(x^{\prime }|x_{-i})$，并且在这种特殊情况下，接受率 $r$ 恒等于1。因此，吉布斯抽样可以看作是一个“从不拒绝”的MH算法。

变分推断
变分推断采用了一种不同的思路：我们不再直接计算复杂的后验分布，而是从一个易于处理的简单分布族 $Q$ 中找到一个分布 $q(Z)\in Q$，使其与真实后验分布 $p(Z|X)$ 尽可能接近。

为了衡量这种接近程度，我们通常使用KL散度作为指标。我们的目标是找到最优的 $q^{\ast }(Z)$，使得KL散度最小：
$$q^{\ast }(Z)=\arg \underset{q(Z)\in Q}{\min }{D}_{KL}(q(Z)||p(Z|X))$$
其中KL散度定义为：
$$D_{KL}(q(Z)||p(Z|X))=\int q(Z)\ln \frac{q(Z)}{p(Z|X)}dZ$$
然而，KL散度本身依赖于 $p(Z|X)$，这恰恰是我们无法直接计算的。为了解决这个问题，我们转而考虑对数证据 $\ln p(X)$，并对其进行如下变换：
$$\begin{array}{rl}\ln p(X)& =\ln \int p(X,Z)dZ\\ & =\ln \int p(X,Z)\frac{q(Z)}{q(Z)}dZ\\ & =\ln {\mathbb{E}}_{q(Z)}\left[\frac{p(X,Z)}{q(Z)}\right]\\ & \ge {\mathbb{E}}_{q(Z)}\left[\ln \frac{p(X,Z)}{q(Z)}\right]\end{array}$$
这里最后一步使用了琴生不等式。我们定义不等式右侧的项为证据下界（ELBO），记作 $\mathcal{L}(q)$：
$$\mathcal{L}(q)={\mathbb{E}}_{q(Z)}\left[\ln \frac{p(X,Z)}{q(Z)}\right]=\int q(Z)\ln \frac{p(X,Z)}{q(Z)}dZ$$

接下来，我们通过代数运算将KL散度与ELBO联系起来：
$$\begin{array}{rl}{D}_{KL}(q(Z)||p(Z|X))& =\int q(Z)\ln \frac{q(Z)}{p(Z|X)}dZ\\ & =\int q(Z)\ln \left(\frac{q(Z)}{p(X,Z)/p(X)}\right)dZ\\ & =\int q(Z)\left(\ln q(Z)-\ln p(X,Z)+\ln p(X)\right)dZ\\ & =\int q(Z)\ln \left(\frac{q(Z)}{p(X,Z)}\right)dZ+\int q(Z)\ln p(X)dZ\\ & =-\mathcal{L}(q)+\ln p(X)\int q(Z)dZ\\ & =-\mathcal{L}(q)+\ln p(X)\end{array}$$
于是我们得到：
$$\ln p(X)=\mathcal{L}(q)+D_{KL}(q(Z)||p(Z|X))$$

由于在给定观测数据 $X$ 时，$\ln p(X)$ 是一个常数，且KL散度非负，因此最小化KL散度等价于最大化ELBO $\mathcal{L}(q)$。而 $\mathcal{L}(q)$ 只涉及联合概率 $p(X,Z)$（在推断问题中已知）和近似分布 $q(Z)$，这使得问题变得可行。

接下来，我们需要选择近似分布族 $Q$。最常用的方法是平均场假设，该假设认为隐变量 Z={z1,z2,…,zM}Z=\{z_1, z_2, \dots, z_M\}Z={z1​,z2​,…,zM​} 在近似分布 $q(Z)$ 中是相互独立的，因此 $q(Z)$ 可以分解为：
$$q(Z)=\prod _{i=1}^M{q}_i(z_i)$$
在平均场假设下，我们通过坐标上升法来最大化 $\mathcal{L}(q)$。具体步骤如下：固定除了一个因子 $q_j(z_j)$ 之外的所有其他因子 $q_{i\ne j}(z_i)$，然后优化 $\mathcal{L}(q)$ 以找到最优的 $q_j^{\ast }(z_j)$。我们对所有 $j$ 轮流进行这一过程，直到收敛。

1. 将 $q(Z)=q_j(z_j){\prod }_{i\ne j}{q}_i(z_i)$ 代入 $\mathcal{L}(q)$。
2. 将 $\mathcal{L}(q)$ 中与 $q_j(z_j)$ 无关的项合并为常数 $\text{const}$。
3. 经过整理，与 $q_j(z_j)$ 相关的部分为：
   $$\mathcal{L}(q)=\int q_j(z_j)\left({\mathbb{E}}_{i\ne j}[\ln p(X,Z)]\right)d{z}_j-\int q_j(z_j)\ln q_j(z_j)d{z}_j+\text{const}$$
   其中 ${\mathbb{E}}_{i\ne j}[\cdot ]$ 表示对 ${\prod }_{i\ne j}{q}_i(z_i)$ 求期望。我们定义 $\ln \tilde{p}(z_j)={\mathbb{E}}_{i\ne j}[\ln p(X,Z)]$。
4. 上式可以表示为 $-D_{KL}(q_j(z_j)||\tilde{p}(z_j))+\text{const}$。为了最大化 $\mathcal{L}(q)$，我们需要最小化KL散度，这当且仅当 $q_j(z_j)\propto \tilde{p}(z_j)$ 时成立。
   最终，我们得到更新方程：
   $$\ln q_j^{\ast }(z_j)={\mathbb{E}}_{i\ne j}[\ln p(X,Z)]+\text{const}$$

14.6 话题模型

<优化后的文本>
现在我们来看一个具体的领域——主题模型。我们每天面对成千上万的文章，有没有办法通过统计手段自动了解某篇文章的主要内容呢？
最朴素的方法就是统计词频，但遇到多领域混合的文章，例如关于“基因编辑伦理”的文章，这种方法可能只识别出伦理相关的词汇，而忽略基因编辑部分。
主题模型就是为了解决这个问题而提出的。它的核心思想是：一篇文档的内容不是由一个单一主题决定的，而是由多个主题以不同的比例混合而成。而每个主题本质上是一个关于词汇的概率分布。
其中，主题模型中最经典的方法称为隐狄利克雷分配 (Latent Dirichlet Allocation, LDA)。LDA 的目标是，在只给定文档和其中词语的情况下，自动推断出这些隐藏的主题结构：

• 每个主题是什么？即每个主题下，哪些词的出现概率高。
• 每篇文档是怎样由这些主题构成的？即每篇文档中，各个主题的比例是多少。
  在深入模型细节之前，我们必须明确几个基本定义：
• 词：文本处理的基本单元。
• 词汇表：数据集中所有不重复的词构成的集合。
• 文档：一篇由多个词构成的文本。
• 语料库：所有文档的集合。
• 主题：一个核心的隐变量。它不是单个词，而是一个覆盖整个词汇表的概率分布。例如，“科技”主题可能表现为 {p("芯片")=0.05,p("算法")=0.04,…,p("香蕉")=0.00001,… }\{ p(\text{"芯片"})=0.05, p(\text{"算法"})=0.04, \dots, p(\text{"香蕉"})=0.00001, \dots \}{p("芯片")=0.05,p("算法")=0.04,…,p("香蕉")=0.00001,…}。

此外，LDA 采用词袋假设，这意味着我们忽略一篇文章中词语的顺序和语法，只关心每个词出现了多少次。

现在假设我们有一个固定的词汇表和预设的 $K$ 个主题。要生成语料库中的第 $d$ 篇文档，步骤如下：

• 确定文档的主题分布：为这篇文档 $d$ 从狄利克雷分布中抽取一个 $K$ 维的主题比例向量 ${\theta }_d$。例如，${\theta }_d=(0.7,0.2,0.1,\dots )$ 表示这篇文档 70% 关于主题1，20% 关于主题2，等等。
• 生成文档中的每个词：对于文档 $d$ 中的第 $n$ 个词（假设文档长度为 $N_d$）：
  • 为这个词选择一个主题：根据主题分布 ${\theta }_d$，抽取一个具体的主题 $z_{d,n}$。
  • 根据选定的主题生成一个词：假设在上一步抽到了主题 $k$（即 $z_{d,n}=k$），则根据主题 $k$ 的词汇分布 ${\beta }_k$ 抽取一个具体的词 $w_{d,n}$。这里，${\beta }_k$ 是一个定义在整个词汇表上的概率分布，它也是从一个狄利克雷分布 $\text{Dir}(\eta )$ 中预先生成的。

这样就能生成一篇文档了。如果要生成整个语料库的所有文档，就需要重复上述过程 $M$ 次（假设有 $M$ 篇文档）。现在我们的目标是计算后验分布：
$$p(W,Z,B,\Theta |\alpha ,\eta )=\underset{\text{步骤 0: 生成所有主题}}{\underbrace{\left(\prod _{k=1}^{K}p({\beta }_k|\eta )\right)}}\times \underset{\text{步骤 1: 循环生成所有文档及其中的词}}{\underbrace{\left(\prod _{d=1}^{M}p({\theta }_d|\alpha )\left(\prod _{n=1}^{N_d}p(z_{d,n}|{\theta }_d)p(w_{d,n}|z_{d,n},B)\right)\right)}}$$
其中：

• $W$：观测到的语料库中的所有词（数据）。
• $\Theta$：所有文档的主题分布集合 {θd}d=1M\{\theta_d\}_{d=1}^M{θd​}d=1M​。
• $B$：所有主题的词汇分布集合 {βk}k=1K\{\beta_k\}_{k=1}^K{βk​}k=1K​。
• $Z$：所有词的主题分配集合 {zd,n}\{z_{d,n}\}{zd,n​}。
• $\alpha ,\eta$：模型的超参数。

而这个后验分布的难点在于边缘化，因此我们采用变分推断来近似它。我们使用平均场假设来构建一个简单的近似后验分布 $q$：
$$q(\Theta ,B,Z|\gamma ,\lambda ,\varphi )=\prod _{d=1}^{M}q({\theta }_d|{\gamma }_d)\prod _{k=1}^{K}q({\beta }_k|{\lambda }_k)\prod _{d=1}^M\prod _{n=1}^{N_d}q(z_{d,n}|{\varphi }_{d,n})$$
其中 $\lambda$、$\gamma$、$\varphi$ 是三种因子对应的变分参数。我们的目标是通过调整这些变分参数 $\gamma ,\lambda ,\varphi$ 来最大化 ELBO（证据下界），后续操作已在变分推断一节中说明。