看完这个你就基本理解HMM是神马东西了
在此之前可以了解下维特比算法
Vertibe
用python代码实现HMM模型:
utility.py
#-*- coding:utf-8
import numpy as np
"""
generateIndex from hiddenState or ovservationsState
Attributes
----------
states : hiddenStates = ("Healthy","Fever")
Returns
----------
labelIndex : {'Healthy': 0, 'Fever': 1}
"""
def generateStatesIndex(lables):
labelIndex = {}
i = 0
for l in lables:
labelIndex[l] = i
i += 1
return labelIndex
"""
generate matrix from map
Attributes
----------
map : A = {
"Healthy":{"Healthy":0.7,"Fever":0.3},
"Fever":{"Healthy": 0.4, "Fever": 0.6}
}
Returns
----------
matrixA : [[ 0.7 0.3]
[ 0.4 0.6]]
"""
def generateMatrix(map, labelIndex1, labelIndex2):
m = np.zeros((len(labelIndex1),len(labelIndex2)),dtype=float)
for row in map:
for col in map[row]:
m[labelIndex1[row]][labelIndex2[col]] = map[row][col]
return m
"""
generate matrix from map
Attributes
----------
map : pi = {"Healthy":0.6,"Fever":0.4}
}
Returns
----------
list : pi = [0.6,0.4]
"""
def generatePiVector(map,labelIndex):
pi = np.zeros((len(labelIndex)),dtype=float)
for key in map:
pi[labelIndex[key]] = map[key]
return pi
hmm.py
#-*- coding:utf-8
import utility
import numpy as np
class HMM:
"""
Order 1 Hidden Markov Model
Attributes
----------
A : numpy.ndarray
State transition probability matrix
B: numpy.ndarray
Output emission probability matrix with shape(N, number of output types)
pi: numpy.ndarray
Initial state probablity vector
"""
def __init__(self, A, B, pi):
self.A = A
self.B = B
self.pi = pi
"""
generate matrix from map
Attributes
----------
T : the length of the output
Returns
----------
hiddenStates,observationStates
"""
def generateData(self,T):
#根据分布列表,返回可能返回的Index
def _getFromProbs(probs):
return np.where(np.random.multinomial(1, probs) == 1)[0][0]
hiddenStates = np.zeros(T,dtype=int)
observationsStates = np.zeros(T,dtype=int)
hiddenStates[0] = _getFromProbs(self.pi) #产生第一个hiddenStates
observationsStates[0] = _getFromProbs(self.B[hiddenStates[0]]) #产生第一个observationStates
for t in range(1,T):
hiddenStates[t] = _getFromProbs(self.A[hiddenStates[t-1]])
observationsStates[t] = _getFromProbs((self.B[hiddenStates[t]]))
return hiddenStates,observationsStates
#维特比算法进行预测,即解码
def viterbi(self,observationsSeq):
T = len(observationsSeq)
N = len(self.pi)
prePath = np.zeros((T,N),dtype=int)
dpMatrix = np.zeros((T,N),dtype=float)
dpMatrix[0,:] = self.pi * self.B[:,observationsSeq[0]]
for t in range(1,T):
for n in range(N):
probs = dpMatrix[t-1,:] * self.A[:,n] * self.B[n,observationsSeq[t]]
prePath[t,n] = np.argmax(probs)
dpMatrix[t,n] = np.max(probs)
maxProb = np.max(dpMatrix[T-1,:])
maxIndex = np.argmax(dpMatrix[T-1,:])
path = [maxIndex]
for t in reversed(range(1,T)):
path.append(prePath[t,path[-1]])
path.reverse()
return maxProb,path
#计算公式中的alpha二维数组
def _forward(self,observationsSeq):
T = len(observationsSeq)
N = len(self.pi)
alpha = np.zeros((T,N),dtype=float)
alpha[0,:] = self.pi * self.B[:,observationsSeq[0]] #numpy可以简化循环
for t in range(1,T):
for n in range(0,N):
alpha[t,n] = np.dot(alpha[t-1,:],self.A[:,n]) * self.B[n,observationsSeq[t]] #使用内积简化代码
return alpha
#计算公式中的beita二维数组
def _backward(self,observationsSeq):
T = len(observationsSeq)
N = len(self.pi)
beta = np.zeros((T,N),dtype=float)
beta[T-1,:] = 1
for t in reversed(range(T-1)):
for n in range(N):
beta[t,n] = np.sum(self.A[n,:] * self.B[:,observationsSeq[t+1]] * beta[t+1,:])
return beta
#前后向算法学习参数
def baumWelch(self,observationsSeq,criterion=0.001):
T = len(observationsSeq)
N = len(self.pi)
while True:
# alpha_t(i) = P(O_1 O_2 ... O_t, q_t = S_i | hmm)
# Initialize alpha
alpha = self._forward(observationsSeq)
# beta_t(i) = P(O_t+1 O_t+2 ... O_T | q_t = S_i , hmm)
# Initialize beta
beta = self._backward(observationsSeq)
#根据公式求解XIt(i,j) = P(qt=Si,qt+1=Sj | O,λ)
xi = np.zeros((T-1,N,N),dtype=float)
for t in range(T-1):
denominator = np.sum( np.dot(alpha[t,:],self.A) * self.B[:,observationsSeq[t+1]] * beta[t+1,:])
for i in range(N):
molecular = alpha[t,i] * self.A[i,:] * self.B[:,observationsSeq[t+1]] * beta[t+1,:]
xi[t,i,:] = molecular / denominator
#根据xi就可以求出gamma,注意最后缺了一项要单独补上来
gamma = np.sum(xi,axis=2)
prod = (alpha[T-1,:] * beta[T-1,:])
gamma = np.vstack((gamma, prod /np.sum(prod)))
newpi = gamma[0,:]
newA = np.sum(xi,axis=0) / np.sum(gamma[:-1,:],axis=0).reshape(-1,1)
newB = np.zeros(self.B.shape,dtype=float)
for k in range(self.B.shape[1]):
mask = observationsSeq == k
newB[:,k] = np.sum(gamma[mask,:],axis=0) / np.sum(gamma,axis=0)
if np.max(abs(self.pi - newpi)) < criterion and \
np.max(abs(self.A - newA)) < criterion and \
np.max(abs(self.B - newB)) < criterion:
break
self.A,self.B,self.pi = newA,newB,newpi
if __name__ == '__main__':
hiddenStates = ("Healthy", "Fever")
observationsStates = ("normal", "cold", "dizzy")
pi = {"Healthy": 0.8, "Fever": 0.2}
A = {
"Healthy": {"Healthy": 0.6, "Fever": 0.4},
"Fever": {"Healthy": 0.3, "Fever": 0.7}
}
B = {
"Healthy": {"normal": 0.6, "cold": 0.2, "dizzy": 0.2},
"Fever": {"normal": 0.1, "cold": 0.5, "dizzy": 0.4},
}
hStatesIndex = utility.generateStatesIndex(hiddenStates)
oStatesIndex = utility.generateStatesIndex(observationsStates)
A = utility.generateMatrix(A,hStatesIndex,hStatesIndex)
B = utility.generateMatrix(B,hStatesIndex,oStatesIndex)
pi = utility.generatePiVector(pi,hStatesIndex)
h = HMM(A=A,B=B,pi=pi)
states_data,observations_data = h.generateData(20)
print (states_data)
print (observations_data)
# h.baumWelch(observations_data)
# prob,path = h.viterbi(observations_data)
# p = 0.0
#
# print path
# print states_data
运行结果: