维特比算法

维特比算法（Viterbi algorithm）是一种动态规划算法。它用于寻找最有可能产生观测事件序列的维特比路径——隐含状态序列，特别是在马尔可夫信息源上下文和隐马尔可夫模型中。

术语“维特比路径”和“维特比算法”也被用于寻找观察结果最有可能解释相关的动态规划算法。例如在统计句法分析中动态规划算法可以被用于发现最可能的上下文无关的派生（解析）的字符串，有时被称为“维特比分析”。

维特比算法由安德鲁·维特比于1967年提出，用于在数字通信链路中解卷积以消除噪音。此算法被广泛应用于CDMA和GSM数字蜂窝网络、拨号调制解调器、卫星、深空通信和802.11无线网络中解卷积码。现今也被常常用于语音识别、关键字识别、计算语言学和生物信息学中。例如在语音（语音识别）中，声音信号做为观察到的事件序列，而文本字符串，被看作是隐含的产生声音信号的原因，因此可对声音信号应用维特比算法寻找最有可能的文本字符串。

## 算法

假设给定隐式马尔可夫模型（HMM）状态空间 $S$，共有k个状态，初始状态 $i$ 的概率为 $\pi_i$，从状态 $i$ 到状态 $j$ 的转移概率为 $a_{i,j}$。令观察到的输出为 $y_1,\dots, y_T$。产生观察结果的最有可能的状态序列 $x_1,\dots,x_T$ 由递推关系给出：

$$\begin{array}{rcl} V_{1,k} &=& \mathrm{P}\big( y_1 \ | \ k \big) \cdot \pi_k \\ V_{t,k} &=& \max_{x \in S} \left( \mathrm{P}\big( y_t \ | \ k \big) \cdot a_{x,k} \cdot V_{t-1,x}\right) \end{array}$$

此处 $V_{t,k}$ 是前 $t$ 个最终状态为 $k$ 的观测结果最有可能对应的状态序列的概率。通过保存向后指针记住在第二个等式中用到的状态 $x$ 可以获得维特比路径。声明一个函数 $\mathrm{Ptr}(k,t)$，它返回若 $t > 1$时计算 $V_{t,k}$ 用到的 $x$ 值 或若 $t=1$时的$k$ . 这样:

$$\begin{array}{rcl} x_T &=& \arg\max_{x \in S} (V_{T,x}) \\ x_{t-1} &=& \mathrm{Ptr}(x_t,t) \end{array}$$

这里我们使用了arg max的标准定义算法复杂度为 $O(T\times\left|{S}\right|^2)$

## 例子

想象一个乡村诊所。村民有着非常理想化的特性，要么健康要么发烧。他们只有问诊所的医生的才能知道是否发烧。聪明的医生通过询问病人的感觉诊断他们是否发烧。村民只回答他们感觉正常、头晕或冷。

假设一个病人每天来到诊所并告诉医生他的感觉。医生相信病人的健康状况如同一个离散马尔可夫链。病人的状态有两种“健康”和“发烧”，但医生不能直接观察到，这意味着状态对他是“隐含”的。每天病人会告诉医生自己有以下几种由他的健康状态决定的感觉的一种：正常、冷或头晕。这些是观察结果。整个系统为一个隐马尔可夫模型(HMM)。

医生知道村民的总体健康状况，还知道发烧和没发烧的病人通常会抱怨什么症状。换句话说，医生知道隐马尔可夫模型的参数。这可以用Python语言表示如下:

```python
states = ('Healthy', 'Fever')
 
observations = ('normal', 'cold', 'dizzy')
 
start_probability = {'Healthy': 0.6, 'Fever': 0.4}
 
transition_probability = {
   'Healthy' : {'Healthy': 0.7, 'Fever': 0.3},
   'Fever' : {'Healthy': 0.4, 'Fever': 0.6},
   }
 
emission_probability = {
   'Healthy' : {'normal': 0.5, 'cold': 0.4, 'dizzy': 0.1},
   'Fever' : {'normal': 0.1, 'cold': 0.3, 'dizzy': 0.6},
   }
```

在这段代码中, 起始概率`start_probability` 表示病人第一次到访时医生认为其所处的HMM状态，他唯一知道的是病人倾向于是健康的。这里用到的特定概率分布不是均衡的,如转移概率大约是`{'Healthy': 0.57, 'Fever': 0.43}`。转移概率`transition_probability`表示潜在的马尔可夫链中健康状态的变化。在这个例子中，当天健康的病人仅有30%的机会第二天会发烧。放射概率`emission_probability`表示每天病人感觉的可能性。假如他是健康的，50%会感觉正常。如果他发烧了，有60%的可能感觉到头晕。

<image>

病人连续三天看医生，医生发现第一天他感觉正常，第二天感觉冷，第三天感觉头晕。于是医生产生了一个问题：怎样的健康状态序列最能够解释这些观察结果。维特比算法解答了这个问题。

```python
def viterbi(obs, states, start_p, trans_p, emit_p):
    V = [{}]
    path = {}

    # Initialize
    for st in states:
        V[0][st] = start_p[st] * emit_p[st][obs[0]]
        path[st] = [st]

    # Run Viterbi when t > 0
    for t in range(1,len(obs)):
        V.append({})
        newpath = {}

        for curr_st in states:
            paths_to_curr_st = []
            for prev_st in states:
                paths_to_curr_st.append((V[t-1][prev_st] * trans_p[prev_st][curr_st] * emit_p[curr_st][obs[t]], prev_st))
            curr_prob, prev_state = max(paths_to_curr_st)
            V[t][curr_st] = curr_prob
            newpath[curr_st] = path[prev_state] + [curr_st]

        # No need to keep the old paths
        path = newpath

    for line in dptable(V):
        print(line)
    prob, end_state = max([(V[-1][st], st) for st in states])
    return prob, path[end_state]

def dptable(V):
    # Print a table of steps from dictionary
    yield ' ' * 4 + '    '.join(states)
    for t in range(len(V)):
        yield '{}   '.format(t) + '    '.join(['{:.4f}'.format(V[t][state]) for state in V[0]])
```

函数`viterbi` 具有以下参数: `obs` 为观察结果序列, 例如 `['normal', 'cold', 'dizzy']`； `states` 为一组隐含状态； `start_p` 为起始状态概率; `trans_p` 为转移概率; 而 `emit_p` 为放射概率。为了简化代码，我们假设观察序列 `obs` 非空且 `trans_p[i][j]` 和 `emit_p[i][j]` 对所有状态 i,j 有定义。

在运行的例子中正向/维特比算法使用如下:

```python
def example():
    return viterbi(observations,
                   states,
                   start_probability,
                   transition_probability,
                   emission_probability)
print(example())
```

维特比算法揭示了观察结果 `['normal', 'cold', 'dizzy']` 最有可能由状态序列 `['Healthy', 'Healthy', 'Fever']`产生。换句话说，对于观察到的活动, 病人第一天感到正常，第二天感到冷时都是健康的，而第三天发烧了。

维特比算法的计算过程可以直观地由格图表示。维特比路径本质上是穿过格式结构的最长路径。诊所例子的格式结构如下, 黑色加粗的是维特比路径：

<image>

在实现维特比算法时需注意许多编程语言使用浮点数计算，当 p 很小时可能会导致结果算术下溢。避免这一问题的常用技巧是在整个计算过程中使用对数概率，在对数系统中也使用了同样的技巧。当算法结束时，可以通过适当的幂运算获得精确结果。

## 伪代码

首先是一些问题必要的设置。设观察值空间为 $O=\{o_1,o_2,\dots,o_N\}$、状态空间为 $S=\{s_1,s_2,\dots,s_K\}$、观察值序列为 $Y=\{y_1,y_2,\ldots, y_T\}$, $A$ 为 $K\times K$转移矩阵，其中 $A_{ij}$ 为从状态 $s_i$ 转移到 $s_j$ 的转移概率、 $B$ 为 $K\times N$ 放射矩阵(emission matrix)，其中 $B_{ij}$ 为在状态 $s_i$ 观察到 $o_j$ 的概率、大小为 $K$ 的初始概率数组 $\pi$， $\pi_i$ 为 $x_1 == s_i$ 的概率。我们称路径 $X=\{x_1,x_2,\ldots,x_T\}$ 为生成观察值 $Y=\{y_1,y_2,\ldots, y_T\}$ 的状态序列。

在这个动态规划问题中, 我们构造了两个大小为 $K\times T$ 的二维表 $T_1, T_2$。$T_1$ 的每个元素, $T_1[i,j]$，保存生成 $Y=\{y_1,y_2,\ldots, y_j\}$ 时最有可能的路径 $\hat{X}=\{\hat{x}_1,\hat{x}_2,\ldots,\hat{x}_j\}$， $\hat{x}_j=s_i$ 的概率。$T_2$ 的每个元素, $T_2[i,j]$，保存最有可能路径 $\hat{X}=\{\hat{x}_1,\hat{x}_2,\ldots,\hat{x}_{j-1},\hat{x}_j\}$ , $\forall j, 2\leq j \leq T$ 的 $\hat{x}_{j-1}$。

我们按 $K\cdot j+i$ 增序填充两个表 $T_1[i,j],T_2[i,j]$。

$T_1[i,j]=\max_{k}{(T_1[k,j-1]\cdot A_{ki}\cdot B_{iy_j})}$,

$$T_2[i,j]=\arg\max_{k}{(T_1[k,j-1]\cdot A_{ki}\cdot B_{iy_j})}$$

输入:

- 观察空间 $O=\{o_1,o_2,\dots,o_N\}$。
- 状态 $S=\{s_1,s_2,\dots,s_K\}$。
- 观察序列 $Y=\{y_1,y_2,\ldots, y_T\}$ 若在$t$ 时间观察值为 $o_i$，则 $y_t==i$,
- 大小为 $K\cdot K$ 的转移矩阵 $A$，$A_{ij}$ 为从状态 $s_i$ 到 $s_j$ 的转移概率。
- 大小为 $K\cdot N$ 的放射矩阵 $B$，$B_{ij}$ 为状态 $s_i$ 观察到 $o_j$ 的概率。
- 大小为 $K$ 的初始概率数组 $\pi$，$\pi_i$ 为 $x_1 == s_i$ 的概率。

输出:

- 最有可能的隐含状态序列 $X=\{x_1,x_2,\ldots,x_T\}$

```
 '''function''' ''VITERBI''( ''O'', ''S'', ''π'', ''Y'', ''A'', ''B'' ) : X
     '''for''' each state ''s<sub>i</sub>'' '''do'''

     '''end for'''
     '''for''' ''i'' ← ''2'',''3'',...,''T'' '''do'''
         '''for''' each state ''s<sub>j</sub>'' '''do'''
             <math>T_1[j,i] \gets \max_{k}{(T_1[k,i-1]\cdot A_{kj}\cdot B_{jy_i})} </math>
             <math>T_2[j,i] \gets \arg\max_{k}{(T_1[k,i-1]\cdot A_{kj}\cdot B_{jy_i})} </math>
         '''end for'''
     '''end for'''
     <math>z_T \gets \arg\max_{k}{(T_1[k,T])} </math>

     '''for''' ''i'' ← ''T'',''T-1'',...,''2'' '''do'''

     '''end for'''
     '''return''' X
 '''end function'''
```

## 使用动态规划的算法

- 最长公共子序列
- Floyd-Warshall算法

## 参考资料

- (note: the Viterbi decoding algorithm is described in section IV.) Subscription required.