From de52ab6941d3d2463cfcc55afa462e36e3425753 Mon Sep 17 00:00:00 2001 From: Guy Zyskind Date: Wed, 14 Nov 2012 09:40:57 +0200 Subject: [PATCH] HMM, DiscreteHMM and Gaussian Mixtures HMM --- .gitignore | 3 + .project | 17 ++ .pydevproject | 10 ++ hmm/_BaseHMM.py | 286 +++++++++++++++++++++++++++++++ hmm/__init__.py | 0 hmm/continuous/GMHMM.py | 29 ++++ hmm/continuous/_ContinuousHMM.py | 175 +++++++++++++++++++ hmm/continuous/__init__.py | 0 hmm/discrete/DiscreteHMM.py | 69 ++++++++ hmm/discrete/__init__.py | 0 hmm/examples/__init__.py | 0 hmm/examples/tests.py | 109 ++++++++++++ 12 files changed, 698 insertions(+) create mode 100644 .gitignore create mode 100644 .project create mode 100644 .pydevproject create mode 100644 hmm/_BaseHMM.py create mode 100644 hmm/__init__.py create mode 100644 hmm/continuous/GMHMM.py create mode 100644 hmm/continuous/_ContinuousHMM.py create mode 100644 hmm/continuous/__init__.py create mode 100644 hmm/discrete/DiscreteHMM.py create mode 100644 hmm/discrete/__init__.py create mode 100644 hmm/examples/__init__.py create mode 100644 hmm/examples/tests.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b642889 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +*.class +*.pyc +*.pyo diff --git a/.project b/.project new file mode 100644 index 0000000..a77e0d2 --- /dev/null +++ b/.project @@ -0,0 +1,17 @@ + + + HMM + + + + + + org.python.pydev.PyDevBuilder + + + + + + org.python.pydev.pythonNature + + diff --git a/.pydevproject b/.pydevproject new file mode 100644 index 0000000..ee70624 --- /dev/null +++ b/.pydevproject @@ -0,0 +1,10 @@ + + + + + +/HMM + +python 2.7 +Default + diff --git a/hmm/_BaseHMM.py b/hmm/_BaseHMM.py new file mode 100644 index 0000000..d399066 --- /dev/null +++ b/hmm/_BaseHMM.py @@ -0,0 +1,286 @@ +''' +Created on Oct 31, 2012 + +@author: GuyZ + +The code is based on hmm.py, from QSTK. See below for license details and information. +''' +import numpy + +class _BaseHMM(object): + ''' + classdocs + ''' + + def __init__(self,n,m,precision=numpy.double,verbose=False): + self.n = n + self.m = m + + self.precision = precision + self.verbose = verbose + self.eta = self._eta1 + + def _eta1(self,t,T): + return 1. + + """ + Forward-Backward procedure is used to efficiently calculate the probability of the observation, given the model - P(O|model) + alpha_t(x) = P(O1...Ot,qt=Sx|model) - The probability of state x and the observation up to time t, given the model. + NOTE: only alpha (the forward variable) is actually required to get P(O|model) + NOTE 2: The returned value is the log-likelihood of the model. It can be greater than one + when using a continous HMM since pdfs are used, and can provide unnormalized values. + """ + def forwardbackward(self,observations): + alpha = self._calcalpha(observations) + return numpy.log(sum(alpha[-1])) + + """ + Calculates 'alpha' the forward variable. + + The alpha variable is a numpy array indexed by time, then state (TxN). + alpha[t][i] = the probability of being in state 'i' after observing the + first t symbols. + """ + def _calcalpha(self,observations): + alpha = numpy.zeros((len(observations),self.n),dtype=self.precision) + + # init stage - alpha_1(x) = pi(x)b_x(O1) + for x in xrange(self.n): + alpha[0][x] = self.pi[x]*self.B_map[x][0] + + # induction + for t in xrange(1,len(observations)): + for j in xrange(self.n): + for i in xrange(self.n): + alpha[t][j] += alpha[t-1][i]*self.A[i][j] + alpha[t][j] *= self.B_map[j][t] + + return alpha + + """ + Calculates 'beta' the backward variable. + + The beta variable is a numpy array indexed by time, then state (TxN). + beta[t][i] = the probability of being in state 'i' and then observing the + symbols from t+1 to the end (T). + """ + def _calcbeta(self,observations): + beta = numpy.zeros((len(observations),self.n),dtype=self.precision) + + # init stage + for s in xrange(self.n): + beta[len(observations)-1][s] = 1. + + # induction + for t in xrange(len(observations)-2,-1,-1): + for i in xrange(self.n): + for j in xrange(self.n): + beta[t][i] += self.A[i][j]*self.B_map[j][t+1]*beta[t+1][j] + + return beta + + """ + Find the best state sequence (path), given the model and an observation. I.E: max(P(Q|O,model)). + + This method is usually used to predict the next state after training. + """ + def decode(self, observations): + # use Viterbi's algorithm. It is possible to add additional algorithms in the future. + return self._viterbi(observations) + + """ + Find the best state sequence (path) using viterbi algorithm - a method of dynamic programming, + very similar to the forward-backward algorithm, with the added step of maximization and eventual + backtracing. + + delta[t][i] = max(P[q1..qt=i,O1...Ot|model] - the path ending in Si and until time t, + that generates the highest probability. + + psi[t][i] = argmax(delta[t-1][i]*aij) - the index of the maximizing state in time (t-1), + i.e: the previous state. + + """ + def _viterbi(self, observations): + delta = numpy.zeros((len(observations),self.n),dtype=self.precision) + psi = numpy.zeros((len(observations),self.n),dtype=self.precision) + + # init + for x in xrange(self.n): + delta[0][x] = self.pi[x]*self.B_map[x][0] + psi[0][x] = 0 + + # induction + for t in xrange(1,len(observations)): + for j in xrange(self.n): + for i in xrange(self.n): + if (delta[t][j] < delta[t-1][i]*self.A[i][j]): + delta[t][j] = delta[t-1][i]*self.A[i][j] + psi[t][j] = i + delta[t][j] *= self.B_map[j][t] + + # termination: find the maximum probability for the entire sequence (=highest prob path) + p_max = 0 # max value in time T (max) + path = numpy.zeros((len(observations)),dtype=self.precision) + for i in xrange(self.n): + if (p_max < delta[len(observations)-1][i]): + p_max = delta[len(observations)-1][i] + path[len(observations)-1] = i + + # path backtracing + path = numpy.zeros((len(observations)),dtype=self.precision) + for i in xrange(1, len(observations)): + path[len(observations)-i-1] = psi[len(observations)-i][ path[len(observations)-i] ] + + return path + + """ + Calculates 'xi', a joint probability from the 'alpha' and 'beta' variables. + + The xi variable is a numpy array indexed by time, state, and state (TxNxN). + xi[t][i][j] = the probability of being in state 'i' at time 't', and 'j' at + time 't+1' given the entire observation sequence. + """ + def _calcxi(self,observations,alpha=None,beta=None): + if alpha is None: + alpha = self._calcalpha(observations) + if beta is None: + beta = self._calcbeta(observations) + xi = numpy.zeros((len(observations),self.n,self.n),dtype=self.precision) + + for t in xrange(len(observations)-1): + denom = 0.0 + for i in xrange(self.n): + for j in xrange(self.n): + thing = 1.0 + thing *= alpha[t][i] + thing *= self.A[i][j] + thing *= self.B_map[j][t+1] + thing *= beta[t+1][j] + denom += thing + for i in xrange(self.n): + for j in xrange(self.n): + numer = 1.0 + numer *= alpha[t][i] + numer *= self.A[i][j] + numer *= self.B_map[j][t+1] + numer *= beta[t+1][j] + xi[t][i][j] = numer/denom + + return xi + + """ + Calculates 'gamma' from xi. + + Gamma is a (TxN) numpy array, where gamma[t][i] = the probability of being + in state 'i' at time 't' given the full observation sequence. + """ + def _calcgamma(self,xi,seqlen): + gamma = numpy.zeros((seqlen,self.n),dtype=self.precision) + + for t in xrange(seqlen): + for i in xrange(self.n): + gamma[t][i] = sum(xi[t][i]) + + return gamma + + """ + Updates this HMMs parameters given a new set of observed sequences. + + observations can either be a single (1D) array of observed symbols, or a 2D + matrix, each row of which is a seperate sequence. The Baum-Welch update + is repeated 'iterations' times, or until the sum absolute change in + each matrix is less than the given epsilon. If given multiple + sequences, each sequence is used to update the parameters in order, and + the sum absolute change is calculated once after all the sequences are + processed. + """ + def train(self, observations, iterations=1,epsilon=0.0001,thres=-0.001): + self._mapB(observations) + + for i in xrange(iterations): + prob_old, prob_new = self.trainiter(observations) # train iter should also update model and cacheB + + if (self.verbose): + print "iter: ", i, ", L(model|O) =", prob_old, ", L(model_new|O) =", prob_new, ", converging =", ( prob_new-prob_old > thres ) + + if ( abs(prob_new-prob_old) < epsilon ): + # converged + break + + def _updatemodel(self,new_model): + self.pi = new_model['pi'] + self.A = new_model['A'] + + def trainiter(self,observations): + # call the EM algorithm + new_model = self._baumwelch(observations) + + # calculate the log likelihood of the previous model + prob_old = self.forwardbackward(observations) + + # update the model with the new estimation + self._updatemodel(new_model) + + # map observable probabilities + self._mapB(observations) + + # calculate the log likelihood of the new model + prob_new = self.forwardbackward(observations) + + return prob_old, prob_new + + """ + new transitions probability matrix A is set to: expected_transitions(i->j)/expected_transitions(i) + """ + def _reestimateA(self,observations,xi,gamma): + A_new = numpy.zeros((self.n,self.n),dtype=self.precision) + for i in xrange(self.n): + for j in xrange(self.n): + numer = 0.0 + denom = 0.0 + for t in xrange(len(observations)-1): + numer += (self.eta(t,len(observations)-1)*xi[t][i][j]) + denom += (self.eta(t,len(observations)-1)*gamma[t][i]) + A_new[i][j] = numer/denom + return A_new + + def _calcstats(self,observations): + stats = {} + + stats['alpha'] = self._calcalpha(observations) + stats['beta'] = self._calcbeta(observations) + stats['xi'] = self._calcxi(observations,stats['alpha'],stats['beta']) + stats['gamma'] = self._calcgamma(stats['xi'],len(observations)) + + return stats + + def _reestimate(self,stats,observations): + new_model = {} + + # new init vector is set to the frequency of being in each step at t=0 + new_model['pi'] = stats['gamma'][0] + new_model['A'] = self._reestimateA(observations,stats['xi'],stats['gamma']) + + return new_model + + + """ + An EM(expectation-modification) algorithm devised by Baum-Welch. Finds a local maximum + that outputs the model that produces the highest probability, given a set of observations. + + xi[t][i][j] = P(qt=Si, qt+1=Sj|O,model) - the probability of being in state i at time t, + and in state j at time t+1, given the ENTIRE observation sequence. + + gamma[t][i] = sum(xi[i][j]) - the probability of being in state i at time t, given the ENTIRE + observation sequence. + """ + def _baumwelch(self,observations): + # E step - calculate statistics + stats = self._calcstats(observations) + + # M step + return self._reestimate(stats,observations) + + def _mapB(self,observations): + raise NotImplementedError("a mapping function for B(observable probabilities) must be implemented") + \ No newline at end of file diff --git a/hmm/__init__.py b/hmm/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/hmm/continuous/GMHMM.py b/hmm/continuous/GMHMM.py new file mode 100644 index 0000000..8e693dd --- /dev/null +++ b/hmm/continuous/GMHMM.py @@ -0,0 +1,29 @@ +''' +Created on Nov 12, 2012 + +@author: GuyZ +''' + +from _ContinuousHMM import _ContinuousHMM +import numpy + +class GMHMM(_ContinuousHMM): + ''' + classdocs + ''' + + def __init__(self,n,m,d=1,A=None,means=None,covars=None,w=None,pi=None,min_std=0.01,init_type='uniform',precision=numpy.double,verbose=False): + ''' + See ContinuousHMM constructor for more information + ''' + _ContinuousHMM.__init__(self,n,m,d,A,means,covars,w,pi,min_std,init_type,precision,verbose) #@UndefinedVariable + + def _pdf(self,x,mean,covar): + ''' + Gaussian PDF function + ''' + covar_det = numpy.linalg.det(covar); + + c = (1 / ( (2.0*numpy.pi)**(float(self.d/2.0)) * (covar_det)**(0.5))) + pdfval = c * numpy.exp(-0.5 * numpy.dot( numpy.dot((x-mean),covar.I), (x-mean)) ) + return pdfval \ No newline at end of file diff --git a/hmm/continuous/_ContinuousHMM.py b/hmm/continuous/_ContinuousHMM.py new file mode 100644 index 0000000..9b9fee8 --- /dev/null +++ b/hmm/continuous/_ContinuousHMM.py @@ -0,0 +1,175 @@ +''' +Created on Nov 12, 2012 + +@author: GuyZ +''' + +from hmm._BaseHMM import _BaseHMM +import numpy + +class _ContinuousHMM(_BaseHMM): + ''' + classdocs + ''' + + def __init__(self,n,m,d=1,A=None,means=None,covars=None,w=None,pi=None,min_std=0.01,init_type='uniform',precision=numpy.double,verbose=False): + ''' + Constructor + ''' + _BaseHMM.__init__(self,n,m,precision,verbose) #@UndefinedVariable + + self.d = d + self.A = A + self.pi = pi + self.means = means + self.covars = covars + self.w = w + self.min_std = min_std + + self.reset(init_type=init_type) + + def reset(self,init_type='uniform'): + if init_type == 'uniform': + self.pi = numpy.ones( (self.n), dtype=self.precision) *(1.0/self.n) + self.A = numpy.ones( (self.n,self.n), dtype=self.precision)*(1.0/self.n) + self.w = numpy.ones( (self.n,self.m), dtype=self.precision)*(1.0/self.m) + self.means = numpy.zeros( (self.n,self.m,self.d), dtype=self.precision) + self.covars = [[ numpy.matrix(numpy.ones((self.d,self.d), dtype=self.precision)) for j in xrange(self.m)] for i in xrange(self.n)] + elif init_type == 'user': + covars_tmp = [[ numpy.matrix(numpy.ones((self.d,self.d), dtype=self.precision)) for j in xrange(self.m)] for i in xrange(self.n)] + for i in xrange(self.n): + for j in xrange(self.m): + if type(self.covars[i][j]) is numpy.ndarray: + covars_tmp[i][j] = numpy.matrix(self.covars[i][j]) + else: + covars_tmp[i][j] = self.covars[i][j] + self.covars = covars_tmp + + """ + b[j][Ot] = sum(1...M)w[j][m]*b[j][m][Ot] + Returns b[j][Ot] based on the current model parameters (means, covars, weights) for the mixtures. + - j - state + - Ot - the current observation + Note: there's no need to get the observation itself as it has been used for calculation before. + """ + def _mapB(self,observations): + self.B_map = numpy.zeros( (self.n,len(observations)), dtype=self.precision) + self.Bmix_map = numpy.zeros( (self.n,self.m,len(observations)), dtype=self.precision) + for j in xrange(self.n): + for t in xrange(len(observations)): + self.B_map[j][t] = self._calcbjt(j, t, observations[t]) + + """ + b[j][Ot] = sum(1...M)w[j][m]*b[j][m][Ot] + Returns b[j][Ot] based on the current model parameters (means, covars, weights) for the mixtures. + - j - state + - Ot - the current observation + Note: there's no need to get the observation itself as it has been used for calculation before. + """ + def _calcbjt(self,j,t,Ot): + bjt = 0 + for m in xrange(self.m): + self.Bmix_map[j][m][t] = self._pdf(Ot, self.means[j][m], self.covars[j][m]) + bjt += (self.w[j][m]*self.Bmix_map[j][m][t]) + return bjt + + """ + Calculates 'gamma_mix' from xi. + + Gamma is a (TxNxK) numpy array, where gamma[t][i][m] = the probability of being + in state 'i' at time 't' with mixture 'm' given the full observation sequence. + """ + def _calcgammamix(self,alpha,beta,observations): + gamma_mix = numpy.zeros((len(observations),self.n,self.m),dtype=self.precision) + + for t in xrange(len(observations)): + for j in xrange(self.n): + for m in xrange(self.m): + alphabeta = 0.0 + for jj in xrange(self.n): + alphabeta += alpha[t][jj]*beta[t][jj] + comp1 = (alpha[t][j]*beta[t][j]) / alphabeta + + bjk_sum = 0.0 + for k in xrange(self.m): + bjk_sum += (self.w[j][k]*self.Bmix_map[j][k][t]) + comp2 = (self.w[j][m]*self.Bmix_map[j][m][t])/bjk_sum + + gamma_mix[t][j][m] = comp1*comp2 + + return gamma_mix + + def _updatemodel(self,new_model): + _BaseHMM._updatemodel(self,new_model) #@UndefinedVariable + + self.w = new_model['w'] + self.means = new_model['means'] + self.covars = new_model['covars'] + + def _calcstats(self,observations): + stats = _BaseHMM._calcstats(self,observations) #@UndefinedVariable + stats['gamma_mix'] = self._calcgammamix(stats['alpha'],stats['beta'],observations) + + return stats + + def _reestimate(self,stats,observations): + # re-estimate A, pi + new_model = _BaseHMM._reestimate(self,stats,observations) #@UndefinedVariable + + # re-estimate the continuous probability parameters of the mixtures + w_new, means_new, covars_new = self._reestimateMixtures(observations,stats['gamma_mix']) + + new_model['w'] = w_new + new_model['means'] = means_new + new_model['covars'] = covars_new + + return new_model + + def _reestimateMixtures(self,observations,gamma_mix): + w_new = numpy.zeros( (self.n,self.m), dtype=self.precision) + means_new = numpy.zeros( (self.n,self.m,self.d), dtype=self.precision) + covars_new = [[ numpy.matrix(numpy.zeros((self.d,self.d), dtype=self.precision)) for j in xrange(self.m)] for i in xrange(self.n)] + + for j in xrange(self.n): + for m in xrange(self.m): + numer = 0.0 + denom = 0.0 + for t in xrange(len(observations)): + for k in xrange(self.m): + denom += (self.eta(t,len(observations)-1)*gamma_mix[t][j][k]) + numer += (self.eta(t,len(observations)-1)*gamma_mix[t][j][m]) + w_new[j][m] = numer/denom + w_new[j] = self._normalize(w_new[j]) + + for j in xrange(self.n): + for m in xrange(self.m): + numer = numpy.zeros( (self.d), dtype=self.precision) + denom = numpy.zeros( (self.d), dtype=self.precision) + for t in xrange(len(observations)): + numer += (self.eta(t,len(observations)-1)*gamma_mix[t][j][m]*observations[t]) + denom += (self.eta(t,len(observations)-1)*gamma_mix[t][j][m]) + means_new[j][m] = numer/denom + + cov_prior = [[ numpy.matrix(self.min_std*numpy.eye((self.d), dtype=self.precision)) for j in xrange(self.m)] for i in xrange(self.n)] + for j in xrange(self.n): + for m in xrange(self.m): + numer = numpy.matrix(numpy.zeros( (self.d,self.d), dtype=self.precision)) + denom = numpy.matrix(numpy.zeros( (self.d,self.d), dtype=self.precision)) + for t in xrange(len(observations)): + vector_as_mat = numpy.matrix( (observations[t]-self.means[j][m]), dtype=self.precision ) + numer += (self.eta(t,len(observations)-1)*gamma_mix[t][j][m]*numpy.dot( vector_as_mat.T, vector_as_mat)) + denom += (self.eta(t,len(observations)-1)*gamma_mix[t][j][m]) + covars_new[j][m] = numer/denom + covars_new[j][m] = covars_new[j][m] + cov_prior[j][m] + + return w_new, means_new, covars_new + + def _normalize(self, arr): + summ = numpy.sum(arr) + for i in xrange(len(arr)): + arr[i] = (arr[i]/summ) + return arr + + def _pdf(self,x,mean,covar): + raise NotImplementedError("PDF function must be implemented") + \ No newline at end of file diff --git a/hmm/continuous/__init__.py b/hmm/continuous/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/hmm/discrete/DiscreteHMM.py b/hmm/discrete/DiscreteHMM.py new file mode 100644 index 0000000..c38e5ab --- /dev/null +++ b/hmm/discrete/DiscreteHMM.py @@ -0,0 +1,69 @@ +''' +Created on Nov 12, 2012 + +@author: GuyZ +''' + +from hmm._BaseHMM import _BaseHMM +import numpy + +class DiscreteHMM(_BaseHMM): + ''' + classdocs + ''' + + def __init__(self,n,m,A=None,B=None,pi=None,init_type='uniform',precision=numpy.double,verbose=False): + ''' + Constructor + ''' + _BaseHMM.__init__(self,n,m,precision,verbose) #@UndefinedVariable + + self.A = A + self.pi = pi + self.B = B + + self.reset(init_type=init_type) + + def reset(self,init_type='uniform'): + if init_type == 'uniform': + self.pi = numpy.ones( (self.n), dtype=self.precision) *(1.0/self.n) + self.A = numpy.ones( (self.n,self.n), dtype=self.precision)*(1.0/self.n) + self.B = numpy.ones( (self.n,self.m), dtype=self.precision)*(1.0/self.m) + + def _mapB(self,observations): + self.B_map = numpy.zeros( (self.n,len(observations)), dtype=self.precision) + + for j in xrange(self.n): + for t in xrange(len(observations)): + self.B_map[j][t] = self.B[j][observations[t]] + + def _updatemodel(self,new_model): + _BaseHMM._updatemodel(self,new_model) #@UndefinedVariable + + self.B = new_model['B'] + + def _reestimate(self,stats,observations): + # re-estimate A, pi + new_model = _BaseHMM._reestimate(self,stats,observations) #@UndefinedVariable + + # re-estimate the discrete probability of the observable symbols + B_new = self._reestimateB(observations,stats['gamma']) + + new_model['B'] = B_new + + return new_model + + def _reestimateB(self,observations,gamma): + B_new = numpy.zeros( (self.n,self.m) ,dtype=self.precision) + + for j in xrange(self.n): + for k in xrange(self.m): + numer = 0.0 + denom = 0.0 + for t in xrange(len(observations)): + if observations[t] == k: + numer += gamma[t][j] + denom += gamma[t][j] + B_new[j][k] = numer/denom + + return B_new \ No newline at end of file diff --git a/hmm/discrete/__init__.py b/hmm/discrete/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/hmm/examples/__init__.py b/hmm/examples/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/hmm/examples/tests.py b/hmm/examples/tests.py new file mode 100644 index 0000000..5b42ffa --- /dev/null +++ b/hmm/examples/tests.py @@ -0,0 +1,109 @@ +''' +Created on Nov 13, 2012 + +@author: GuyZ +''' + +from hmm.continuous.GMHMM import GMHMM +from hmm.discrete.DiscreteHMM import DiscreteHMM +import numpy + +def test_simple(): + n = 2 + m = 2 + d = 2 + pi = numpy.array([0.5, 0.5]) + A = numpy.ones((n,n),dtype=numpy.double)/float(n) + + w = numpy.ones((n,m),dtype=numpy.double) + means = numpy.ones((n,m,d),dtype=numpy.double) + covars = [[ numpy.matrix(numpy.eye(d,d)) for j in xrange(m)] for i in xrange(n)] + + w[0][0] = 0.5 + w[0][1] = 0.5 + w[1][0] = 0.5 + w[1][1] = 0.5 + means[0][0][0] = 0.5 + means[0][0][1] = 0.5 + means[0][1][0] = 0.5 + means[0][1][1] = 0.5 + means[1][0][0] = 0.5 + means[1][0][1] = 0.5 + means[1][1][0] = 0.5 + means[1][1][1] = 0.5 + + gmmhmm = GMHMM(n,m,d,A,means,covars,w,pi,init_type='user',verbose=True) + + obs = numpy.array([ [0.3,0.3], [0.1,0.1], [0.2,0.2]]) + + print "Doing Baum-welch" + gmmhmm.train(obs,10) + print + print "Pi",gmmhmm.pi + print "A",gmmhmm.A + print "weights", gmmhmm.w + print "means", gmmhmm.means + print "covars", gmmhmm.covars + +def test_rand(): + n = 5 + m = 4 + d = 2 + atmp = numpy.random.random_sample((n, n)) + row_sums = atmp.sum(axis=1) + a = numpy.array(atmp / row_sums[:, numpy.newaxis], dtype=numpy.double) + + wtmp = numpy.random.random_sample((n, m)) + row_sums = wtmp.sum(axis=1) + w = numpy.array(wtmp / row_sums[:, numpy.newaxis], dtype=numpy.double) + + means = numpy.array((0.6 * numpy.random.random_sample((n, m, d)) - 0.3), dtype=numpy.double) + covars = numpy.zeros( (n,m,d,d) ) + + for i in xrange(n): + for j in xrange(m): + for k in xrange(d): + covars[i][j][k][k] = 1 + + pitmp = numpy.random.random_sample((n)) + pi = numpy.array(pitmp / sum(pitmp), dtype=numpy.double) + + gmmhmm = GMHMM(n,m,d,a,means,covars,w,pi,init_type='user',verbose=True) + + obs = numpy.array((0.6 * numpy.random.random_sample((40,d)) - 0.3), dtype=numpy.double) + + print "Doing Baum-welch" + gmmhmm.train(obs,1000) + print + print "Pi",gmmhmm.pi + print "A",gmmhmm.A + print "weights", gmmhmm.w + print "means", gmmhmm.means + print "covars", gmmhmm.covars + +def test_discrete(): + + ob5 = (3,1,2,1,0,1,2,3,1,2,0,0,0,1,1,2,1,3,0) + print "Doing Baum-welch" + + atmp = numpy.random.random_sample((4, 4)) + row_sums = atmp.sum(axis=1) + a = atmp / row_sums[:, numpy.newaxis] + + btmp = numpy.random.random_sample((4, 4)) + row_sums = btmp.sum(axis=1) + b = btmp / row_sums[:, numpy.newaxis] + + pitmp = numpy.random.random_sample((4)) + pi = pitmp / sum(pitmp) + + hmm2 = DiscreteHMM(4,4,a,b,pi,init_type='user',precision=numpy.longdouble,verbose=True) + hmm2.train(numpy.array(ob5*10),100) + print "Pi",hmm2.pi + print "A",hmm2.A + print "B", hmm2.B + + +#test_simple() +test_rand() +#test_discrete() \ No newline at end of file