1 """
2 This package includes code for representing and learning HMM's.
3
4 Most of the code in this package was derived from the descriptions provided in
5 'A Tutorial on Hidden Markov Models and Selected Applications in Speach
6 Recognition' by Lawence Rabiner.
7
8 Conventions:
9 The keyword argument elem_size will be passed in when creating numpy array
10 objects.
11 """
12 import math,random,sys
13 import numpy
14
15 -def calcalpha(stateprior,transition,emission,observations,numstates,elem_size=numpy.longdouble):
16 """
17 Calculates 'alpha' the forward variable.
18
19 The alpha variable is a numpy array indexed by time, then state (TxN).
20 alpha[t][i] = the probability of being in state 'i' after observing the
21 first t symbols.
22 """
23 alpha = numpy.zeros((len(observations),numstates),dtype=elem_size)
24 for x in xrange(numstates):
25 alpha[0][x] = stateprior[x]*emission[x][observations[0]]
26 for t in xrange(1,len(observations)):
27 for j in xrange(numstates):
28 for i in xrange(numstates):
29 alpha[t][j] += alpha[t-1][i]*transition[i][j]
30 alpha[t][j] *= emission[j][observations[t]]
31 return alpha
32
33 -def forwardbackward(stateprior,transition,emission,observations,numstates,elem_size=numpy.longdouble):
34 """
35 Calculates the probability of a sequence given the HMM.
36 """
37 alpha = calcalpha(stateprior,transition,emission,observations,numstates,elem_size)
38 return sum(alpha[-1])
39
40 -def calcbeta(transition,emission,observations,numstates,elem_size=numpy.longdouble):
41 """
42 Calculates 'beta' the backward variable.
43
44 The beta variable is a numpy array indexed by time, then state (TxN).
45 beta[t][i] = the probability of being in state 'i' and then observing the
46 symbols from t+1 to the end (T).
47 """
48 beta = numpy.zeros((len(observations),numstates),dtype=elem_size)
49 for s in xrange(numstates):
50 beta[len(observations)-1][s] = 1.
51 for t in xrange(len(observations)-2,-1,-1):
52 for i in xrange(numstates):
53 for j in xrange(numstates):
54 beta[t][i] += transition[i][j]*emission[j][observations[t+1]]*beta[t+1][j]
55 return beta
56
57 -def calcxi(stateprior,transition,emission,observations,numstates,alpha=None,beta=None,elem_size=numpy.longdouble):
58 """
59 Calculates 'xi', a joint probability from the 'alpha' and 'beta' variables.
60
61 The xi variable is a numpy array indexed by time, state, and state (TxNxN).
62 xi[t][i][j] = the probability of being in state 'i' at time 't', and 'j' at
63 time 't+1' given the entire observation sequence.
64 """
65 if alpha is None:
66 alpha = calcalpha(stateprior,transition,emission,observations,numstates,elem_size)
67 if beta is None:
68 beta = calcbeta(transition,emission,observations,numstates,elem_size)
69 xi = numpy.zeros((len(observations),numstates,numstates),dtype=elem_size)
70 for t in xrange(len(observations)-1):
71 denom = 0.0
72 for i in xrange(numstates):
73 for j in xrange(numstates):
74 thing = 1.0
75 thing *= alpha[t][i]
76 thing *= transition[i][j]
77 thing *= emission[j][observations[t+1]]
78 thing *= beta[t+1][j]
79 denom += thing
80 for i in xrange(numstates):
81 for j in xrange(numstates):
82 numer = 1.0
83 numer *= alpha[t][i]
84 numer *= transition[i][j]
85 numer *= emission[j][observations[t+1]]
86 numer *= beta[t+1][j]
87 xi[t][i][j] = numer/denom
88 return xi
89
90 -def calcgamma(xi,seqlen,numstates, elem_size=numpy.longdouble):
91 """
92 Calculates 'gamma' from xi.
93
94 Gamma is a (TxN) numpy array, where gamma[t][i] = the probability of being
95 in state 'i' at time 't' given the full observation sequence.
96 """
97 gamma = numpy.zeros((seqlen,numstates),dtype=elem_size)
98 for t in xrange(seqlen):
99 for i in xrange(numstates):
100 gamma[t][i] = sum(xi[t][i])
101 return gamma
102
103 -def baumwelchstep(stateprior,transition,emission,observations,numstates,numsym,elem_size=numpy.longdouble):
104 """
105 Given an HMM model and a sequence of observations, computes the Baum-Welch
106 update to the parameters using gamma and xi.
107 """
108 xi = calcxi(stateprior,transition,emission,observations,numstates,elem_size=elem_size)
109 gamma = calcgamma(xi,len(observations),numstates,elem_size)
110 newprior = gamma[0]
111 newtrans = numpy.zeros((numstates,numstates),dtype=elem_size)
112 for i in xrange(numstates):
113 for j in xrange(numstates):
114 numer = 0.0
115 denom = 0.0
116 for t in xrange(len(observations)-1):
117 numer += xi[t][i][j]
118 denom += gamma[t][i]
119 newtrans[i][j] = numer/denom
120 newemiss = numpy.zeros( (numstates,numsym) ,dtype=elem_size)
121 for j in xrange(numstates):
122 for k in xrange(numsym):
123 numer = 0.0
124 denom = 0.0
125 for t in xrange(len(observations)):
126 if observations[t] == k:
127 numer += gamma[t][j]
128 denom += gamma[t][j]
129 newemiss[j][k] = numer/denom
130 return newprior,newtrans,newemiss
131
133 """
134 A class for modeling and learning HMMs.
135
136 This class conveniently wraps the module level functions. Class objects hold 6
137 data members:
138 - num_states number of hidden states in the HMM
139 - num_symbols number of possible symbols in the observation
140 sequence
141 - precision precision of the numpy.array elements (defaults to
142 longdouble)
143 - prior The prior probability of starting in each state
144 (Nx1 array)
145 - transition_matrix The probability of transitioning between each state
146 (NxN matrix)
147 - emission_matrix The probability of each symbol in each state
148 (NxO matrix)
149 You can set the 3 matrix parameters as you wish, but make sure the shape of
150 the arrays matches num_states and num_symbols, as these are used internally
151
152 Typical usage of this class is to create an HMM with a set number of states
153 and external symbols, train the HMM using addEvidence(...), and then use
154 the sequenceProb(...) method to see how well a specific sequence matches
155 the trained HMM.
156 """
157 - def __init__(self,num_states,num_symbols,init_type='uniform',precision=numpy.longdouble):
158 """
159 Creates a new HMMLearner object with the given number of internal
160 states, and external symbols.
161
162 calls self.reset(init_type=init_type)
163 """
164 self.num_states = num_states
165 self.num_symbols = num_symbols
166 self.precision = precision
167 self.reset(init_type=init_type)
168
169 - def reset(self, init_type='uniform'):
170 """
171 Resets the 3 arrays using the given initialization method.
172
173 Wipes out the old arrays. You can use this method to change the shape
174 of the arrays by first changing num_states and/or num_symbols, and then
175 calling this method.
176
177 Currently supported initialization methods:
178 uniform prior, transition, and emission probabilities are all
179 uniform (default)
180 """
181 if init_type == 'uniform':
182 self.prior = numpy.ones( (self.num_states), dtype=self.precision) *(1.0/self.num_states)
183 self.transition_matrix = numpy.ones( (self.num_states,self.num_states), dtype=self.precision)*(1.0/self.num_states)
184 self.emission_matrix = numpy.ones( (self.num_states,self.num_symbols), dtype=self.precision)*(1.0/self.num_symbols)
185
187 """
188 Returns the probability that this HMM generated the given sequence.
189
190 Uses the forward-backward algorithm. If given an array of
191 sequences, returns a 1D array of probabilities.
192 """
193 if len(newData.shape) == 1:
194 return forwardbackward( self.prior,\
195 self.transition_matrix,\
196 self.emission_matrix,\
197 newData,\
198 self.num_states,\
199 self.precision)
200 elif len(newData.shape) == 2:
201 return numpy.array([forwardbackward(self.prior,self.transition_matrix,self.emission_matrix,newSeq,self.num_states,self.precision) for newSeq in newData])
202
203 - def addEvidence(self, newData, iterations=1,epsilon=0.0):
204 """
205 Updates this HMMs parameters given a new set of observed sequences
206 using the Baum-Welch algorithm.
207
208 newData can either be a single (1D) array of observed symbols, or a 2D
209 matrix, each row of which is a seperate sequence. The Baum-Welch update
210 is repeated 'iterations' times, or until the sum absolute change in
211 each matrix is less than the given epsilon. If given multiple
212 sequences, each sequence is used to update the parameters in order, and
213 the sum absolute change is calculated once after all the sequences are
214 processed.
215 """
216 if len(newData.shape) == 1:
217 for i in xrange(iterations):
218 newp,newt,newe = baumwelchstep( self.prior, \
219 self.transition_matrix, \
220 self.emission_matrix, \
221 newData, \
222 self.num_states, \
223 self.num_symbols,\
224 self.precision)
225 pdiff = sum([abs(np-op) for np in newp for op in self.prior])
226 tdiff = sum([sum([abs(nt-ot) for nt in newti for ot in oldt]) for newti in newt for oldt in self.transition_matrix])
227 ediff = sum([sum([abs(ne-oe) for ne in newei for oe in olde]) for newei in newe for olde in self.emission_matrix])
228 if(pdiff < epsilon) and (tdiff < epsilon) and (ediff < epsilon):
229 break
230 self.prior = newp
231 self.transition_matrix = newt
232 self.emission_matrix = newe
233 else:
234 for i in xrange(iterations):
235 for sequence in newData:
236 newp,newt,newe = baumwelchstep( self.prior, \
237 self.transition_matrix, \
238 self.emission_matrix, \
239 sequence, \
240 self.num_states, \
241 self.num_symbols,\
242 self.precision)
243 self.prior = newp
244 self.transition_matrix = newt
245 self.emission_matrix = newe
246 pdiff = sum([abs(np-op) for np in newp for op in self.prior])
247 tdiff = sum([sum([abs(nt-ot) for nt in newti for ot in oldt]) for newti in newt for oldt in self.transition_matrix])
248 ediff = sum([sum([abs(ne-oe) for ne in newei for oe in olde]) for newei in newe for olde in self.emission_matrix])
249 if(pdiff < eps) and (tdiff < eps) and (ediff < eps):
250 break
251 self.prior = newp
252 self.transition_matrix = newt
253 self.emission_matrix = newe
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290