# Code from Chapter 16 of Machine Learning: An Algorithmic Perspective (2nd Edition)
# by Stephen Marsland (http://stephenmonika.net)
# You are free to use, change, or redistribute the code in any way you wish for
# non-commercial purposes, but please maintain the name of the original author.
# This code comes with no warranty of any kind.
# Stephen Marsland, 2008, 2014
# A basic Hidden Markov Model
import numpy as np
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 Viterbi(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))
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[path[t+1],t+1]
return path,delta, phi
def BaumWelch(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 evenings():
pi = np.array([0.25, 0.25, 0.25, 0.25])
a = np.array([[0.05,0.7, 0.05, 0.2],[0.1,0.4,0.3,0.2],[0.1,0.6,0.05,0.25],[0.25,0.3,0.4,0.05]])
b = np.array([[0.3,0.4,0.2,0.1],[0.2,0.1,0.2,0.5],[0.4,0.2,0.1,0.3],[0.3,0.05,0.3,0.35]])
obs = np.array([3,1,1,3,0,3,3,3,1,1,0,2,2])
print Viterbi(pi,a,b,obs)[0]
alpha,c = HMMfwd(pi,a,b,obs)
print np.sum(alpha[:,-1])
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])
HMMfwd(pi,a,b,obs)
Viterbi(pi,a,b,obs)
print BaumWelch(obs,4)
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)