HMM-EM算法估计参数实现

作者: IntoTheVoid | 来源:发表于2021-03-04 20:16 被阅读0次

Learning Problem: HMM Training

Learning Problem: $P(\theta|D), 其中\theta=\{A, B\}$ , D是训练数据
HMM的训练过程需要使用Forward/Backward算法用于求隐状态的期望
HMM的训练过程也需要使用EM算法，根据F/B算法得到的隐状态期望以及D求关于 $\theta$ 的最大似然估计

极大似然估计示例

假设隐状态已知，取值分别是{H， G}，且观测变量已知，取值分别是{1, 2, 3}, 共有四条数据，如下：

image.png

如何通过上面的示例数据计算最大似然估计呢？？？

1.估计π（隐状态的初始概率分布）

思路：根据上述示例先初始化π，只需要通过统计每个隐状态在序列中出现在首位的次数/总的序列个数，所以例如此处，可以计算得到：

$\pi_{H} = \frac{1}{4} \\ \pi_{G} = \frac{3}{4}$

2.估计A（转移概率矩阵）

image.png

$A=\begin{bmatrix} 2/4&2/4\\3/4&1/4 \end{bmatrix}$

3.估计B（发射概率矩阵）

image.png

$B=\begin{bmatrix} 4/6&2/6&0/6\\1/6&2/6&3/6 \end{bmatrix}$

极大似然估计总结：上述的计算过程，仅当已知隐藏状态序列时，以上最大似然估计才有效。但是，事实并非如此，我们并不知道隐藏状态。因此，我们需要找到另一种方法来估算过A，B。

EM算法（原理）

那么如何使用EM算法来对A， B， π估计呢？

初始化 $[A, B]$ .初始化所有概率相等或者随机初始化
基于1中计算每个过渡/发射使用频率的期望值。用计算结果来估计隐变量 $[\xi,\gamma]$ ，这里面用到了forward/backward算法
基于2中计算的隐变量重新估计概率 $[A, B]$ ，这里面用到了最大似然估计算法
重复2,3直至收敛

概率方法

1.推导 $\hat{a_{ij}}$

idea：如果对于任意一个序列，我们知道每个t时刻，从状态i到状态j的概率，我们对该序列所有时刻的统计加和即可以得到从状态i到状态j的概率

即从数学表述，给定 $\theta, V^T$ ,在t时刻的状态为i，t+1时刻状态为j的概率, 其中 $V^T$ 表示长度为T的观测序列

$p(s(t) = i, s(t+1) = j | V_T, \theta)$

根据条件概率公式，有如下结论

$p(X, Y) = p(X|Y)p(Y)$ ,如果把 $Y|Z$ 看成整体，可以推导 $p(X, Y|Z) = p(X|Y, Z)p(Y|Z)$ ，进一步得到 $p(X|Y,Z) = \frac{p(X,Y|Z)}{p(Y|Z)}$

所以对上式可以进一步改写成

表达方式一：
$p(s(t) = i, s(t+1) = j | V^T, \theta) = \frac{p(s(t)=i, s(t+1)=j, V^T|\theta)}{p(V^T|\theta)}$

当然也可以通过另一种表达方式书写上面的式子,上式中 $V^T$ 即下式中的 $X$ , 上式中 $s(t)=i,s(t+1)=j$ 即下式中的 $z_k$ ,其中下式中忽略了 $\theta$
表达方式二：

image.png

也就是说不管是表达式1还是表达式2，我们都需要分别计算分子部分和分母部分，

1）分子部分是可以通过forward/backward算法进行计算的，即forward的 $\alpha_{i}(t)$ 与 $\beta_{j}(t+1)$ ，因为是两个时刻所以还需要计算两个时刻的转移概率，以及 $t+1$ 时刻的发射概率，那么最终得到的分子部分计算公式如下：

$p(s(t)=i, s(t+1)=j, V^T|\theta) = \alpha_{i}(t)*a_{ij}*b_{jkv(t+1)}*\beta_{j}(t+1)$

image.png

2）分母部分 $p(V^T|\theta)$ 是给定观测序列 $V^T$ 以及模型参数 $\theta$ 的概率，那么对于该观测序列的某个时刻t，假设隐状态一共有M种取值，对于t时刻每个状态i,都有M中可能的下一个时刻t+1的状态j，而每个t时刻都有M种取值，故而根据边缘概率计算得到如下公式

$p(V^T|\theta) = \sum_{i=1}^{M}\sum_{j=1}^{M}\alpha_{i}(t)*a_{ij}*b_{jkv(t+1)}*\beta_{j}(t+1)$

我们定义 $\xi_{ij}(t) = p(s(t) = i, s(t+1) = j | V^T, \theta)$ , 那么可以得到

$\xi_{ij}(t) = \frac{\alpha_{i}(t)*a_{ij}*b_{jkv(t+1)}*\beta_{j}(t+1)}{\sum_{i=1}^{M}\sum_{j=1}^{M}\alpha_{i}(t)*a_{ij}*b_{jkv(t+1)}*\beta_{j}(t+1)}$

注意：上述的 $\xi_{ij}(t)$ 只是对于某个观测序列而言的某一个时刻，我们需要对该观测序列所有的时刻都执行上面的操作，然后求和, 但是求和后需要进行归一化使其称为概率，我们知道了分子是从状态i到状态j,那么分母可以是状态i确定的情况下，下一时刻出现的不同的状态j

进一步推导可得：

$\hat{a_{ij}}=\frac{\sum_{t=1}^{T-1}\xi_{ij}(t)}{\sum_{t=1}^{T-1}\sum_{j=1}^{M}\xi_{ij}(t)} \tag{1}$

3）分母的另一种解释
假设给定某个序列 $V^T$ 以及模型参数，我们可以表示在t时刻隐状态为 $i$ 的概率为

注：下面的公式中用到了条件概率和联合概率公式的转换，注意区分
$\begin{equation} \label{eq01} \begin{split} p(s(t)=i|V^T, \theta) & = \frac{p(s(t)=i, V^T|\theta)}{p(V^T|\theta)} \\ & = \frac{p(v(1)...v(t), s(t)=i|\theta)p(v(t+1)...v(T)|s(t)=i,\theta)}{p(V^T|\theta)} \\ & = \frac{\alpha_{i}(t)\beta_{i}(t)}{p(V^T|\theta)} \\ & = \frac{\alpha_{i}(t)\beta_{i}(t)}{\sum_{i=1}^{M}\alpha_{i}(t)\beta_{i}(t)}=\gamma_{i}(t) \end{split} \end{equation}$

image.png

如果使用上述的表述重新表达 $\hat{a_{ij}}$ ,可以得到
$\hat{a_{ij}}=\frac{\sum_{t=1}^{T-1}\xi_{ij}(t)}{\sum_{t=1}^{T-1}\gamma_{i}(t)}\tag{2}$

通过公式（1）和公式（2）对比可以得到 $\gamma_{i}(t)=\sum_{j=1}^{M}\xi_{ij}(t)$

也就是说两种角度都可以求得我们最终的目标 $\hat{a_{ij}}$

2.推导 $\hat{b_{jk}}$

给定隐状态 $j$ 时，生成 $v_k$ 的概率可以用 $b_{jk}$ 表示,同时我们通过上面的分析已知，
在t时刻状态为j的概率为
$\gamma_{j}(t)=\frac{\alpha_{j}(t)\beta_{j}(t)}{\sum_{t=1}^{T}\gamma_{j}(t)}$
那么要计算 $\hat{b_{jk}}$ 相当于求给定状态为j时的观测变量为 $v(t)=k$ 的概率，即如下表示

$\hat{b_{jk}}=\frac{\sum_{t=1}^{T}\gamma_{j}(t)I(v(t)=k)}{\sum_{t=1}^{T}\gamma_{j}(t)}\tag{3}$

其中 $I(v(t)=k)$ 表示当t时刻为k时才为1，否则为0，即示性函数

EM算法（伪代码）

initialize A and B(random initialize or all equal)
while not convergence:
- E-Step
  - $\xi_{ij}(t) = \frac{\alpha_{i}(t)*a_{ij}\;*b_{jkv(t+1)}\;\;\;\;*\beta_{j}(t+1)}{\sum_{i=1}^{M}\;\;\sum_{j=1}^{M}\;\;\alpha_{i}(t)*a_{ij}\;*b_{jkv(t+1)}\;\;\;\;*\beta_{j}(t+1)}$
  - $\gamma_{i}(t)=\sum_{j=1}^{M}\xi_{ij}(t)$
- M-Step
  - $\hat{a_{ij}}=\frac{\sum_{t=1}^{T-1}\;\;\xi_{ij}(t)}{\sum_{t=1}^{T-1}\;\;\sum_{j=1}^{M}\;\xi_{ij}(t)}$
  - $\hat{b_{jk}}=\frac{\sum_{t=1}^{T}\;\;\gamma_{j}(t)I(v(t)=k)}{\sum_{t=1}^{T}\;\;\gamma_{j}(t)}$
return A, B

代码实现


import pandas as pd
import numpy as np
 
 
def forward(V, a, b, initial_distribution):
    alpha = np.zeros((V.shape[0], a.shape[0]))
    alpha[0, :] = initial_distribution * b[:, V[0]]
 
    for t in range(1, V.shape[0]):
        for j in range(a.shape[0]):
            # Matrix Computation Steps
            #                  ((1x2) . (1x2))      *     (1)
            #                        (1)            *     (1)
            alpha[t, j] = alpha[t - 1].dot(a[:, j]) * b[j, V[t]]
 
    return alpha
 
 
def backward(V, a, b):
    beta = np.zeros((V.shape[0], a.shape[0]))
 
    # setting beta(T) = 1
    beta[V.shape[0] - 1] = np.ones((a.shape[0]))
 
    # Loop in backward way from T-1 to
    # Due to python indexing the actual loop will be T-2 to 0
    for t in range(V.shape[0] - 2, -1, -1):
        for j in range(a.shape[0]):
            beta[t, j] = (beta[t + 1] * b[:, V[t + 1]]).dot(a[j, :])
 
    return beta
 
 
def em_learning(V, a, b, initial_distribution, n_iter=100):
    M = a.shape[0]
    T = len(V)
 
    for n in range(n_iter):
        alpha = forward(V, a, b, initial_distribution)
        beta = backward(V, a, b)
 
        xi = np.zeros((M, M, T - 1))
        for t in range(T - 1):
            denominator = np.dot(np.dot(alpha[t, :].T, a) * b[:, V[t + 1]].T, beta[t + 1, :])
            for i in range(M):
                numerator = alpha[t, i] * a[i, :] * b[:, V[t + 1]].T * beta[t + 1, :].T
                xi[i, :, t] = numerator / denominator
 
        gamma = np.sum(xi, axis=1)
        a = np.sum(xi, 2) / np.sum(gamma, axis=1).reshape((-1, 1))
 
        # Add additional T'th element in gamma
        gamma = np.hstack((gamma, np.sum(xi[:, :, T - 2], axis=0).reshape((-1, 1))))
 
        K = b.shape[1]
        denominator = np.sum(gamma, axis=1)
        for l in range(K):
            b[:, l] = np.sum(gamma[:, V == l], axis=1)
 
        b = np.divide(b, denominator.reshape((-1, 1)))
 
    return {"a":a, "b":b}

data = pd.read_csv('data/data_python.csv')
 
V = data['Visible'].values
 
# Transition Probabilities
a = np.ones((2, 2))
A = a / np.sum(a, axis=1)
 
# Emission Probabilities
b = np.array(((1, 3, 5), (2, 4, 6)))
B = b / np.sum(b, axis=1).reshape((-1, 1))
 
# Equal Probabilities for the initial distribution
π = np.array((0.5, 0.5))

print(em_learning(V, A, B, π, n_iter=100))

{'a': array([[0.53816345, 0.46183655],
       [0.48664443, 0.51335557]]), 'b': array([[0.16277513, 0.26258073, 0.57464414],
       [0.2514996 , 0.27780971, 0.47069069]])}

HMM-EM算法估计参数实现

Learning Problem: HMM Training

极大似然估计示例

1.估计π（隐状态的初始概率分布）

2.估计A（转移概率矩阵）

3.估计B（发射概率矩阵）

EM算法（原理）

概率方法

1.推导 $\hat{a_{ij}}$

2.推导 $\hat{b_{jk}}$

EM算法（伪代码）

代码实现

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读