隐马尔可夫模型是关于时序的概率模型，描述由一个隐藏的马尔可夫链生成不可观测的状态随机序列，再由各个状态生成一个观测而产生观测随机序列的过程。隐藏的马尔可夫链随机生成的状态的序列，称为状态序列（state sequence);每个状态生成的一个观测，而由此产生的观测的随机序列，称为观测序列（observation sequence）。序列的每一个位置可以看作是一个时刻。
隐马尔可夫模型由初始状态概率分布向量（ $\pi$），状态转移概率矩阵（$A$）和状态产生的观测概率矩阵（$B$）决定，$\pi$和$A$决定状态序列，$B$决定状态产生的观测序列，因此假设马尔可夫模型为$\lambda$，则： $$\lambda=(A,B,\pi)$$
假设
$Q$是所有可能隐藏状态的集合 ,其中N为状态个数，则 $$Q=\{q_1,q_2,...,q_N\}$$ $V$是所有可能观测的集合，其中M为可能观测数，则 $$V=\{v_1,v_2,...,v_M\}$$ $I$是长度为T的状态序列，$O$是对应的观测序列，则 $$I=\{i_1,i_2,...,i_T\},O=\{o_1,o_2,...,o_T\}$$
那么隐马尔可夫的参数$A$、$B$、$\pi$的定义如下： $$A=[a_{ij}]_{NxN}$$ 其中$a_{ij}=P(i_{t+1}=q_j|i_t=q_i)$,i=1,2,...,N;j=1,2,...,N，表示时刻t处于状态$q_i$的条件下在时刻t+1转移到$q_j$的概率 $$B=[b_j(k)]_{NxM}$$ 其中$b_j(k)=P(o_t=v_k|i_t=q_j)$,k=1,2,...,M;j=1,2,...,N,表示时刻t处于状态$q_j$的条件下生成观测$v_k$的概率 $$\pi=(\pi_i)$$ 其中$\pi_i=P(i_1=q_i)$，i=1,2,...,N,是时刻t=1处于状态$q_i$的概率
具体请参见斯坦福的材料
预测问题就是寻找最有可能产生可观擦序列的隐藏状态序列，维特比算法和forward算法比较相似，维特比算法在每一步取最大的转移概率，而foward在每一步对概率求和。 假设： $$\delta_t(i)=\underset{q_1,q_2,...,q_{t-1}}{\max}P(q_1q_2...q_t=s_i,o_1,o_2,...o_t|\lambda)$$
维特比（Viterbi）算法：
1.初始化 $$\delta_1(i)=\pi_ib_i(o_1),1\leq i \leq N,\psi_1(i)=0$$
2.递归 $$\delta_t(j)=\underset{1 \leq i \leq N}{\max}[\delta_{t-1}(i)a_{ij}]b_j(o_t),2 \leq t \leq T,1 \leq j \leq N$$
$$\psi_t(j)=\arg\underset{1 \leq i \leq N}{\max}[\delta_{t-1}(i)a_{ij}],2 \leq t \leq T,1 \leq j \leq N$$
3.中止 $$P=\underset{1 \leq i \leq N}{\max}[\delta_T(i)]$$ $$q_T=\arg \underset{1 \leq i \leq N}{\max}[\delta_T(i)]$$
4.回溯(backtracking) $$q_t=\psi_{t+1}(q_{t+1}),t=T-1,T-2,...,1$$
import numpy as np
def Viterbi(pi,a,b,obs):
"""
pi:初始概率
a:状态转移矩阵
b:发射概率矩阵
obs：可观察系列
"""
#隐藏状态
nStates = np.shape(b)[0]
#时间步长
T = np.shape(obs)[0]
#隐藏状态序列
path = np.zeros(T)
delta = np.zeros((nStates,T))
phi = np.zeros((nStates,T))
#obs[0]表示第一个可观擦的状态
#b[:,obs[0]]表示每种隐藏状态转为第一可观察状态的概率
delta[:,0] = pi * b[:,obs[0]]
phi[:,0] = 0
for t in range(1,T):
for s in range(nStates):
#取每步状态的概率值，然后转到下一步的可观察状态
delta[s,t] = np.max(delta[:,t-1]*a[:,s])*b[s,obs[t]]
#保留每步最大概率的状态
phi[s,t] = np.argmax(delta[:,t-1]*a[:,s])
path[T-1] = np.argmax(delta[:,T-1])
for t in range(T-2,-1,-1):
path[t] = phi[int(path[t+1]),t+1]
return path,delta, phi
def test():
np.random.seed(4)
pi = np.array([0.25,0.25,0.25,0.25])
aLast = np.array([0.25,0.25,0.25,0.25])
#a = np.array([[.7,.3],[.4,.6]] )
a = np.array([[.4,.3,.1,.2],[.6,.05,.1,.25],[.7,.05,.05,.2],[.3,.4,.25,.05]])
#b = np.array([[.2,.4,.4],[.5,.4,.1]] )
b = np.array([[.2,.1,.2,.5],[.4,.2,.1,.3],[.3,.4,.2,.1],[.3,.05,.3,.35]])
obs = np.array([0,0,3,1,1,2,1,3])
#obs = np.array([2,0,2])
path,delta,phi=Viterbi(pi,a,b,obs)
print("path=",path)
print("delta=",delta)
print("phi=",phi)
test()
forward算法计算过程：
1.初始化 $$\alpha_1(i)=\pi_ib_i(o_1),1 \leq i \leq N$$
2.归纳 $$\alpha_{t+1}(j)=[\sum_{i=1}^N\alpha_t(i)a_{ij}]b_j(o_{t+1}),1 \leq t \leq T-1,1 \leq j \leq N$$
3.中止 $$P(O|\lambda)=\sum_{i=1}^N\alpha_T(i)$$
scaling = False
def HMMfwd(pi,a,b,obs):
nStates = np.shape(b)[0]
T = np.shape(obs)[0]
alpha = np.zeros((nStates,T))
alpha[:,0] = pi*b[:,obs[0]]
for t in range(1,T):
for s in range(nStates):
alpha[s,t] = b[s,obs[t]] * np.sum(alpha[:,t-1] * a[:,s])
c = np.ones((T))
if scaling:
for t in range(T):
c[t] = np.sum(alpha[:,t])
alpha[:,t] /= c[t]
return alpha,c
def HMMbwd(a,b,obs,c):
nStates = np.shape(b)[0]
T = np.shape(obs)[0]
beta = np.zeros((nStates,T))
beta[:,T-1] = 1.0 #aLast
for t in range(T-2,-1,-1):
for s in range(nStates):
beta[s,t] = np.sum(b[:,obs[t+1]] * beta[:,t+1] * a[s,:])
for t in range(T):
beta[:,t] /= c[t]
#beta[:,0] = b[:,obs[0]] * np.sum(beta[:,1] * pi)
return beta
def BaumWelch(obs,nStates):
"""
obs:可观察序列
nStates：状态序列
"""
T = np.shape(obs)[0]
xi = np.zeros((nStates,nStates,T))
# Initialise pi, a, b randomly
pi = 1./nStates*np.ones((nStates))
a = np.random.rand(nStates,nStates)
b = np.random.rand(nStates,np.max(obs)+1)
tol = 1e-5
error = tol+1
maxits = 100
nits = 0
while ((error > tol) & (nits < maxits)):
nits += 1
oldpi = pi.copy()
olda = a.copy()
oldb = b.copy()
# E step
alpha,c = HMMfwd(pi,a,b,obs)
beta = HMMbwd(a,b,obs,c)
for t in range(T-1):
for i in range(nStates):
for j in range(nStates):
xi[i,j,t] = alpha[i,t]*a[i,j]*b[j,obs[t+1]]*beta[j,t+1]
xi[:,:,t] /= np.sum(xi[:,:,t])
# The last step has no b, beta in
for i in range(nStates):
for j in range(nStates):
xi[i,j,T-1] = alpha[i,T-1]*a[i,j]
xi[:,:,T-1] /= np.sum(xi[:,:,T-1])
# M step
for i in range(nStates):
pi[i] = np.sum(xi[i,:,0])
for j in range(nStates):
a[i,j] = np.sum(xi[i,j,:T-1])/np.sum(xi[i,:,:T-1])
for k in range(max(obs)):
found = (obs==k).nonzero()
b[i,k] = np.sum(xi[i,:,found])/np.sum(xi[i,:,:])
error = (np.abs(a-olda)).max() + (np.abs(b-oldb)).max()
print(nits, error, 1./np.sum(1./c), np.sum(alpha[:,T-1]))
return pi, a, b
def biased_coins():
a = np.array([[0.4,0.6],[0.9,0.1]])
b = np.array([[0.49,0.51],[0.85,0.15]])
pi = np.array([0.5,0.5])
obs = np.array([0,1,1,0,1,1,0,0,1,1,0,1,1,1,0,0,1,0,0,1,1,0,1,1,1,1,0,1,0,0,1,0,1,0,0,1,1,1,0])
print( Viterbi(pi,a,b,obs)[0])
print (BaumWelch(obs,2))
biased_coins()