HMM

模型结构

隐马尔科夫模型如下图所示，

$s$我们是无法直接观测到的，我们得到的是对state的一种度量，记为$y$。模型涉及两个部分：一个是从$s_{t-1}\to s_t$的转移，另一个是$s_t\to y_t$由state生成观测变量$y$的部分。

HMM主要涉及以下三个问题：

1. 在已知$s_{t-1}\to s_t$和$s_t\to y_t$两个分布的情况下计算$P(X)$，即似然
2. 根据得到的观测数据计算$s_{t-1}\to s_t$和$s_t\to y_t$两个分布的模型参数
3. 预测，即已知模型分布参数和观测数据情况下，推测对应的state信息

文章基本是参考资料^[1]的概括和代码的简化整理。

似然计算

如果直接用$P(X)=\sum _{Z}P(X,Z)$来计算的话，时间复杂度大约是$O(N^T)$，其中$X$表示观测变量，$Z$表示隐变量，N表示$Z$的可能取值，$T$表示时间长度。

一般都会用$\alpha -\beta$算法来简化计算，$\alpha ,\beta$定义如下

$$\begin{array}{rl}{\alpha }_t& =P(X_{1:t},Z_t)\\ {\beta }_t& =P(X_{t+1:N}|Z_t)\end{array}$$

可以看到$\sum _{Z_t=i}{\alpha }_t\times {\beta }_t=P(X)$

${\alpha }_t$和${\beta }_t$的递推关系如下

$$\begin{array}{rl}{\alpha }_t& =\sum _{z_{t-1}}{\alpha }_{t-1}P(Z_t|Z_{t-1})P(X_t|Z_t)\\ {\beta }_t& =\sum _{z_{t+1}}{\beta }_{t+1}P(Z_{t+1}|Z_t)P(X_{t+1}|Z_{t+1})\end{array}$$

可以看到$\alpha$是从前往后，而$\beta$是从后往前，这样给定初始值的情况下，就可以计算下去。

$$\begin{array}{rl}{\alpha }_0& ={\pi }_0\ast P(X_0|Z_0)\\ {\beta }_N& =\text{bold}1\end{array}$$

代码如下:

def evaluate_alpha(self, X):
    '''
    alpha_t = P(X_{1:t}, s_t)
    :param X: observed list
    :return: likelihood of X given parameters
    '''
    alpha = self.pi * self.B[:, X[0]]

    for t in range(1, len(X)):
        alpha = np.sum(alpha.reshape(-1, 1) * self.A, axis=0) * self.B[:, X[t]]

    return  alpha.sum()

def evaluate_beta(self, X):
    '''
    beta_t = P(X_{t+1:N} | s_t)
    :param X: observed list
    :return: beta_1
    '''
    beta = np.ones(self.N)
    for t in X[:0:-1]:
        beta = np.sum(self.A * beta * self.B[:, t], axis=1)
    return np.sum(beta * self.pi * self.B[:, X[0]])

参数推断

关于参数推断部分，套用$EM$框架，再利用Lagrange即可，详细步骤可以查阅参考资料^[1]，参数更新公式如下

$$\begin{array}{rl}{\pi }_i^{t+1}& =\frac{{\alpha }_1(i){\beta }_1(i)}{P(X|)}\\ a_{i,j}& =\frac{\sum _{t=1}^{T-1}{a}_t(i){\beta }_{t+1}(j)a_{i,j}{b}_j(x_{t+1})}{\sum _{t=1}^{T-1}{\alpha }_t(i){\beta }_t(i)}\\ b_{j,k}& =\frac{\sum _{t=1}^T{\alpha }_t(j){\beta }_t(j)I(x_t=k)}{\sum _{t=1}^T{\alpha }_t(j){\beta }_t(j)}\end{array}$$

文章代码里关于$\beta$的更新是element-wise的形式，其实可以写成向量形式的。推断部分的代码如下

def fit(self, X):
    '''
    learn parameter
    :param X: observed list
    :return: updated parameter
    '''
    self.pi = np.random.rand(self.N)
    self.pi = self.pi / self.pi.sum()
    self.A = np.ones((self.N, self.N)) / self.N
    self.B = np.ones((self.N, self.M)) / self.M
    T = len(X)

    for _ in range(50):
        alpha, beta = self.alpha_beta(X)
        gamma = alpha * beta
        # update A
        for i in range(self.N):
            for j in range(self.N):
                self.A[i, j] =  np.sum(alpha[:-1, i] * beta[1:, j] * self.A[i, j] * self.B[j, X[1:]]) / gamma[:-1, i].sum()
        # update B
        for i in range(self.N):
            for j in range(self.M):
                self.B[i, j] =  np.sum(gamma[:, i] * (X == j)) / gamma[:, i].sum()
        self.pi = (alpha[0] * beta[0]) / np.sum(alpha[0] * beta[0])


def alpha_beta(self, X):
    T = len(X)
    alpha = np.zeros((T, self.N))
    beta = np.ones((T, self.N))
    alpha[0, :]  = self.pi * self.B[:, X[0]]
    for t in range(1, T):
        alpha[t, :] = np.sum(alpha[t-1, :].reshape(-1, 1) * self.A, axis=0) * self.B[:, X[t]]
    for t in range(T-1, 0, -1):
        beta[t-1, :] = np.sum(self.A * beta[t, :] * self.B[:, X[t]], axis=1)
    return (alpha, beta)

预测

在预测state状态时，要计算全概率$P(X,Z)$的概率，代码如下

def decode(self, X):
    '''
    :param X: observation list
    :return:  hidden state list
    '''
    T = len(X)
    x = X[0]
    delta = self.pi * self.B[:,x]
    varphi = np.zeros((T, self.N), dtype=int)
    path = [0] * T
    for i in range(1, T):
        delta = delta.reshape(-1,1)     # 转成一列方便广播
        tmp = delta * self.A
        varphi[i,:] = np.argmax(tmp, axis=0)
        delta = np.max(tmp, axis=0) * self.B[:,X[i]]
    path[-1] = np.argmax(delta)
    # 回溯最优路径
    for i in range(T-1,0,-1):
        path[i-1] = varphi[i, path[i]]
    return path

全部代码如下

import os
import numpy as np

class HMM():
    def __init__(self, N, M, pi=None, A=None, B=None):
        self.N = N
        self.M = M
        self.pi = pi
        self.A = A
        self.B = B

    def get_data_with_distribution(self, dist):
        '''
        :param dist: distribution of data
        :return: sample index of the distribution
        '''
        r = np.random.rand()
        assert type(r) is not np.ndarray
        for i, p in enumerate(dist):
            print(r, p)
            if r < p:
                return i
            else:
                r -= p

    def generate(self, T:int):
        '''

        :param T: number of samples
        :return: simulated dataset
        '''
        z = self.get_data_with_distribution(self.pi)  # state
        x = self.get_data_with_distribution(self.B[z])  # observation
        result = [x]
        for _ in range(T-1):
            z = self.get_data_with_distribution(self.A[z])  # next state
            x = self.get_data_with_distribution(self.B[z])  # next observation
            result.append(x)
        return result

    def evaluate_alpha(self, X):
        '''
        alpha_t = P(X_{1:t}, s_t)
        :param X: observed list
        :return: likelihood of X given parameters
        '''
        alpha = self.pi * self.B[:, X[0]]

        for t in range(1, len(X)):
            alpha = np.sum(alpha.reshape(-1, 1) * self.A, axis=0) * self.B[:, X[t]]

        return  alpha.sum()

    def evaluate_beta(self, X):
        '''
        beta_t = P(X_{t+1:N} | s_t)
        :param X: observed list
        :return: beta_1
        '''
        beta = np.ones(self.N)
        for t in X[:0:-1]:
            beta = np.sum(self.A * beta * self.B[:, t], axis=1)
        return np.sum(beta * self.pi * self.B[:, X[0]])

    def fit(self, X):
        '''
        learn parameter
        :param X: observed list
        :return: updated parameter
        '''
        self.pi = np.random.rand(self.N)
        self.pi = self.pi / self.pi.sum()
        print(self.pi)
        self.A = np.ones((self.N, self.N)) / self.N
        self.B = np.ones((self.N, self.M)) / self.M
        T = len(X)

        for _ in range(50):
            alpha, beta = self.alpha_beta(X)
            gamma = alpha * beta
            # update A
            for i in range(self.N):
                for j in range(self.N):
                    self.A[i, j] =  np.sum(alpha[:-1, i] * beta[1:, j] * self.A[i, j] * self.B[j, X[1:]]) / gamma[:-1, i].sum()
            # update B
            for i in range(self.N):
                for j in range(self.M):
                    self.B[i, j] =  np.sum(gamma[:, i] * (X == j)) / gamma[:, i].sum()
            self.pi = (alpha[0] * beta[0]) / np.sum(alpha[0] * beta[0])


    def alpha_beta(self, X):
        T = len(X)
        alpha = np.zeros((T, self.N))
        beta = np.ones((T, self.N))
        alpha[0, :]  = self.pi * self.B[:, X[0]]
        for t in range(1, T):
            alpha[t, :] = np.sum(alpha[t-1, :].reshape(-1, 1) * self.A, axis=0) * self.B[:, X[t]]
        for t in range(T-1, 0, -1):
            beta[t-1, :] = np.sum(self.A * beta[t, :] * self.B[:, X[t]], axis=1)
        return (alpha, beta)

    def decode(self, X):
        '''

        :param X:
        :return:
        '''
        T = len(X)
        x = X[0]
        delta = self.pi * self.B[:,x]
        varphi = np.zeros((T, self.N), dtype=int)
        path = [0] * T
        for i in range(1, T):
            delta = delta.reshape(-1,1)     # 转成一列方便广播
            tmp = delta * self.A
            varphi[i,:] = np.argmax(tmp, axis=0)
            delta = np.max(tmp, axis=0) * self.B[:,X[i]]
        path[-1] = np.argmax(delta)
        # 回溯最优路径
        for i in range(T-1,0,-1):
            path[i-1] = varphi[i, path[i]]
        return path


if __name__ == "__main__":
    import matplotlib.pyplot as plt


    def triangle_data(T):  # 生成三角波形状的序列
        data = []
        for x in range(T):
            x = x % 6
            data.append(x if x <= 3 else 6 - x)
        return data


    data = np.array(triangle_data(30))
    hmm = HMM(10, 4)
    hmm.fit(data)  # 先根据给定数据反推参数
    gen_obs = hmm.generate(30)  # 再根据学习的参数生成数据
    x = np.arange(30)
    plt.scatter(x, gen_obs, marker='*', color='r')
    plt.plot(x, data, color='g')
    plt.show()

参考资料

1. 一站式解决：隐马尔可夫模型（HMM）全过程推导及实现，永远在你身后