HMM

模型结构

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

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

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

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

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

似然计算

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

一般都会用(alpha-eta)算法来简化计算，(alpha,eta)定义如下

[egin{align*} alpha_t &= P(X_{1:t},Z_t)\ eta_t &= P(X_{t+1:N}vert Z_t) end{align*} ]

可以看到(sumlimits_{Z_t=i}alpha_t imeseta_t=P(X))

(alpha_t)和(eta_t)的递推关系如下

[egin{align*} alpha_t &= sum_{z_{t-1}} alpha_{t-1} P(Z_tvert Z_{t-1}) P(X_tvert Z_t)\ eta_t & =sum_{z_{t+1}} eta_{t+1} P(Z_{t+1}vert Z_t)P(X_{t+1}vert Z_{t+1}) end{align*} ]

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

[egin{align*} alpha_0 &= pi_0 * P(X_0vert Z_0) \ eta_N &= old 1 end{align*} ]

代码如下:

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]，参数更新公式如下

[egin{align*} pi_i^{t+1} &= frac{alpha_1(i)eta_1(i)}{P(Xvert)} \ a_{i,j}&= frac{sum_{t=1}^{T-1}a_t(i)eta_{t+1}(j)a_{i,j}b_j(x_{t+1})}{sum_{t=1}^{T-1}alpha_t(i)eta_t(i)}\ b_{j,k}&=frac{sum_{t=1}^Talpha_t(j)eta_t(j)I(x_t=k)}{sum_{t=1}^T alpha_t(j)eta_t(j)} end{align*} ]

文章代码里关于(eta)的更新是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）全过程推导及实现，永远在你身后