隐马尔可夫模型

隐马尔可夫模型（hidden Markov model，HMM）用一个不可直接观测的马尔可夫状态序列，产生随时间出现的观测。状态如何转移与状态如何发出观测是两套概率；看到同一个观测不一定意味着处于同一个状态。

设一台合成装置有状态0和1，传感器每步只报告符号0或1。初始状态概率为 $\pi =(0.6,0.4)$，转移矩阵和发射概率矩阵分别为 $$A=\left(\begin{array}{cc}0.8& 0.2\\ 0.3& 0.7\end{array}\right),B=\left(\begin{array}{cc}0.8& 0.2\\ 0.2& 0.8\end{array}\right).$$ A的第i行第j列表示从状态i转到状态j的概率；B的第j行第a列表示状态j发出符号a的概率。本例观察到序列 $(0,1,0)$，目标是计算整个序列的概率，并逐步更新当前状态的判断。

隐状态与条件独立

令隐状态 $Z_t\in \{0,\dots ,m-1\}$、观测 $Y_t\in \{0,\dots ,s-1\}$，t=1到T。π非负且总和一，A与B各行非负且各行和一。模型假定 $Pr(Z_{t+1}=j\mid Z_1,\dots ,Z_t)=a_{Z_{t}j}$；并且给定整个隐状态序列后，观测相互独立，$Y_t$ 的条件分布只依赖 $Z_t$。因此联合质量分解为 $$Pr(z_{1:T},y_{1:T})={\pi }_{z_1}{b}_{z_1}(y_1)\prod _{t=2}^T{a}_{z_{t-1}{z}_t}{b}_{z_t}(y_t).$$ 对观测和状态逐步求和，每一步对应一个概率行的和一，证明这定义了总质量为一的联合分布。状态的一阶记忆与观测的条件独立都是假设；观测序列本身一般既不独立，也不必是一阶马尔可夫序列。

箭头表示模型的条件依赖；每个观测由同一列的隐状态发出。没有直接观测到隐状态。

前向递推及其正确性

定义 ${\alpha }_t(j)=Pr(Y_{1:t}=y_{1:t},Z_t=j)$。初始值为 ${\alpha }_1(j)={\pi }_j{b}_j(y_1)$。将t+1时刻的联合事件按此前状态i划分，利用上述条件独立性，得到 $${\alpha }_{t+1}(j)=b_j(y_{t+1})\sum _{i=0}^{m-1}{\alpha }_t(i)a_{ij},Pr(Y_{1:T}=y_{1:T})=\sum _j{\alpha }_T(j).$$ 每条到达j的隐路径都有唯一的前一状态；递推把这些互斥路径相加，并没有选出其中最大的一条。这说明每一层保持α的联合概率含义，也证明算法正确。直接枚举有 $m^T$ 条隐路径，递推每步只做m乘m次合并。

主例第一步为 ${\alpha }_1=(0.6\cdot 0.8,0.4\cdot 0.2)=(0.48,0.08)$。第二个观测是1，所以 $${\alpha }_2(0)=(0.48\cdot 0.8+0.08\cdot 0.3)\cdot 0.2=0.0816,$$ $${\alpha }_2(1)=(0.48\cdot 0.2+0.08\cdot 0.7)\cdot 0.8=\mathrm{0.1216.}$$ 第三个观测为0，类似得到 ${\alpha }_3=(0.081408,0.020288)$，最终序列概率为 $0.101696$。在t时刻除以当层概率和，得到仅使用截至当时观测的滤波概率： $$Pr(Z_t=0\mid y_{1:t})\approx 0.857142857,0.401574803,0.800503461$$ 依次对应t=1、2、3。第一次看到0提高状态0的概率，下一次看到1又让状态1更可能；转移概率保留了此前状态信息。

点表示三个离散时刻的滤波概率，连线只帮助追踪变化；每个时刻两状态概率和为一。

归一化计算与可运行程序

长序列的联合概率连乘可能小到浮点下溢。可令 $f_t(j)=Pr(Z_t=j\mid y_{1:t})$，先从 $f_{t-1}$ 预测，再乘本次发射概率，最后归一化： $${\tilde{f}}_t(j)=b_j(y_t)\sum _i{f}_{t-1}(i)a_{ij},c_t=\sum _j{\tilde{f}}_t(j),f_t(j)={\tilde{f}}_t(j)/c_t.$$ 第一步的预测改用π。只要此前观测有正概率，$c_t=Pr(y_t\mid y_{1:t-1})$，所以 $\log Pr(y_{1:T})=\sum _t\log c_t$。本例三个缩放因子为0.56、约0.362857143、约0.500472441，它们的乘积仍是0.101696。

Python 3

import math

def forward_filter(initial, transition, emission, observations):
    initial = tuple(initial)
    transition = tuple(tuple(row) for row in transition)
    emission = tuple(tuple(row) for row in emission)
    observations = tuple(observations)
    states = len(initial)
    if states == 0 or len(transition) != states or len(emission) != states:
        raise ValueError("inconsistent state dimensions")
    symbols = len(emission[0])
    if symbols == 0 or any(len(row) != symbols for row in emission):
        raise ValueError("inconsistent emission dimensions")
    if any(len(row) != states for row in transition):
        raise ValueError("transition must be square")
    for row in (initial,) + transition + emission:
        if (any(not math.isfinite(p) or not 0 <= p <= 1 for p in row)
                or abs(sum(row) - 1) > 1e-12):
            raise ValueError("probability rows must sum to one")
    if any(type(o) is not int or not 0 <= o < symbols for o in observations):
        raise ValueError("invalid observation symbol")
    def multiply(a, b):
        value = a * b
        if a > 0 and b > 0 and value == 0:
            raise ArithmeticError("positive probability product underflows")
        return value

    filtered, history, log_likelihood = list(initial), [], 0.0
    for t, symbol in enumerate(observations):
        predicted = (filtered if t == 0 else
                     [sum(multiply(filtered[i], transition[i][j]) for i in range(states))
                      for j in range(states)])
        weights = [multiply(predicted[j], emission[j][symbol]) for j in range(states)]
        scale = sum(weights)
        if scale == 0:
            return -math.inf, None  # conditioning on an impossible prefix
        filtered = [w / scale for w in weights]
        if any(w > 0 and f == 0 for w, f in zip(weights, filtered)):
            raise ArithmeticError("positive filtering probability underflows")
        history.append(tuple(filtered))
        log_likelihood += math.log(scale)
    return log_likelihood, history

initial = [0.6, 0.4]
A = [[0.8, 0.2], [0.3, 0.7]]
B = [[0.8, 0.2], [0.2, 0.8]]
logp, belief = forward_filter(initial, A, B, [0, 1, 0])
assert abs(math.exp(logp) - 0.101696) < 1e-12
assert abs(belief[-1][0] - 0.8005034612964129) < 1e-12
assert forward_filter(initial, A, B, [])[0] == 0.0
assert forward_filter([1.0], [[1.0]], [[1.0, 0.0]], [1]) == (-math.inf, None)
print(math.exp(logp), belief[-1])

输入必须是有限的概率行与合法整数符号，代码不改写输入，返回对数似然及每一步滤波概率。空观测序列的概率取一，故对数为零。若模型使某个前缀真正不可能，返回负无穷及 None：这时条件状态分布没有定义，不能用均匀向量悄悄替代。缩放减少了长连乘下溢；程序还逐次检查正因子的乘积和归一化比值，若它们下溢成零，抛出 ArithmeticError，避免把数值范围不足误报为概率为零。遇到这种情况需改用对数域算法或更高精度。

给定已经存储的模型，递推时间为 $O(T{m}^2)$、单步工作空间为 $O(m)$；这份程序保留所有滤波向量，输出与历史另外占 $O(Tm)$，复制和验证模型还需 $O(m^2+ms)$。这不是训练程序，π、A、B都由调用者给定。

次序与不同推断目标

变式。 把相同的两个0和一个1改排为 $(1,0,0)$，前向三层依次为 $(0.12,0.32)$、$(0.1536,0.0496)$、$(0.110208,0.013088)$，序列概率变为0.123296。总计数相同、顺序不同，似然可以不同；无条件把观测当独立同分布会丢掉这种结构。

滤波只使用过去和当前观测；利用完整序列判断较早状态叫平滑，需要额外的后向信息。寻找后验概率最大的整条隐路径又是另一个目标，通常使用Viterbi算法中的最大化递推。前向算法将路径概率相加，不能把其和当作最大路径概率。即使逐时刻最可能的状态各自容易求，它们拼在一起也不保证是最可能的合法路径。

在实际解释中，隐状态名称来自建模约定。将状态0、1同时交换，并相应交换π、A、B的行列，观测分布不变；若两个状态的发射与转移行为无法区分，更多数据也不自动让它们获得唯一物理含义。