模型结构
隐马尔科夫模型如下图所示，
(s)我们是无法直接观测到的，我们得到的是对state的一种度量，记为(y)。模型涉及两个部分：一个是从(s_{t-1} o s_t)的转移，另一个是(s_t o y_t)由state生成观测变量(y)的部分。
HMM主要涉及以下三个问题：
- 在已知(s_{t-1} o s_t)和(s_{t} o y_t)两个分布的情况下计算(P(X))，即似然
- 根据得到的观测数据计算(s_{t-1} o s_t)和(s_{t} o y_t)两个分布的模型参数
- 预测，即已知模型分布参数和观测数据情况下，推测对应的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()