HMM经典介绍论文【Rabiner 1989】翻译（十六）——放大

来源：互联网发布：java temp dir 编辑：程序博客网时间：2024/05/16 19:35

5 HMM的实现问题

前面两节的讨论主要是关于HMM的理论以及模型的变体。这一节我们会讨论HMM的实现问题，包括放大、多观测序列、初始参数估计、数据丢失、模型大小以及类型的选择。对其中一些实现问题，我们可得到精确解析解；而对于其他问题，我们只能给出一些经验建议。

为了理解在HMM参数估计过程中为什么需要放大，考虑(18)中定义的αt(i)。可以看到αt(i)包含很多项的和，每一项的形式为

(\prod s = 1 t - 1 a q s q s + 1 \prod s = 1 t b q s (O s)),

其中qt=Si。由于每个a和b项是小于1的（一般是远小于1的），所以当t变得比较大时（比如10或者更多），αt(i)的每一项开始趋向于0。对非常大的t（100或者更多），αt(i)计算的动态范围会超出机器的精度范围（甚至是双精度）。所以，唯一有效的计算方式是结合一个放大过程。

基本的放大步骤是把αt(i)乘以一个放大系数并且这个放大系数是独立于i的（即只依赖于t），使得放大后的αt(i)（1≤t≤T）位于计算机的动态范围中。对βt(i)也要进行一个相似的放大步骤，在计算的最后，需要消掉放大系数。

为了更好地理解放大过程，考虑状态转移系数aij的估计公式。如果我们把估计公式(41)直接写成前向变量和后项变量的形式。可以得到

a ¯ i j = \sum T - 1 t = 1 α t ( i ) a i j b j ( O t + 1 ) β t + 1 ( j ) \sum T t = 1 \sum N j = 1 α t ( i ) a i j b j ( O t + 1 ) β t + 1 ( j ) . (90)

（原文中的符号表示有些混乱，下面做了些修改。）

首先令α¯1(i)=α1(i)。

对每个t，我们首先根据递推公式(20)计算α¯t(i)，然后乘以一个放大系数ct：

c t = 1 \sum N i = 1 α ¯ t ( i ) . (91)

于是，对固定的t，我们首先计算

α ¯ t (i) = \sum j = 1 N α^t - 1 (j) a j i b i (O t) . (92a)

然后计算放大α^t(i)：

α^t (i) = α ¯ i \sum N i = 1 α ¯ i = \sum N j = 1 α ^ t - 1 ( j ) a j i b i ( O t ) \sum N i = 1 \sum N j = 1 α ^ t - 1 ( j ) a j i b i ( O t ) . (92b)

通过递推，可以把α^t−1(j)写成

α^t - 1 (j) = [\prod τ = 1 t - 1 c τ] α t - 1 (j) . (93a)

于是有

α^t (i) = \sum N j = 1 α t - 1 ( j ) ( \prod t - 1 τ = 1 c τ ) a i j b j ( O t ) \sum N i = 1 \sum N j = 1 α t - 1 ( j ) ( \prod t - 1 τ = 1 c τ ) a i j b j ( O t ) = α t ( i ) \sum N i = 1 α t ( i ) . (93b)

下一步我们用后向地推计算βt(i)。这里我们在每个t上仍然使用计算α时用的放大因子。于是放大后的β是

β^t (i) = c t β ¯ t (i) . (94)

因为每个放大因子把α放大到和为1，由于α和β项的值是相近的，使用同样的放大因子可以把β放大到有效范围。而且，参数估计公式(90)现在变成

a ¯ i j = \sum T - 1 t = 1 α ^ t ( i ) a i j b j ( O t + 1 ) β ^ t + 1 ( j ) \sum T - 1 t = 1 \sum N j = 1 α ^ t ( i ) a i j b j ( O t + 1 ) β ^ t + 1 ( j ), (95)

而每个α^t(i)可以写成

α^t (i) = [\prod s = 1 t c s] α t (i) = C t α t (i) (96)

每个β^t+1(j)可以写成

β^t + 1 (j) = [\prod s = t + 1 T c s] β t + 1 (j) = D t + 1 β t + 1 (j) . (97)

于是(95)可以写成

a ¯ i j = \sum T - 1 t = 1 C t α t ( i ) a i j b j ( O t + 1 ) D t + 1 β t + 1 ( j ) \sum T - 1 t = 1 \sum N j = 1 C t α t ( i ) a i j b j ( O t + 1 ) D t + 1 β t + 1 ( j ), (98)

而CtDt+1项可以写作

C t D t + 1 = \prod s = 1 t c s \prod s = t + 1 T c s = \prod s = 1 T c s = C T (99)

和t无关。于是(98)分子分母中的CtDt+1抵消了，得到的结果就是(90)对应的准确结果。（虽然是用放大后的α^进行计算的，但得到的结果竟然和真实值一样！）

显然，上面的放大过程同样适用于π和B的估计。缩放过程(92)不需要在每个时刻t都要执行，而是可以在需要或必要（比如防止下溢）的时候才执行。如果某时刻t没有进行放大，只要把放大系数ct设置成1，然后上面讨论的所有条件都满足。

由放大导致的唯一变化是P(O|λ)的计算过程。我们不能只对α^T(i)进行求和，因为这些值是放大后的值。但是我们可以利用性质

\prod t = 1 T c t \sum i = 1 N α T (i) = C T \sum i = 1 N α T (i) = 1. (100)

于是有

\prod t = 1 T c t \cdot P (O | λ) = 1 (101)

P (O | λ) = 1 \prod T t = 1 c t (102)

log [P (O | λ)] = - \sum t = 1 T log c t . (103)

于是得到P的log值，但是不是P，这是因为后者可能会超出计算机的动态范围。

最后我们注意到当使用Viterbi算法来得到似然最大的状态序列时，如果我们按下面的算法计算那么不需要放大过程。定义

ϕ t (i) = max q 1, q 2, \dots, q t log P [q 1 q 2 \dots q t, O 1 O 2 \dots O t | λ], (104)

初始值

ϕ 1 (i) = log (π i) + log [b i (O 1)], (105a)

递推步骤

ϕ t (j) = max 1 \leq I l e q N [ϕ t - 1 (i) + log a i j] + log [b j (O t)] (105b)

以及终止步骤

log P * = max 1 \leq i \leq N [ϕ T (i)] . (105c)

这样我们得到logP∗而非P∗，而计算量大大减少了并且没有数值问题。读者应该注意到(105b)中的logaij项可以事先计算。并且，当观测值有限时log[bj(Ot)]项也可以事先计算。

阅读全文

0 0