使用 Perl 语言开始隐藏马尔可夫模型

马尔可夫模型（以数学家安德烈·马尔可夫命名）用于随机变化系统中的预测。马尔可夫的洞察力在于，在这种情况下，好的预测只能从事件最近一次的发生中得出，而忽略当前事件之前的任何发生。这种方法可以被描述为无记忆或与历史无关的预测。

马尔可夫的第一个例子（1913 年）预测了普希金诗歌《叶甫盖尼·奥涅金》中的元音出现。今天的挑战是找到一个马尔可夫模型不起作用的研究领域。这些模型用于研究热力学和统计力学；生物信息学、酶活性和种群动态；太阳辐射和风力；价格趋势；语音识别和生成；数据压缩和模式识别；强化学习和手势识别。该列表还在不断增加。

本文介绍了隐藏马尔可夫模型（HMM），它是原始模型的改进。本文首先介绍 HMM 和适用于 HMM 的 Viterbi 算法，然后提供一个完整的 Perl 代码示例。

1. 具有马尔可夫性质的随机游走

想象一下一只正在寻找食物的狐狸，目前位于 C 位置（例如，在溪流旁边的灌木丛旁）。狐狸搜索路径上的先前位置是 P1、P2、P3 等。现在的问题是如何预测狐狸的下一个位置。

一种方法是使用整个搜索历史 P1、P2、...、C 来预测下一个位置。例如，已经访问过的狐狸搜索位置可能会被赋予非常低的概率成为下一个位置，理由是狐狸足够聪明，不会重复失败的搜索位置。相比之下，附近但未访问过的位置（例如，溪流对面的露头）可能会被赋予很高的概率。

一种更简单的方法是假设可以仅从当前位置 C 就能够像从整个搜索历史 P1、P2、...、C 中一样准确地预测狐狸的下一个位置 N。如果是这样，狐狸的搜索可以被建模为一个具有马尔可夫性质的过程。在本例中，该性质意味着狐狸的下一个搜索位置的概率分布仅取决于狐狸的当前位置，而不取决于它之前的位置。因此，马尔可夫过程可以被描述为表现出马尔可夫性质的随机游走。

马尔可夫性质可以使用条件概率的符号简洁地表示：P(A|B) 是在给定 B 的情况下 A 的概率。在狐狸的例子中，马尔可夫性质支持一个等式：P(N|P1,P2,…,C) = P(N|C)，其中C和N分别是当前位置和下一个位置。

2. HMM 和 Viterbi 算法

在狐狸的例子中，过程中每个感兴趣的状态（在特定位置觅食）对于观察者来说是可见的：没有状态是隐藏的。在 HMM 改进中，状态是隐藏的，但每个状态都具有一个或多个相关的输出标记（又名观察结果），这些标记不是隐藏的。目标是从观察结果推断状态序列。一个经典的例子是 瓮问题。

想象一个隐藏的房间，里面装有多个瓮，每个瓮里都有标记为 B1、B2 等的球。不知道每个瓮是否都包含相同的球。一个机器人在秘密中选择一个瓮，然后从中随机选择一个球，并将球的标签副本放在观察者看到的传送带上。因此，观察者知道拣选球的序列（例如，B9、B4、B13、...），但不知道从哪个瓮中拣选给定的球。从瓮中拣选一个球并复制标签后，机器人在下次拣选之前将球返回到其所在的瓮。

在此过程中，拣选的瓮属于隐藏状态，放置在传送带上的球标签副本是输出标记。机器人有一个用于拣选每个瓮的配方（尽管是秘密的）；因此，这种拣选瓮可以被建模为具有隐藏状态的马尔可夫过程——HMM。

来自 20 世纪 60 年代后期的 Viterbi 算法可用于从观察到的事件（在本例中为传送带上的标签）推断隐藏状态的可能序列（在本例中为拣选瓮的序列）。推断的序列称为 Viterbi 路径，以纪念 Viterbi。

作为一个例子，考虑这个输出标记的切片

B9,B4,B13,B27,B6,B13,B29 ## sequence of 7 observations

底层 Viterbi 路径可能是

  B9     B4    B13     B27    B6     B13    B29  ## observations
  /      /      /      /      /      /      /
UrnD-->UrnA-->UrnF-->UrnA-->UrnG-->UrnF-->UrnN   ## hidden states

Viterbi 算法从传送带上观察到的标签推断机器人的隐藏状态转换（拣选瓮）。当然，推断本质上是统计性的。存在猜测，但 Viterbi 算法应该使这种猜测具有教育意义。

3. Viterbi 算法的详细信息

Viterbi 算法使用六个输入来生成 Viterbi 路径

• 观察空间由可能的输出标记或观察结果组成（例如，B1、B2、B3 等）。
• 排序的观察结果是观察到标记的顺序（例如，B9、B4、B13 等）。
• 状态空间由可能的隐藏状态组成（例如，UrnA、UrnB 等）。
• 初始概率是特定瓮代表机器人选择过程中的起始状态的可能性。
• 转移概率是从一个瓮到另一个瓮的转移（包括返回到自身）的可能性。
• 发射概率是观察到的球来自特定瓮的可能性。

输出是

• 推断的隐藏状态序列（在本例中为选定的瓮序列），即 Viterbi 路径。

这些输入可以从示例中学习，而不是任意规定。因此，以下简短的代码示例分为两部分。第一部分模拟从示例中进行的机器学习，第二部分使用此学习的结果来推断观察到的球标签序列背后的 HMM。

4. hmm 程序

下面的 hmm 程序可从 我的网站 获取。Unix 类型系统预装了 Perl，但是该程序需要 Algorithm::Viterbi 软件包，该软件包可在 CPAN 和其他 Perl 存储库中找到。可选的软件包 Data::Dumper 可用于完整地查看 Viterbi 输入；为此，请取消注释程序的最后一行。

示例 1. hmm 程序

use strict;
use warnings;
use Algorithm::Viterbi;
use Data::Dumper; ## to view the inputs in full

use constant HowMany    => 64_000; # 64K
use constant UrnCount   =>      5; # five urns
use constant TokenCount =>     64; # observations

my @urns   =  ('UrnA', 'UrnB', 'UrnC', 'UrnD', 'UrnE');

my %balls = ('UrnA' => ['B1', 'B2', 'B4', 'B5', 'B7', 'B8'],
             'UrnB' => ['B1', 'B3', 'B5', 'B7'],
             'UrnC' => ['B2', 'B3', 'B4', 'B6', 'B7'],
             'UrnD' => ['B1', 'B2', 'B3', 'B5', 'B7'],
             'UrnE' => ['B2', 'B3', 'B5', 'B8']);

my @data = (); # list of ball-urn pairs for training

sub get_ball_urn_pair { ## return a ball-urn pair or just a ball
    my $both = shift;
    my $urn = $urns[int(rand(UrnCount))];
    my $list = $balls{$urn};
    my $i = int(rand(scalar @$list));
    return $both ? ($list->[$i], $urn) : $list->[$i];
}

#### execute
## randomly pick an urn and a ball from it
srand(time());
for (my $i = 0; $i < HowMany; $i++) {
    my @add = get_ball_urn_pair(1); # random pick
    push(@data, \@add);
}

## train the algorithm
my $viterbi = Algorithm::Viterbi->new();
$viterbi->train(\@data);

## generate observations
my @tokens = (); ## observed labels
for (my $i = 0; $i < TokenCount; $i++) {
    my $ball = get_ball_urn_pair(0); # random pick
    push(@tokens, $ball);
}

## output Viterbi path
my $v_path = ($viterbi->forward_viterbi(\@tokens))[1]; ## Viterbi path
print join('->', @$v_path), "\n";

# print Dumper($viterbi);  ## uncomment to see the probabilities in full

hmm 程序使用五个瓮（UrnA 到 UrnE），但每个瓮中都装有不同的球

my %balls = ('UrnA' => ['B1', 'B2', 'B4', 'B5', 'B7', 'B8'],
             'UrnB' => ['B1', 'B3', 'B5', 'B7'],
             'UrnC' => ['B2', 'B3', 'B4', 'B6', 'B7'],
             'UrnD' => ['B1', 'B2', 'B3', 'B5', 'B7'],
             'UrnE' => ['B2', 'B3', 'B5', 'B8']);

这些球存储在一个映射中，其中瓮作为键，包含的球列表作为该键的值。

hmm 程序中的学习部分使用随机生成的测试数据，这些数据由 64,000 个球-瓮对组成，例如

B5, UrnD ## evidence

此对确认 UrnD 包含一个标记为 B5 的球，但其他对确认球 B5 可能来自其他瓮

B5, UrnA ## more evidence

这两个示例共同表明 B5 球可能来自 UrnA 或 UrnD。学习对没有提供关于机器人的任何信息，只是它可能从特定的瓮中拣选一个特定的球。

来自学习部分的测试数据存储在一个列表中，该列表传递给 Algorithm::Viterbi 实例的 train 方法

my $viterbi = Algorithm::Viterbi->new(); ## instantiate
$viterbi->train(\@data);                 ## pass training data

接下来将阐明对此 train 方法调用的结果。

5. 训练的结果

train 方法检测到底层 HMM 中的五个可能状态（瓮），因为每个瓮在随机生成的测试数据中至少出现一次。以下是按检测顺序排列的状态

states => [UrnA, UrnD, UrnE, UrnB, UrnC]

回想一下，Viterbi 算法所需的三种概率是初始概率（用于起始瓮）；转移概率（机器人从当前瓮到下一个瓮的移动）；和发射概率（观察到的标签标识来自特定瓮的球的可能性）。训练产生这些初始概率的结果

start => {UrnA => 0.199578125,
          UrnD => 0.201703125,
          UrnE => 0.198546875,
          UrnB => 0.200343750,
          UrnC => 0.199828125}

这些值总和为 1，使五个瓮中的每一个瓮大约有五分之一的机会成为机器人的初始选择。UrnD 比 UrnB 略有优势，但结果基本上是随机的。

以下是从测试中得出的转移概率，以 UrnC 为例

UrnC => {UrnA => 0.198230642762076,
         UrnE => 0.192885810970331,
         UrnB => 0.204414287942599,
         UrnD => 0.205128205128205,
         UrnC => 0.198373602314489}

以 UrnC 作为当前状态，转移到任何瓮（包括再次回到 UrnC）的概率约为 20%。

发射概率很有趣，因为并非每个瓮都有每个球。例如，B6 仅出现在 UrnC 中。但是，Viterbi 算法无法根据训练示例排除其他某个瓮也可能包含此球的可能性。训练的结果仅保守地给出 20% 的可能性，即发射的 B6 标签来自 UrnC 中的球

B6 => {UrnC => 0.201110329189147}

相比之下，B7 出现在四个瓮中，训练数据包括这些出现的情况，因此每个瓮都有 B7 源自其中的概率

B7 => {UrnA => 0.160259923275660,
       UrnB => 0.244657619716113,
       UrnD => 0.197923929041754,
       UrnC => 0.201266713581985}

从训练数据来看，UrnB 是 B7 球最有可能的来源，而 UrnA 是此球最不可能的来源。（UrnA 和 UrnB 都包含一个 B7 球。）

回顾一下，训练数据支持对三个关键事项进行统计推断

• 机器人可能最初选择的瓮。在这种情况下，每个瓮被选择的机会大致相同。
• 给定一个选择的瓮，机器人可能接下来选择的瓮，它可能是同一个瓮。训练数据表明，所有瓮到瓮的转移都同样可能。
• 已从中选择特定球（例如，B7）的瓮。这里的可能性差异很大，因为并非每个瓮都包含相同的球。

Viterbi 算法使用三个输入概率，尤其是最后一个关于发射的概率，试图从观察到的标签序列中确定从中拣选带有观察到的标签的球的瓮的序列。

6. 观察结果

Viterbi 算法的最终输入由观察结果组成——传送带上可见的标签。hmm 程序也随机生成这些值

my @tokens = (); ## observed labels
for (my $i = 0; $i < TokenCount; $i++) {
    my $ball = get_ball_urn_pair(0); ## random pick
    push(@tokens, $ball);
}

有了观察结果，最后一步是将它们作为参数传递给 forward_viterbi 方法，该方法推断底层 HMM 并将其作为 Viterbi 路径返回

my $v_path = ($viterbi->forward_viterbi(\@tokens))[1]; ## Viterbi path

然后，程序将结果输出为状态转移序列

UrnB->UrnE->UrnE->UrnB->UrnC->...->UrnB->UrnE->UrnB->UrnA

在此示例运行中，只有 UrnD 未出现在 Viterbi 路径中。UrnD 中的每个球都至少出现在另一个瓮中，并且转移概率不支持 UrnD 作为目标。因此，该算法可以呈现一个合理的状态序列，该序列将 UrnD 从 HMM 中排除。

forward_viterbi 方法还会返回两个概率，它们表示对推断的 Viterbi 路径的置信度。为简单起见，这里仅输出 Viterbi 路径。

7. 分而治之与动态编程

hmm 程序运行速度很快，这引出了 Viterbi 算法如何工作的问题。该算法构建一个网格（trellis），这是一个图，其节点位于垂直切片中：较高的切片表示较早的时间，较低的切片表示较晚的时间。图 1 是一个来自 Wikipedia 的通用示例，其中最短的从左到右的路径以红色显示。

图 1. 一个网格（Trellis）

​

Viterbi 算法从顶部切片中的一个节点开始，计算通过网格的最短路径。 这就是 Viterbi 路径，其节点表示从观测中推断出的 HMM 状态。 该算法使用经典的动态规划技术来提高效率。 一个更简单（也可能更熟悉）的例子可以阐明这项技术，突出显示 Perl 中非常容易实现的“记忆化”。

考虑一种分而治之的方法来计算斐波那契数，使用递归 fib 函数作为实现。

sub fib { ## divide and conquer: exponential time
   my ($n) = shift;                  # read argument
   return $n if $n < 2;              # base case
   return fib($n - 1) + fib($n - 2); # recursive case
}

为了计算诸如 fib(40) 之类的值，该函数将问题分解为 fib(39) 和 fib(38) 的子问题。 这些划分将继续，直到子问题足够小以直接解决为止； 在这种情况下，计算 fib(0) 或 fib(1) 的值。

fib(40) 的值为 102,334,155，但是递归计算需要调用 fib 函数 331,160,281 次，这突出了这种分而治之设计的指数时间复杂度。 困难在于，由于使用 n-1 和 n-2 作为参数的递归调用，太多的值被一遍又一遍地重新计算。 例如，fib(40) 调用 fib(39) 和 fib(38)； 但是 fib(39) 再次调用 fib(38)。 后两个调用中的每一个都会导致调用 fib(37)，依此类推。

需要改变设计。 可以记忆化已编写的 fib 函数，从而将设计从分而治之更改为动态编程。 两种设计都将问题分解为更易于管理的子问题，直到可以简单地解决子问题为止。 这两种方法的不同之处在于，动态编程缓存子问题解决方案，以便可以在需要时重复使用这些解决方案。 简而言之，缓存代替了重新计算。 结果可以显着提高效率。

在 Perl 中很容易获得 fib 的记忆化版本，而无需更改原始源代码。

memoize('fib'); # dynamic-programming: linear time
fib(40);        # 41 calls

使用记忆化版本，fib(40) 仅需要 41 个调用即可计算该值。 在幕后，memoize 调用为 fib 函数提供了一个额外的参数，该参数是一个哈希表，用于存储已经计算出的斐波那契值； 因此，每个所需的斐波那契值仅被计算一次。 分而治之版本的指数时间复杂度降低为线性时间复杂度。

Viterbi 算法在使用相同的缓存方法计算通过网格的最短路径，该路径表示 HMM。 hmm 程序之所以高效，是因为 Viterbi 算法非常高效。