Advanced Topics in Statistical Machine Learning 笔记06：Markov Chains and Hidden Markov Models

我们在前面几篇中已经了解了贝叶斯网络，它是有向图 (Directed Graph) 的概率图模型 (Probabilistic Graphic Model)。本节中，我们主要来了解马尔科夫链 (Markov Chains) 以及隐马尔可夫模型 (Hidden Markov Models, HMM)

马尔科夫链 (Markov Chains) 以及隐马尔可夫模型 (Hidden Markov Models, HMM) 属于动态贝叶斯网络 (Dynamic Bayesian Networks (DBN))。它们都有着随时间或者空间增长的重复结构，同时这种结构很简单并且使用了马尔科夫性质。

所谓的马尔科夫性质 (Markov Property) 声明了在已知现在状态 (Current State) 的情况下，未来的状态 (Future States) 和过去的状态 (Past States) 互相独立。在马尔科夫链 (Markov Chain) 中，所有的状态都是可以被直接观察到的，而隐马尔可夫模型 (HMM) 通过允许存在隐藏状态 (Hidden State) 扩展了马尔科夫链。

到目前为止，我们所讨论的都是静态的贝叶斯网络，也就是网络的结构不会发生变化，但在实际问题中，我们往往需要动态的模型，比如下图中的这样一个序列：

这个序列随着时间或者空间会发生动态的变化，比如不断有新的节点加入，它可以用来为那些动态问题建模，比如语音识别，机器人定位等。换句话说，动态贝叶斯网络（DBN）允许在模型中纳入时间或空间，而最简单的动态贝叶斯网络就是马尔科夫链以及隐马尔可夫模型。

马尔科夫链（Markov Chains）

我们先来看马尔科夫链 (Markov Chains)。马尔科夫链是一个状态机 (State Machine) ，在马尔科夫链中，X 是离散的变量，而变量的每一个取值就被称为状态 (State) ，状态之间的转换是不确定的 (Nondeterministic)。 马尔科夫链本质上是一个非常简单的隐马尔科夫模型。马尔科夫链也有其自己的应用，主要体现在图模型，它是用于近似推理的采样算法（蒙特卡洛马尔科夫链）的一个基本组成部分。

在马尔科夫链中，我们主要有两个参数：

• 先验概率 (Prior Probability): P(X₁)
• 转换概率 (Transition Probability) / 动态 (Dynamic): P(X_t|X_t-1)

这两个参数很好理解，先验概率就是在最初时刻出现某个状态的概率，而转换概率用来表示每个 Time Step 状态发生变化的概率。我们需要注意的是，我们认为转换概率在任何时候都是相同的，因此，这是一个时间均一链 (Time-homogeneous Chain)

下面我们用一个具体的例子来进行理解。我们用马尔科夫链来表示连续几天的天气情况，因此变量 X 的取值（状态）为 {Sun, Rain}。

此外，我们已知第一天为晴天，即 X₁ = Sun，因此初始概率分布 (Initial Distribution) 为：

接下来定义转换概率，即每天天气变化的概率表：

这里我们用 3 中方式 (Table, Matrix, Diagram) 来表示转换概率。简单的来说就是已知前一天为晴天/雨天，第二天为晴天/雨天的概率。

和贝叶斯网络一样，马尔科夫链同样也包含各种独立性关系，我们可以像之前一样用 D-Separation 方法来寻找，但是，我们在之前已经说过，不管是马尔科夫链还是马尔科夫模型，它们都遵循马尔科夫性质，即 “在已知现在状态 (Current State) 的情况下，未来的状态 (Future States) 和过去的状态 (Past States) 互相独立”。因此，我们能直接给出马尔科夫链中的独立关系表达式：X_t+1 ⊥ X_t-1 | X_t 。因此，在上面的例子中，我们很容易就能得到以下独立关系：𝑋₃ ⊥ 𝑋₁|𝑋₂ , 𝑋₄ ⊥ 𝑋₂|𝑋₃。当然，我们也能够使用贝叶斯网络的链式法则和一般链式法则从数学角度去证明它的合理性，这里我们只给出简单的证明示意：

接下来，我们来看马尔科夫链的参数量。首先我们回忆一下，如果要使用一个联合概率表来表示，需要的参数量，如果用 |X| 表示变量的取值数量，那么联合概率表需要的参数量就是 |X|ⁿ ，其中 n 是变量的数量。接下来看马尔科夫链，我们可以用以下式子表示马尔科夫链：

他需要的参数量为 |X| + (n - 1) |X|² 。前半部分表示定义初始分布需要的参数，即 P(X₁) 需要的参数，后半部分表示除了初始分布以外，剩余所有状态需要的参数量。但是因为我们之前已经定义了转换概率不会因时间的变化而变化，所以，实际上参数量应为 |X| + |X|² ，前半部分就是初始分布的参数，后半部分就是转换概率表的参数。 这再一次表明，我们只需要定义了初始分布以及转换概率，就得到了整个马尔科夫链。

下面我们就可以对马尔科夫链进行计算。最简单的计算就是求取整个状态序列的概率，这非常直观，直接根据链式法则就能进行计算。 比如我们想知道上面例子中的序列 {Sun, Rain, Rain, Sun, Sun} 的概率：

下面我们来看稍微复杂一些的计算，保持某种状态的概率 (Staying in a Certain State) ，即在前 d 个 Time Step 都为某种状态的概率。我们来看一个具体一些的例子，比如一个前三天全都是雨天的序列 {Rain, Rain, Rain, Sun}，那么就可以进行如下计算：

为什么这里的指数是 d - 1，因为在 d 个 Time Step 之间，只进行了 d - 1 次转换。所以，上面的例子中，我们能够得到如下结果：

第三类计算是求处于某个状态的预期时间 (Expected number of time step)。比如，在上面的天气模型中，我们想要计算预期的下雨天数是多少？这里我们不做详细的介绍，只是给出计算的公式：

Mini-forward Algorithm

我们怎么计算某一天 t 的状态概率 P(X) 呢？答案是我们可以模拟一个马尔科夫链来回答这个问题。我们已经知道，马尔科夫链一定有一个初始分布 (Initial Distribution) （也就是先验概率）以及转换概率 (Transition Probability)。根据这两个基本信息，可以对 P(x_t) 进行计算：

因此 P(x_t = Sun) = P(x_t = Sun, x_t-1 = Sun) + P(x_t = Sun, x_t-1 = Rain) = P(x_t = Sun|x_t-1 = Sun)*P(x_t-1 = Sun) + P(x_t = Sun|x_t-1 = Rain)*P(x_t-1 = Rain)。也就是说，我们只需要知道前一天的状态概率以及转换概率就能得到某一天的状态，这本质上就是一个状况分析 (Case Analysis)，我们考虑前一天是雨天的情况以及前一天是晴天的情况。

具体一点来看，我们仍用之前的天气例子：

此时如果我们计算 P(x₃ = Sun) = P(x₃ = Sun|x₂ = Sun) * P(x₂ = Sun) + P(x₃= Sun|x₂ = Rain)*P(x₂ = Rain) = 0.9 * 0.9 + 0.3 * 0.1 = 0.84

该算法就被称为 Mini-forward 算法：

该算法会将 Time Step 1 ~ n 的所有状态概率全部计算出来。它的复杂度为 O(n|X|²)

Stationary Distributions

我们仍使用前面的天气例子，不过现在我们考虑第一天是晴天和第一天是雨天两种情况，我们针对这两种情况都去使用 Mini-forward 算法来计算序列的状态，可以得到如下结果

这里不管我们给先验概率分配怎样的数值，最后都会变为 Sun=0.75, Rain=0.25 这个概率分布。这并不是特例，事实上，对于很多马尔科夫链来说，初始概率分布 (Initial Distribution) 的影响都会随着时间越来越小，最终都会向 Stationary Distribution 收敛，这个 Stationary Distribution 与初始分布相互独立。

我们可以这样来定义 Stationary Distribution ：马尔科夫链的 Stationary Distribution π 是最终收敛的结果，它满足：

那么我们如何计算这个 Stationary Distribution？只需要使用上面的公式 2 以及转换概率表即可：

这里再给出一个蚂蚱的例子：

根据上图，我们可以得到转换概率矩阵：

根据矩阵，我们能得到以下方程组：

π(-2) = 0.75 * π(-2) + 0.25 * π(-1)
π(-1) = 0.25 * π(-2) + 0.5 * π(-1) + 0.25 * π(0)
π(0) = 0.25 * π(-1) + 0.5 * π(0) + 0.25 * π(1)
π(1) = 0.25 * π(0) + 0.5 * π(1) + 0.25 * π(2)
π(2) = 0.25 * π(1) + 0.75 * π(2)

因此可以得到它的 Stationary Distribution 为：
π(-2) = π(-1) = π(0) = π(1) = π(2) = 1/5

Irreducible Markov Chains（不可约马尔科夫链）

不可约马尔科夫链指的是，在该马尔科夫链中，所有状态 x’ 都能从另一个状态 x 到达，这表明 P(X_t=x’|X₁=x) > 0。这也被称为正则链 (Regular Chain)。下图中是一个反例，它是可约的马尔科夫链 (Reducible Markov Chains)，因为我们从 State 1 出发，能够到达 2，但是无法到达 3 和 4。

这里我们再次强调，上图不是一个贝叶斯网络，节点并不是变量，而是马尔科夫链中变量的取值，也就是状态 (State) ！！

每一个马尔科夫链都至少有一个 Stationary Distribution，不可约的马尔科夫链 (Irreducible Markov Chains) 能够保证只有唯一的 (Unique) Stationary Distribution 但是，一个不可约的马尔科夫链不一定会收敛到它的 Stationary Distribution。为了保证收敛，我们需要另一个属性：非周期性 (Aperiodic)

Aperiodic Markov Chains（非周期性马尔科夫链）

如果可以随时返回任何状态 (State) ，则马尔可夫链是非周期性 (Aperiodic) 的。 下面看一个反例：

首先这个马尔科夫链确实是 Irreducible 的，因为我们无论从哪个 State 出发，都能到达其他任意的 State。这就意味它确实有一个唯一的 Stationary Dictionary

至于周期性怎么理解，我们需要分开看每个 Time Step。

T / D
1：1
2：2
3：1,3
4：2,4
5：1,3
6：2,4

很显然，这个情况下，若 Time Step 为奇数，都是1,3；若是偶数，都是2,4。这就违反了非周期性的定义：在任何时间 t，都有可能为任何状态。因为这时，凡是奇数 Time Step，都不可能去2, 4

一个不可约且非周期性的马尔科夫链会收敛到唯一的 Stationary Distribution。

• 不可约 (Irreducible)：我们能够从任意状态到达任意状态
• 非周期 (Aperiodic)：避免链条在状态之间永远交替，而不会固定在稳态状态中

Markov chains Convergence（马尔科夫链收敛）

尽管一个不可约且非周期性的马尔科夫链会收敛到唯一的 Stationary Distribution，但收敛的速度可能会非常慢。

上面的这个马尔科夫，它最终的 Stationary Distribution 应该是 <0.5, 0.5>。但是如果 𝜖 很小，收敛的速度会非常慢，这很好理解，𝜖 很小就意味着 1 - 𝜖 会很大，这样就会有很高的概率保持某个状态，很少转变到另一个状态。

Applications of Markov Chains（马尔科夫链的应用）

最著名的应用有以下几个：

• 马尔可夫链蒙特卡罗（MCMC）是一种强大的近似推理算法，用于Stan等统计软件
• 马尔科夫链是7zip中使用的 Lempel-Ziv-Markov压缩算法（LZMA）的一部分
• Google 1.0使用的PageRank算法是马尔科夫链的直接应用

这里对 PageRank 稍作一些介绍。它用状态图对网络进行建模，网页 (Pages) 是状态 (Sate)，超链接 (Hyperlink) 是转换 (Transitions)。从状态 i 的每一次转换概率为 α / k_i ，这里的 α 是一个常数，一般会选择 0.85，表示有效连接的权值是0.85，剩下15%给是 Phantom Link 的权值；k_i 则是节点 i 的出度 (Outgoing Degree)。网页与查询的相关性与指向该网页的边的数量成正比。但是这样的模型很容易被欺骗，比如有一个相关度很低的网页 E，它可以被伪装成相关度很高的样子，可以建立一系列的伪页面，这些页面可能和其他的一些页面也相连，它们同时与 E 相连，这样会使得连向 E 的边数大大增加。所以 PageRank 的核心思想是使结果更加可靠。

隐马尔可夫模型（Hidden Markov Models, HMM）

隐马尔可夫模型 (HMM) 就是状态 (States) 无法被直接观察的马尔科夫链 (Markov Chain)。我们还是用天气的例子来理解，比如我们在一个封闭的实验室，无法观察到室外的情况，只能根据各种传感器来判断具体的情况怎样，比如湿度传感器，温度传感器，或是观察进入实验室的人有没有带雨伞等。因此，从这里我们就可以看出，隐马尔可夫模型由两个部分组成：

• 状态 X 的底层马尔可夫链
• 每个 Time Step 的可观察到的输出 (Output)，也被称为状态的影响 (Effect of state) 或是辐射 (Emissions)

接下来看，隐马尔可夫模型需要哪些参数：

• 初始分布 (Initial Distribution) P(X₁)
• 转换概率 (Transition Probability) P(X_t|X_t-1)
• 辐射概率 (Emission Probability) P(E_t|X_t)

可以看到，比起马尔科夫链，仅多出了辐射概率，这可以理解为，当前状态下，出现这种（被我们观察到的）现象的概率。

独立性（Independencies）

相较于马尔科夫链，隐马尔科夫模型的独立性会更复杂一些，因为现在不仅有 States，还有 emissions. 我们先来看使用贝叶斯网络的链式法则来表示隐马尔科夫模型是怎样的：

有了上面这个式子后，我们可以结合一般链式法则来看具体的独立关系，依然用之前的天气例子，我们之前 X₁ ~ X₃：

我们对独立性进行一个总结：

• 在给定前一个状态 (Previous State) 的情况下，当前的状态 (Current State) 与之前所有的状态 (All past states) 和辐射 (All past emissions) 都相互独立。即 𝑋_t ⊥ 𝑋₁, … , 𝑋_t-2, 𝐸₁, … , 𝐸_t-2|𝑋_t-1
• 给定当前的状态 (Current State) ，当前的辐射 (Current Emission) 与之前所有的辐射 (All past emissions) 和状态 (All past states) 相互独立。即 𝐸_t ⊥ 𝑋₁, … , 𝑋_t-1, 𝐸₁, … , 𝐸_t-1|𝑋_t

和马尔科夫链一样，我们设想转换概率 (Transition Probability) 和辐射概率 (Emission Probability) 都不会改变。

推断（Inference）

在了解了隐马尔可夫模型中的独立性关系之后，我们可以用它们来进行推断。首先，我们来看如何跟踪随时间变化的分布 P(X_t)，这个任务也被称为 Filtering/Monitoring。我们用 B(X_t|e₁, …, e_t) 来表示状态的信念 (Belief of State)，我们从一个正态分布 (Uniform Distribution) B(X₁) 开始，随时间更新信念 B(X_t) 并且得到新的观察 (Observation)。 可以看到，这个推断主要分为两个部分：

• Passage of time（时间推移）
• Observation（观察结果）

Passage of time（时间推移）

这一部分的计算和马尔科夫链很相似，都是在更新状态。假设我们知道当前的信念 B(X_t) = P(X_t|e₁, …, e_t)，随着时间流逝，我们想要更新它，也就是得到 P(X_t+1|e₁, …, e_t)。

根据总概率定理（Law of total probability）以及贝叶斯公式，我们能够得到下列计算公式：

Observation（观察结果）

因为我们想要的是新的状态信念 B(X_t+1) = P(X_t+1|e₁, …, e_t+1)，在得到了 P(X_t+1|e₁, …, e_t) 之后，我们还需要新的观察结果。

在上述的计算过程中，可能有两个地方会让人觉得困惑，这实际上就是我在 笔记01 中介绍的计算多变量条件概率的一个 Trick，通过移动 ‘|’ 来进行计算。根据以上的计算，我们就能够得到更新后的信念应该为：

也就是辐射概率和时间推移的状态变换结果的乘积，但是要注意的的是，这里的结果并不是 ‘=’，而是正比于 ‘∝’，因此需要对结果进行标准化。

下面我们看一个具体的例子，还是用天气来作为例子，要注意的是，现在我们关注的是一个隐马尔可夫模型，因此，天气的状态 (State) 我们是无法直接观察到的，这里我们智能观察到别人是否带了雨伞。

首先，一开始我们不确定是晴天还是雨天，所以最开始的初始分布是一个标准正态分布 <0.5, 0.5>。在此基础上的时间推移 (Passage of time) 是不会有变化的，所以我们直接来看 X₂ 对应的 Observation。因为我们能观察到结果是“带了雨伞”，因此此时再计算天气状态的信念为：

B(Sun) = B(Sun) * P(Umb|Sun) = 0.5 * 0.1 = 0.05 / 0.45 = 0.182
B(Rain) = B(Rain) * P(Umb|Rain) = 0.5 * 0.8 = 0.4 / 0.45 = 0.818

接下来，我们需要再次进行时间推移，也就是计算在 X₃ 的状态信念：

B(Sun) = P(Sun|Sun) * B(Sun) + P(Sun|Rain) * B(Rain) = 0.7 * 0.182 + 0.3 * 0.818 = 0.373
B(Rain) = P(Rain|Sun) * B(Sun) + P(Rain|Rain) * B(Rain) = 0.3 * 0.182 + 0.7 * 0.818 = 0.627

后续的计算相同。

Forward Algorithm

上述的这部分推断过程，即“给定一个序列的观察结果，我们想得到该序列最终的状态信念”，就是 Forward 算法的定义。

下面给出 Forward 算法的伪代码：

该算法的复杂度为 O(n|X|²)

Example: Robot Localization（机器人定位）

接下来我们来看一个具体的例子 - Robot Localization（机器人定位）。我们在之前已经了解了 Forward 算法，这个算法就是给定所有过去的 Observation 来估计当前的状态信念 (Belief of state)。我们就使用这个算法去理解这个例子。

在本例中，我们使用一个机器人，这个机器人由两个部分组成：传感器模块 (Sensor Model) 和移动模块 (Motion Model)。后者负责进行移动，该模块有很小的概率会出现错误，比如输入了移动指令，但是没有执行；前者的传感器，用于测量它与附近障碍物之间的距离，我们可以看做就是用来检测机器人的四周是否有障碍物。我们将机器人活动的地图进行简化，可以用离散值表示每个位置，同时我们设定每次最多只会有一个方向的传感器产生错误。如此一来，我们能够得到如下的模型：

这是最初时刻的情况，在这个例子中，我们使用不同深浅的颜色来表示概率，所以在一开始，我们完全不知道机器人在哪里，因此所有位置的概率都相同，是一个标准正态分布。在这个例子中，机器人具体所在的位置就是隐藏状态 (Hidden State)，我们无法直接观测到，我们所能知道的，只有传感器的读数，因此传感器的数据就是观察结果 (Observation)。

在 t = 1 时刻，根据 First Observation（传感器的数据），我们了解到 Robot 上下都有墙壁，而左右都没有墙壁。在有这个观察结果的情况下，我们去更新 State 的概率，这时，有一些 State 会有更高的可能性，一些 State 可能性会比较低，还有一些或许根本不可能。我们的 Sensor 最多只会出现一个错误。如果要去到位置 1，这需要 2 个错误一起发生才行：一个是左侧没有障碍（但实际上有墙壁），二是下方有墙（但实际上没有）。因此位置 1 是不可能的。而在位置 6，如果下方的传感器出现错误，是有可能的，但因为概率较小，所以我们没有完全涂白。

在 t = 2 时刻，我们输入了“右移”的指令。根据传感器的数据，我们了解到 Robot 上下都有墙壁，而左右都没有墙壁。但是，移动模块可能会发生错误，所以仍有可能在原地不动，但是概率会比较小。也有其他的一些位置因为“右移”这个指令，相比上一个时刻可能性会变得更低。接下来，根据每个时刻的观察（传感器数据），我们可以不断更新对状态 (机器人位置) 的信念，越来越确信机器人的位置。

最可能的解释（Most Probable Explanation - MPE）

Forward 算法用于跟踪状态的概率 (Probability of states) 。这些概率随着时间的流逝和我们观察到的结果而更新。现在我们来看另一种任务，找到最可能的解释 (MPE)。意思是我们要考虑所有可能的状态组合 (State combinations) ，找到使得观察结果可能性最高的那一组。换言之，给定所有的 Emissions e₁, e₂, …, e_n，找到对于所有变量 x₁, x₂, …, x_n 最有可能的赋值。即计算：

我们可以看作有一个联合概率表，我们想要找到最大的
P(X₁, X₂, X₃|E₁,E₂,E₃) 以及使得这个概率最大的对于变量 X 的一组赋值。

Viterbi Algorithm

State Trellis

State Trellis 是用来表示状态随着时间转换的图。这个图中的每一个节点表示一个状态 (变量的取值)，每一条边表示可能的转换，每条边上权值表示为 P(x_t|x_t-1) P(e_t|x_t)

因此，图中的一条通路 (Path) 就表示一个状态序列 (Sequence of states)。这条通路上的权值乘积即为根据 Observation 得到的序列概率。Forward 算法就是把指向同一个 State 的路径权值相加，比如计算 P(X_n = Sun) = P(Sun|Sun) * P(e_n|Sun) + P(Sun|Rain) * P(e_n|Sun)。接下来我们要介绍 Viterbi 算法，用来计算具有最高概率的路径。

算法理解

接下来我们用一个具体的例子来理解 Viterbi 算法。假设我们有两个概率不等的硬币 c₁ 和 c₂。某个人不断掷硬币，但我们不清楚他每次掷的是哪一枚硬币，只能观察到每次掷出来的结果是 Head 还是 Tail。但是我们知道 c₁ 掷出 Head 的概率更高， c₂ 掷出 Tail 的概率要高。同样，这个掷硬币的人更喜欢掷相同的硬币（这意味着他每次换硬币的概率比较低），此人一开始随机选择一枚硬币开始实验。 因此，我们能得到下列参数：

接下来，我们来看具体如何计算：

首先，每一次我们都选择可能性最高的情况。

在 X₂ ，如果掷的硬币是 c₁, 它有两种可能：一开始就是 c₁ ，保持选择；一开始 c₂ ，改变为 c₁。我们选择可能性最高的一种：

MAX(keep：0.5 * 0.7，change：0.5 * 0.3)

显然，Keep 的可能性更高，之后，因为我们观察到的结果是 Head，所以要乘上对应的 Emission 概率：

0.35 * P(E₂ = Head|X₂ = c₁) = 0.35 * 0.8 = 0.28

如果在 X₂ 选择的是 c2，同样有两种情况：

MAX(Keep：0.5 * 0.75，Change：0.5 * 0.3)

依旧选择 Keep，之后，乘以 Emission 概率：

0.35 * P(E₂ = Head|X₂ = c₂) = 0.35 * 0.2 = 0.07

因此，在当前 Time Step 最有可能的结论就是一开始选 c₁，之后保持选择，依旧是 c₁，再掷一次，得到 Head。

后面的分析一样。在 X3,如果是 c1：

MAX(Keep: 0.28 * 0.7, Change: 0.07 * 0.3)

依旧 Keep，乘以 Emission 概率：0.196 * 0.2 = 0.0392

如果是 c2：

MAX( Keep: 0.07 * 0.7, Change: 0.28 * 0.3,)

选择 Change，乘以 Evidence：0.084 * 0.8 = 0.0672。因此最有可能的是改换 c₂ 硬币，投出 Tail。 所以如果我们的序列就到此为止，那么根据计算可以得到的 MPE 就是 <c₁, c₁, c₂>

而当我们完整计算之后，可以得到如下结果：

最终的答案应该为 <c₁, c₁, c₂, c₂, c₂>。下面给出 Viterbi 算法的伪代码：

该算法的复杂度为 O(n|X|²)。但需要注意的是，该算法最终输出的只是最可能状态序列的概率，也就是 max(P(x_1:t-1, x_t|e_1:t))。但我们真正感兴趣的是 argmax(P(x_1:t-1, x_t|e_1:t))，也就是具体的状态序列（每个 x_i 的取值）。对此，一般有两种比较常见的解决方式，一是使用一个额外的数据结构来保存路径，二是从最后一个节点进行回溯计算 (Backtracking Computation).

回溯计算（Backtracking Computation）

所谓的回溯计算，就是把我们在进行正向计算的步骤反转过来。我们之前的正向计算实际上就是一直在使用下面的公式：

因此，在进行回溯计算的时候，我们需要以下 3 个步骤：

1. 除以 Emission 概率
2. 对 x_t-1 的每一个状态 (State)，我们除以对应的转换概率 P(x_t|x_t-1)
3. 查看 x_t-1 的哪个取值/状态符合结果

我们使用之间计算的结果来进行上述的回溯计算：

首先，我们已经得到了 MPE 为 0.02107，因此，我们就从这里开始进行回溯：

Step 1. 除以 Emission 概率
已知 X₅ = c₂，观察结果为 Tail，所以 0.02107 / 0.8 = 0.0263375

Step 2. 对 x_t-1 的每一个状态 (State)，我们除以对应的转换概率 P(x_t|x_t-1)
x₄ = c₁ = 0.0263375 / 0.3 = 0.08779
x₄ = c₂ = 0.0263375 / 0.7 = 0.03763

Step 3. 查看 x_t-1 的哪个取值/状态符合结果
显然，x₄ = c₂ 是符合正向计算结果的

我们在正向计算的过程中，实际上可以发现，概率随着序列越来越小，这是因为我们一直在做乘法，这有可能会导致概率消失 (Probability Vanishing) 或者叫下溢 (Underflow) ，也就是没有充足的位数去表示一个特别小的小数，因此会直接变为 0。这和机械最小值 (Machine Epsilon) 密切相关，对32bits 的数字，其机械最小值为10^-7，对 64bits 的数字，则为10^-16。一个长序列（例如100个步骤）将很容易导致下溢，并且概率将变为 0。为了避免这个现象，可以使用对数概率 (Log Probability) （与 NBC 朴素贝叶斯分类器类似），也就是求取概率的对数，这样一来，概率之间的相乘也会变为相加，进一步降低计算压力。

我们感兴趣的范围 [0, 1] 内的数字会变为负数，但是增减性不会改变。我们可以把所有的概率都转换为 log 形式，这样在计算的过程中，不需要再保持书写 log 式。

粒子过滤（Particle Filtering）

根据这个算法的名字我们就能够知道，它和 Forward 的作用一样，同样是跟踪随时间变化的分布 P(X_t)，区别在于该算法给出的是一个近似解 (Approximate Solution)。之所以会有这个算法，是因为有时 |X|，变量的取值/状态 (State) 数量太多了或者状态取值时连续值 (Continuous Value)，不好进行准确推断。因此，我们使用近似推断 (Approximate Inference)，该方法的核心思想在于采样 (Sample)。因此：

• 跟踪 X 的样本，而不是所有的值
• 样本就被称为粒子 (Particles)
• 每个 Time Step 与样本数成线性关系
• 样本的数量需要足够多
• 在内存 (Memory) 中存储的是粒子的列表而不是 States 的列表

这里我们给出一个例子，图 1 表示我们一共有 9 个 States，每个 State 的出现概率如图所示。这里我们一共使用了 10 个 Particles，我们所需要保存的就是一系列 Particles 的一个列表，比如它们的坐标 [(1, 2), (2, 3), (2, 3), … ]。这样一来，就不需要再保存所有的状态对应的概率表以及映射图，那样太耗费空间。同时，这也保证了整个算法的时间复杂度与样本的数量呈线性关系，因为我们所有的操作只针对这个列表中的这些 Particles，与 |X| 不相关，一般来说 N << |X|。 这个方法就是机器人定位现实中采用的方案。

我们会使用值为 x 的粒子的数量来近似 P(x)。目前，我们把所有粒子的权值设为 1。我们要牢记，HMM 的改变由两个部分组成：Passage of time 以及 Observe。

我们先来看 Passage of time。

在原本的计算中我们需要计算 P(Xt+1|Xt)，现在我们只需要计算 P(Xt+1|Xt = (1, 3))，即Particle从 (1, 3) 位置移动到其他位置的概率，因为大多数粒子都是顺时针顺序移动的，比如这里我们可能会有
x2,3 = 0.5
x2,2 = 0.2
x3,3 = 0.2
x3,2 = 0.1

我们可以想象一条 Tape，被这些概率分为很多的区间，我们随机生成一个 [0,1] 的数，根据它落在哪个区间决定粒子会移动到哪个位置

在完成了转换之后，我们会得到新的观察结果。然后根据观察结果去改变粒子的权值。我们还用机器人定位的例子来理解，我们可以假想一个机器人，它具备一个感应器来判断它所处的位置，但这个感应器并不完美，因此，它只会对当前最有可能的位置给予最高的概率，但是其他的位置仍具备一些可能

这里会存在一个问题，随着Transition不断进行，权值可能会越来越小，最终接近0，也就是有 Weight Vanishing 的风险。所以这里提出重采样 (Resample) 来避免这个问题，每经过 N 个 Time Step，进行一次重采样。比如这里，我们对所有权值求和，得到 Z，然后所有权值都除以 Z，得到一个标准化概率，在这个概率的基础上，再次使用之前的 Tape 思想，分布节点，即为重采样的过程

所以整体的流程如下图所示：采样 - 变换 - 更新权值 - 重采样