HMM学习笔记

HMM学习笔记

0 前言#

隐马尔可夫模型（Hidden Markov Model,HMM）是一种经典的统计学模型，根本上属于传统机器学习中的概率图模型。学习该模型的前置知识为随机过程。

1 定义#

1.1 概念引入#

上图展示了隐马尔可夫链的基本元素和基本结构。隐马尔可夫的基本元素包含状态序列(state sequence)和观测序列(obervation sequence)，状态序列一般不可知，或者说被隐藏起来了，故名隐马尔可夫链，观测序列一般可知。从结构上可知，状态变量和观测变量会随着时间而变化。描述这个变化的过程需要引入状态转移概率矩阵A和观测概率矩阵B。

1.2 变量定义#

设状态变量 $s_t$ 的值域为 $Q=\{q_1,q_2,\cdots,q_t,\cdots,q_N\}$ ，观测变量 $o_t$ 的值域为 $V=\{v_1,v_2,\cdots,v_t,\cdots,v_M\}$ ， $N=length(Q)$ ， $M=length(V)$ 。

设序列长度为 $T$ ，则状态序列 $S=(s_1,s_2,\cdots,s_t,\cdots,s_T)$ ，观测序列 $O=(o_1,o_2,\cdots,o_t,\cdots,o_T)$ 。

定义状态转移概率矩阵

A = [a_{i j}]_{N \times N}

$A=[a_{ij}]_{N\times N}$

其中， $a_{ij}=p(s_{t+1}=q_j|s_t=q_i),\qquad i=1,2,\cdots,N;j=1,2,\cdots,N$ 。 $a_{ij}$ 是在时刻 $t$ 处于状态 $q_i$ 条件下在时刻 $t+1$ 转移到状态 $q_j$ 的概率。

定义观测概率矩阵

B = [b_{j} (k)]_{N \times M}

$B=[b_j(k)]_{N\times M}$

其中， $b_j(k)=P(o_t=v_k|s_t=q_j),\qquad k=1,2,\cdots,M;j=1,2,\cdots,N$ ， $b_j(k)$ 是在时刻 $t$ 处于状态 $q_j$ 的条件下观测 $v_k$ 的概率。

定义初始概率向量

π = [π_{i}]_{N \times 1}

$\pi=[\pi_i]_{N\times 1}$

其中， $\pi_i=P(s_1=q_i),\qquad i=1,2,\cdots,N$ ， $\pi_i$ 是时刻 $t=1$ 处于状态 $q_i$ 的概率。

由此，可构建隐马尔可夫模型参数 $\lambda=(A,B,\pi)$

1.3 两个基本假设#

(i)齐次马尔可夫性假设，即假设隐藏的马尔可夫链在任意时刻 $t$ 的状态只依赖于前一时刻的状态，于其他时刻的状态及观测无关，也与时刻 $t$ 无关。

$p(s_{t}|s_{t-1},\cdots,s_1,o_t,o_{t-1},\cdots,o_1)=p(s_{t}|s_{t-1})$
(ii)观测独立性假设，即假设任意时刻的观测只依赖于该时刻的马尔可夫链的状态，与其他观测及状态无关。

$p(o_t|s_T,s_{T-1},\cdots,s_1,o_T,o_{T-1},\cdots,o_1)=p(o_t|s_t)$

1.4 举例解释#

假设现有{盒子1，盒子2，盒子3，盒子4}，每个盒子内含有若干数量不等的红球和白球。现在，观测者首先选择一个盒子，然后从该盒子里抽取一个球，获得该球的颜色，然后操作员会更换盒子，继续让观测者抽球。说明：盒子的更换规则观测者不知情。

那么，在这个例子中，观测者会获得一个关于小球颜色的序列，该序列就是观测序列；操作员依次更换盒子，操作员会得到一个关于盒子序号的序列，该序列就是对于观测者来讲就是状态序列，各个盒子中抽取红白球的概率会组成观测概率矩阵，盒子的更换方式会组成状态转移概率矩阵，观测者选择第一个盒子的概率分布称为初始状态概率向量。

详细例子见《统计学习方法》，李航，清华大学出版社，P173

1.5 三个基本问题#

Evaluation：已知 $O,\lambda$ ，求解 $p(O|\lambda)$ ；
Learning：已知 $O$ ，求 $\lambda$ ，使 $\lambda_{MLE}=\mathop{argmax}\limits_{\lambda}p(O|\lambda)$
Decoding：已知 $O,\lambda$ ，求 $S$ ，使 $\hat{S}=\mathop{argmax}\limits_{S}p(S|O,\lambda)$

2 Evaluation Problem#

2.1 直接计算法#

\begin{aligned} p (O | λ) & = \sum_{S} p (S, O | λ) \\ = \sum_{S} p (O | S, λ) p (S | λ) \end{aligned}

$\begin{aligned} p(O|\lambda) &= \sum\limits_{S}p(S,O|\lambda) \\ &= \sum\limits_{S}p(O|S,\lambda)p(S|\lambda) \end{aligned}$

注：
边缘概率可以展开为联合概率的和；
联合概率可以展开为含有条件概率的形式；如 $P(AB)=P(A|B)P(B)$

根据基本假设(i)，可得

\begin{aligned} p (S | λ) & = p (s_{1}, s_{2}, \dots, s_{t} | λ) \\ = p (s_{t} | s_{1}, s_{2}, \dots, s_{t - 1}, λ) p (s_{1}, s_{2}, \dots, s_{t - 1} | λ) \\ = a_{s_{t - 1} s_{t}} a_{s_{t - 2} s_{t - 1}} \dots a_{s_{1} s_{2}} \\ = π_{s_{1}} \prod_{t = 2}^{T} a_{s_{t - 1}, s_{t}} \end{aligned}

$\begin{aligned} p(S|\lambda)&= p(s_1,s_2,\cdots,s_t|\lambda)\\ &= p(s_t|s_1,s_2,\cdots,s_{t-1},\lambda)p(s_1,s_2,\cdots,s_{t-1}|\lambda)\\ &= a_{s_{t-1}s_t} a_{s_{t-2}s_{t-1}} \cdots a_{s_1s_2}\\ &= \pi_{s_1}\prod\limits_{t=2}^Ta_{s_{t-1},s_t} \end{aligned}$

注：其中， $a_{s_{t-1},s_t}$ 表示从状态 $s_{t-1}$ 到状态 $s_t$ 的转移概率，为了和1.2中的变量a_{ij}相对应，可以令 $i=locate(s_{t-1})$ ， $j=locate(s_t)$ ，对应 $s_{t-1}$ 和 $s_t$ 在 $Q$ (定义见1.2)中的位置序号，则 $a_{s_{t-1},s_t}$ 可以用 $a_{ij}$ 表示，上式为了表述清晰，便直接使用状态变量代替了位置序号，与定义中的公式表达略有不同。

根据基本假设(ii)，对于状态序列 $S=(s_1,s_2,\cdots,s_t,\cdots,s_T)$ ，有

p (O | S, λ) = \prod_{t = 1}^{T} b_{s_{t}} (o_{t})

$p(O|S,\lambda)=\prod\limits_{t=1}^Tb_{s_t}(o_t)$

注：上式与 $a_{ij}$ 类似，使用状态变量 $s_t$ 表示状态变量在 $Q$ 中的位置序号，使用观测变量 $o_t$ 表示观测变量在 $V$ (定义见1.2)中的位置序号。

于是有

p (O | λ) = \sum_{S} π_{s_{1}} \prod_{t = 2}^{T} a_{s_{t - 1}, s_{t}} \prod_{t = 1}^{T} b_{s_{t}} (o_{t})

$p(O|\lambda)=\sum\limits_{S}\pi_{s_1}\prod\limits_{t=2}^Ta_{s_{t-1},s_t}\prod\limits_{t=1}^Tb_{s_t}(o_t)$

对上述 $\sum\limits_{S}$ 做解释：这里的意思是对后面的含有 $s_1,s_2,\cdots,s_T$ 这些 $T$ 个随机变量分别求均值再求和，即对状态变量求T次均值，每个状态变量都有 $N$ 个取值，也就是说，计算上式需要进行 $N^T$ 个均值计算再求和，计算时间复杂度为 $O(N^T)$ 。随着 $T$ 的增加，时间计算复杂度指数级增长。

2.2 前向算法#

给定 $\lambda$ ， $O=(o_1,o_2,\cdots,o_t)$ ，定义：到时刻 $t$ 观测序列为 $o_1,o_2,\cdots,o_t$ 且状态为 $q_i$ 的概率为前向概率，记为

α_{t} (i) = p (o_{1}, o_{2}, \dots, o_{t}, s_{t} = q_{i} | λ)

$\alpha_t(i)=p(o_1,o_2,\cdots,o_t,s_t=q_i|\lambda)$

所以，

p (O | λ) = \sum_{i = 1}^{N} p (O, s_{t} = q_{i} | λ) = \sum_{i = 1}^{N} α_{t} (i)

$p(O|\lambda)=\sum\limits_{i=1}^Np(O,s_t=q_i|\lambda)=\sum\limits_{i=1}^N\alpha_t(i)$

因此，欲求时刻t的 $p(O|\lambda)$ ，即求 $\alpha_t(i)$ 。可以通过递归的思想求解 $\alpha_t(i)$

\begin{aligned} α_{t} (j) & = p (o_{1}, o_{2}, \dots, o_{t}, s_{t} = q_{j} | λ) \\ = \sum_{i = 1}^{N} p (o_{1}, o_{2}, \dots, o_{t}, s_{t - 1} = q_{i}, s_{t} = q_{j} | λ) \\ = \sum_{i = 1}^{N} p (o_{t} | o_{1}, o_{2}, \dots, s_{t - 1} = q_{i}, s_{t} = q_{j} | λ) p (o_{1}, \dots, o_{t - 1}, s_{t - 1} = q_{i}, s_{t} = q_{j} | λ) \\ = \sum_{i = 1}^{N} p (o_{t} | s_{t} = q_{j}) p (o_{1}, \dots, o_{t}, s_{t - 1} = q_{i}, s_{t} = q_{j} | λ) \\ = \sum_{i = 1}^{N} p (o_{t} | s_{t} = q_{j}) p (s_{t} = q_{j} | o_{1}, \dots, o_{t}, s_{t - 1} = q_{i}, λ) p (o_{1}, \dots, o_{t}, s_{t - 1} = q_{i} | λ) \\ = \sum_{i = 1}^{N} p (o_{t} | s_{t} = q_{j}) p (s_{t} = q_{j} | s_{t - 1} = q_{i}) p (o_{1}, \dots, o_{t}, s_{t - 1} = q_{i} | λ) \\ = \sum_{i = 1}^{N} b_{j} (o_{t}) a_{i j} α_{t - 1} (i) \end{aligned}

$\begin{aligned} \alpha_{t}(j)&=p(o_1,o_2,\cdots,o_{t},s_{t}=q_j|\lambda)\\ &=\sum\limits_{i=1}^Np(o_1,o_2,\cdots,o_{t},s_{t-1}=q_i,s_t=q_j|\lambda)\\ &=\sum\limits_{i=1}^Np(o_{t}|o_1,o_2,\cdots,s_{t-1}=q_i,s_{t}=q_j|\lambda)p(o_1,\cdots,o_{t-1},s_{t-1}=q_i,s_{t}=q_j|\lambda) \\ &=\sum\limits_{i=1}^Np(o_{t}|s_{t}=q_j)p(o_1,\cdots,o_t,s_{t-1}=q_i,s_{t}=q_j|\lambda) \\ &=\sum\limits_{i=1}^Np(o_{t}|s_{t}=q_j)p(s_{t}=q_j|o_1,\cdots,o_t,s_{t-1}=q_i,\lambda)p(o_1,\cdots,o_t,s_{t-1}=q_i|\lambda) \\ &=\sum\limits_{i=1}^Np(o_{t}|s_{t}=q_j)p(s_{t}=q_j|s_{t-1}=q_i)p(o_1,\cdots,o_t,s_{t-1}=q_i|\lambda) \\ &=\sum\limits_{i=1}^Nb_{j}(o_t)a_{ij}\alpha_{t-1}(i) \end{aligned}$

可以看出来，上式依旧是使用了2.1使用过技巧

将联合概率拆解为条件概率及其条件的概率的乘积
利用基本假设(i)(ii)

因此，给定 $\lambda=(A,B,\pi), O=(o_1,o_2,\cdots,o_T)$ ，可得 $\alpha_t(i)$ 的初始条件 $\alpha_1(i)$

α_{1} (i) = π_{i} b_{i} (o_{1})

$\alpha_1(i)=\pi_ib_i(o_1)$

然后可以利用推导出的递归表达式计算出 $\alpha_t(1),\alpha_t(2), \cdots,\alpha_t(T)$ ，然后求和即可得到 $p(O|\lambda)$ 。在每个递归式中，需要计算 $N$ 个 $\alpha_{t-1}(i),\{i=1,2,\cdots,N\}$ ，获得 $N$ 个 $\alpha_{t}(i),\{i=1,2,\cdots,N\}$ ，则每个递归式的时间复杂度为 $O(N^2)$ ，又观测序列的长度为 $T$ ，即需要进行T次递归，则该算法的时间复杂度为 $O(N^2T)$ 。

3 Decoding Problem#

decoding 问题就是根据观测序列求解最有可能的状态序列。

该问题的求解主要涉及到Viterbi Algorithm（维特比算法）。维特比算法是一种动态规划算法。它用于寻找最有可能产生观测事件序列的维特比路径——隐含状态序列。

定义：观测序列为 $o_1,o_2,\cdots,o_t$ ，模型参数为 $\lambda$ ，到时刻 $t$ ，状态变量 $s_t$ 为 $q_i$ 的概率最大，概率大小记为 $\delta_t(i)$ 。

δ_{t} (i) = max_{s_{1}, \dots, s_{t - 1}} p (o_{1}, \dots, o_{t}, s_{1}, \dots, s_{t - 1}, s_{t} = q_{i} | λ)

$\delta_{t}(i)=\max\limits_{s_1,\cdots,s_{t-1}}p(o_1,\cdots,o_t,s_1,\cdots,s_{t-1},s_t=q_i|\lambda)$

相应的，可以推出

δ_{t + 1} (j) = max_{1 \leq i \leq N} [δ_{t} (i) a_{i j}] b_{j} (o_{t + 1})

$\delta_{t+1}(j)=\max\limits_{1\le i\le N}[\delta_t(i)a_{ij}]b_j(o_{t+1})$

定义在时刻 $t+1$ 状态为 $j$ 的所有单个路径 $(j_1,j_2,\cdots,j_t,j_{t+1})$ 中概率最大的路径的第 $t+1$ 个节点为

ψ_{t + 1} (j) = \underset{1 \leq i \leq N}{a r g m a x} [δ_{t} (i) a_{i j}]

$\psi_{t+1}(j)=\mathop{argmax}\limits_{1\le i\le N}[\delta_t(i)a_{ij}]$

代码示例

物理背景同1.4

clc
clear 
close all

A = [0.5 0.2 0.3;
     0.3 0.5 0.2;
     0.2 0.3 0.5];
B = [0.5 0.5;
     0.4 0.6;
     0.7 0.3];
PI = [0.2; 0.4; 0.4];

lambda = {A, B, PI};

%红-1 白-2
O = [1 2 1];
%代码中的输入观测序列和输出状态序列需进行合理的数学编码，例如上述O代表{红，白，红}，输出状态序列同理
my_viterbi(lambda, O)

function S_decode = my_viterbi(lambda, O)
    A = lambda{1, 1};
    B = lambda{1, 2};
    PI = lambda{1, 3};
    numO = length(O);
    S_decode = zeros(1, numO);

    % initialize
    N = length(PI);% number of states
    Psi = zeros(N, numO);
    Delta = zeros(N, numO);
    for i = 1:N
        Delta(i, 1) = PI(i) * B(i, O(1));
    end

    % recursive
    for k = 1:numO-1
        for i = 1:N
            temp = zeros(N,1);
            for j = 1:N
                temp(j) = Delta(j, k)*A(j,i);
            end
            [temp_max, max_index] = max(temp);
            Delta(i, k+1) = temp_max*B(i, O(k+1));
            Psi(i, k+1) = max_index;
        end
    end

    % end
    [temp_max, max_index] = max(Delta(:, end));
    P_decode_path = temp_max;
    tail_decode_point = max_index;

    % path traceback
    S_decode(numO) = tail_decode_point;
    for k = numO-1:-1:1
        S_decode(k) = Psi(S_decode(k+1), k+1);
    end
end