半隐马尔可夫模型（semi-hidden Markov model）详解

目录

1. 引言
2. 为什么需要半隐马尔可夫模型
3. 半隐马尔可夫模型的形式化定义
   1. 状态与持续时间
   2. 观测生成模型
4. 半隐马尔可夫模型的联合概率
5. 半隐马尔可夫模型的前向-后向算法
   1. 前向变量及其公式推导
   2. 后向变量及其公式推导
   3. 观测序列概率的计算
6. 持续时间分布与参数学习
7. 小结与应用场景
8. 代码示例

---

引言

在标准的隐马尔可夫模型（HMM）里，状态序列 { X t } \{X_t\} {Xt​} 每一步都可能发生转移或保持在同一状态。这种“自环”带来的状态持续时间服从几何分布，即停留 $d$ 步的概率为
$$(a_{ii}{)}^{d-1}(1-a_{ii}).$$
在很多现实场景下，这种几何分布的限制较强，我们需要对状态持续时间（duration）有更加灵活的建模。半隐马尔可夫模型（Semi-HMM），又称显式持续时间HMM（Explicit Duration HMM），就是为此而设计。

---

为什么需要半隐马尔可夫模型

• 标准HMM的局限：
  对于离散时间步长，每一步都存在转移概率 $a_{ij}$。若自环概率 $a_{ii}$ 较大，虽然可以实现“较长”地停留在状态 $i$，但其分布必然是几何型，未必与现实数据相符。

• 半隐马尔可夫模型的改进：
  引入显式持续时间分布 $p_i(d)$，表示在状态 $i$ 上一次性停留 $d$ 步的概率。这样就能灵活地对“停留时长”进行刻画，比如可以使用离散分布、高斯分布、Gamma分布等，更贴合具体领域（例如语音识别、基因片段识别等）。

---

半隐马尔可夫模型的形式化定义

状态与持续时间

我们仍然有隐藏状态集合 { s 1 , s 2 , … , s N } \{s_1, s_2, \dots, s_N\} {s1​,s2​,…,sN​}。但每个状态 $s_i$ 会带有一个显式的持续时间分布：
$$p_i(d)=P(\text{停留时长}=d|X_t=s_i),$$
其中 $d\ge 1$ 通常是离散整数，也可以在连续情况下定义相应的概率密度。若是离散时长，则要求
$$\sum _{d=1}^{\infty }{p}_i(d)=1.$$
在实际实现中，也可限定最大持续时长 $D_{\max }$，则
$$\sum _{d=1}^{D_{\max }}{p}_i(d)=1.$$

观测生成模型

如果某一段处于状态 $s_i$ 并持续 $d$ 步，那么这段时间内产生了 { o t , o t + 1 , … , o t + d − 1 } \{o_t, o_{t+1}, \dots, o_{t+d-1}\} {ot​,ot+1​,…,ot+d−1​} 这 $d$ 个观测。常见假设有：

1. 独立观测：在状态 $s_i$ 内的每个时刻都独立地生成观测，即
   $$P(o_t,o_{t+1},\dots ,o_{t+d-1}|X_{t:t+d-1}=s_i)=\prod _{\tau =0}^{d-1}{b}_i(o_{t+\tau }).$$
2. 段级观测：将一段(状态$s_i$持续$d$步)看成一个整体，使用 ${\varphi }_i\left(o_{t:t+d-1}\right)$ 来刻画生成概率。

下文主要示例采用独立观测的形式，便于与标准HMM做类比。

---

半隐马尔可夫模型的联合概率

如果把时间 $1$ 到 $T$ 的整个隐藏状态过程，拆成若干段：
$$(s_{i_1},d_1)\to (s_{i_2},d_2)\to \dots \to (s_{i_M},d_M),$$
其中 ${\sum }_{m=1}^M{d}_m=T$，那么这代表：

• 第1段在状态 $s_{i_1}$ 上停留 $d_1$ 步，产生前 $d_1$ 个观测；
• 第2段在状态 $s_{i_2}$ 上停留 $d_2$ 步，产生接下来 $d_2$ 个观测；
• … 直到总长度 $T$。

如果 $\pi =({\pi }_1,\dots ,{\pi }_N)$ 是初始状态分布，$a_{ij}$ 是从状态 $s_i$ 转到状态 $s_j$ 的概率，$p_i(d)$ 是状态 $i$ 的时长分布，$b_i(o)$ 是状态 $i$ 下观测概率（离散）或密度（连续），则该段式序列对应的联合概率：

$$\begin{array}{rl}& P(O,(s_{i_1},d_1),(s_{i_2},d_2),\dots ,(s_{i_M},d_M)\mid \lambda )\\ =& \underset{\text{初始状态}}{\underbrace{{\pi }_{i_1}}}\times \underset{\text{持续时间分布}}{\underbrace{p_{i_1}(d_1)}}\times \prod _{\tau =1}^{d_1}\underset{\text{观测概率}}{\underbrace{b_{i_1}(o_{\tau })}}\\ & \times \underset{\text{状态转移}}{\underbrace{a_{i_1{i}_2}}}\times p_{i_2}(d_2)\times \prod _{\tau =d_1+1}^{d_1+d_2}{b}_{i_2}(o_{\tau })\times \dots \\ & \times a_{i_{M-1}{i}_M}\times p_{i_M}(d_M)\times \prod _{\tau =T-d_M+1}^T{b}_{i_M}(o_{\tau }).\end{array}$$

对所有可能的分段组合进行求和(或积分)，即可得到 $P(O\mid \lambda )$，但直接枚举通常会有指数复杂度，因此提出了半隐马尔可夫模型的前向-后向算法来高效求解。

---

半隐马尔可夫模型的前向-后向算法

前向变量及其公式推导

在标准HMM中，前向变量 ${\alpha }_t(i)$ 表示观测到 $o_1,\dots ,o_t$ 且时刻 $t$ 处于状态 $i$ 的联合概率。
在Semi-HMM中，因为状态 $i$ 一旦进入，有可能停留多步，必须记录“已经在 $i$ 上停留了多少步”。因此，我们定义一个带持续时间标记的前向量：

$${\alpha }_t(i,d)=P(o_1,\dots ,o_t,X_t=s_i,\text{该状态已持续了 }d\text{ 步}|\lambda ).$$
其中 $1\le d\le D_{\max }$或直到 $t$。在具体实现中为了索引方便，$d$ 也可能从0开始。

初始化

对于 $t=1$（第一个观测），若我们要表示 ${\alpha }_1(i,1)$，即可写为
$${\alpha }_1(i,1)={\pi }_i{p}_i(1)b_i(o_1).$$
因为第一时刻就进入状态 $i$，并且停留时长为1的概率是 ${\pi }_i\cdot p_i(1)$，还要乘以观测 $o_1$ 的生成概率 $b_i(o_1)$。

其它 ${\alpha }_1(i,d)$ 对于 $d>1$ 则为 0，因为在第1时刻不可能说“已经持续 $>1$ 步”了。

递推

当 $t\ge 2$ 时，需要分两种情形来更新 ${\alpha }_t(i,d)$：

1. 继续留在相同状态 $i$ 并且持续时间数加1：
   $${\alpha }_{t-1}(i,d-1)⟶{\alpha }_t(i,d),d>1.$$
   如果在时刻 $t-1$ 仍是状态 $i$ 并且已连续 $(d-1)$ 步，则下一步就变成“在同一状态 $i$ 持续 $d$ 步”。因为这里是同一段内部继续，所以只需乘上一帧(或一步)的观测概率 $b_i(o_t)$。
   因此对这一条更新：
   $${\alpha }_t(i,d)+={\alpha }_{t-1}(i,d-1)\times b_i(o_t).$$

2. 从其他状态 $j$ 结束其持续后，转移到新状态 $i$ 的持续时间 $d=1$：
   $$\sum _{j=1}^N\sum _{r=1}^{D_{\max }}{\alpha }_{t-1}(j,r)\times a_{ji}\times p_i(1)\times b_i(o_t).$$
   意思是：若在上一时刻 $t-1$ 处于状态 $j$ 并且已经在 $j$ 上持续了 $r$ 步，那么此时段结束，转移到 $i$ ，于是现时刻 $t$ 刚开始进入状态 $i$（持续时长=1）。需要乘以转移概率 $a_{ji}$，再乘上 $p_i(1)$（开始一段在 $i$ 上持续时长为1的概率），并且当前观测是 $o_t$，故还要乘 $b_i(o_t)$。
   对应数学写作：
   $${\alpha }_t(i,1)=\sum _{j=1}^N\sum _{r=1}^{D_{\max }}{\alpha }_{t-1}(j,r)a_{ji}{p}_i(1)b_i(o_t).$$

归纳起来，可以用以下公式统一表示：

$${\alpha }_t(i,d)=\{\begin{array}{llc}\sum _{j=1}^N\sum _{r=1}^{D_{\max }}{\alpha }_{t-1}(j,r)a_{ji}{p}_i(1)b_i\left(o_t\right),& d=1,{\alpha }_{t-1}(i,d-1)b_i\left(o_t\right),& d>1.\end{array}$$

边界条件

当 $t<d$ 时，“已持续 $d$ 步”就不可能，所以一般在实现中要对下标进行相应检查。或设定 $d\le t$ 以保证逻辑一致。

后向变量及其公式推导

类似地可以定义后向变量 ${\beta }_t(i,d)$，表示：“在时刻 $t$ 处于状态 $i$ 并已持续 $d$ 步，且从 $t+1$ 到 $T$ 的观测序列为 $(o_{t+1},\dots ,o_T)$ 的概率(条件于$\lambda$并给定到 $t$ 的信息)”。

可以通过在 $t=T$ 时做初始化：${\beta }_T(i,d)=1$，然后从后往前递推，类似标准HMM。但要分情况“继续在同一状态”或者“转移到下一个状态并相应乘以持续时间分布”等。
举例而言，若我们当前在 $i$ 上还可以再持续 $d^{\prime }$ 步，不断乘以观测概率；或者从 $i$ 转移到别的状态 $j$ 并开始新的持续段 $\dots$ 相应写出递推公式。由于与前向类似，只是逆向，本文不详述全部推导，思路类似。

观测序列概率的计算

一旦拿到前向量 ${\alpha }_t(i,d)$ 或后向量 ${\beta }_t(i,d)$，在 $t=T$ 处把所有可能的“状态+时长”可能性都求和，即可得到 $P(O\mid \lambda )$：

$$P(O\mid \lambda )=\sum _{i=1}^N\sum _{d=1}^{\min (T,D_{\max })}{\alpha }_T(i,d).$$
如果使用后向变量，则从初始时刻相结合来计算也可以得到同样结果。

---

持续时间分布与参数学习

半隐马尔可夫模型在持续时间分布方面非常灵活，可以使用：

• 离散分布 $p_i(d)$
• 高斯分布、Gamma分布、Poisson分布等连续分布，对应在具体实施时要处理积分或者离散化
• 混合分布 等

学习（训练）

当要从数据中估计半隐马尔可夫模型参数 $\lambda =(A,p_i(\cdot ),b_i(\cdot ),\pi )$ 时，可以利用与 HMM 类似的EM 算法(或Baum-Welch思路)，但需要在“E步”中借助前向-后向变量（带duration标识），计算相应的后验概率，然后在“M步”中更新：

1. 转移概率 $a_{ij}$。
2. 持续时间分布 $p_i(d)$。
3. 观测概率 $b_i(o)$。
4. 初始分布 $\pi$。

更新公式与标准HMM类似，只是多了对“时长”以及段内观测的统计计数。每次迭代会使对数似然函数 $ \ln P(O \mid \lambda)$ 单调上升，最终收敛到局部极值。

---

小结与应用场景

1. 小结

   • 半隐马尔可夫模型在状态停留时长上更加灵活，可显式定义分布 $p_i(d)$。
   • 与标准HMM相比，它摆脱了几何分布的限制，往往更符合很多现实序列中“停留时间不一定是几何形”的情况。
   • 前向-后向算法、解码算法（类似维特比）都需要多一维“持续时间”标记，公式推导与实现相对复杂。

2. 典型应用场景

   • 语音识别：音素常持续多帧，且时长有一定分布特性。
   • 生物信息学：基因/蛋白质片段中功能区段的长度分布。
   • 故障检测：故障状态持续时长往往不是几何分布。

---

代码示例

下面是一段简化的Python示例，用离散观测 + 离散时长分布展示如何实现一个半隐马尔可夫模型的前向算法，并据此计算 $P(O|\lambda )$。该示例未包含后向算法、解码（维特比）或参数学习（Baum-Welch），仅示范核心的“带持续时间”前向过程。

import numpy as np

class SemiHMM:
    def __init__(self, A, D, B, pi):
        """
        A: (N x N) 状态转移矩阵
           A[i, j] = P(从状态 i 转移到状态 j)
        D: (N x D_max) 持续时间分布
           D[i, d] = P(在状态 i 上停留 (d+1) 步)
        B: (N x M) 观测概率矩阵
           B[i, k] = P(观测值 k | 状态 i)
        pi: (N,) 初始状态概率向量
        """
        self.A = np.array(A)
        self.D = np.array(D)
        self.B = np.array(B)
        self.pi = np.array(pi)
        
        self.N = self.A.shape[0]         # 状态数 N
        self.D_max = self.D.shape[1]    # 最大持续时间
        self.M = self.B.shape[1]        # 观测值种类数
        
    def forward(self, O):
        """
        计算 P(O | 模型) 的简化前向算法, O 是长度 T 的离散观测序列索引
        返回: P(O | 模型)
        """
        T = len(O)
        
        # alpha[t, i, d] = 在时刻 t 处于状态 i, 且在 i 上已持续 d 步 的联合概率
        # d 范围 1..D_max(或不超过t本身)
        alpha = np.zeros((T, self.N, self.D_max))
        
        # 初始化 (t=0)
        # t=0 代表第一个观测, 对应 alpha[0, i, d=1]
        for i in range(self.N):
            # (d=1) => D[i, 0], 并乘以 pi[i] 和 B[i, O[0]]
            alpha[0, i, 0] = self.pi[i] * self.D[i, 0] * self.B[i, O[0]]
        
        # 递推
        for t in range(1, T):
            for i in range(self.N):
                for d in range(self.D_max):
                    # 如果 d=0 => 表示刚切换到 i 的持续时间=1
                    if d == 0:
                        # 累加来自所有 (j,r), r=1..D_max, alpha[t-1, j, r] * A[j,i] 
                        temp_sum = 0.0
                        for j in range(self.N):
                            for r in range(self.D_max):
                                temp_sum += alpha[t-1, j, r] * self.A[j, i]
                        # 刚进入i => p_i(1)=D[i,0], 再乘 B[i,O[t]]
                        alpha[t, i, d] = temp_sum * self.D[i, 0] * self.B[i, O[t]]
                    else:
                        # d>0 => 在 i 上继续停留
                        # alpha[t, i, d] = alpha[t-1, i, d-1] * B[i,O[t]]
                        alpha[t, i, d] = alpha[t-1, i, d-1] * self.B[i, O[t]]
        
        # 终止: 在 t=T-1 时, 所有 (i, d) 的概率和
        return np.sum(alpha[T-1, :, :])


if __name__ == '__main__':
    # 示例: 2个状态, 最大持续时间D_max=3, 观测值有2种(索引0或1)
    A = [
        [0.6, 0.4],  # 状态0 -> 状态0 / 1
        [0.3, 0.7]   # 状态1 -> 状态0 / 1
    ]
    # D[i, d] = P(在状态 i 上停留 d+1 步)
    # 例如 d=0 => 停留1步; d=1 => 停留2步; d=2 => 停留3步
    D = [
        [0.5, 0.3, 0.2],  # 状态0
        [0.2, 0.3, 0.5]   # 状态1
    ]
    # B[i, k] = P(观测值 k | 状态 i)
    B = [
        [0.6, 0.4],  # 状态0
        [0.3, 0.7]   # 状态1
    ]
    # 初始状态概率
    pi = [0.7, 0.3]
    
    model = SemiHMM(A, D, B, pi)
    
    # 观测序列
    O = [0, 1, 1, 0]  # 示例: 长度4
    
    prob_O = model.forward(O)
    print(f"观测序列 {O} 的概率 P(O|模型) = {prob_O}")