Python的隐马尔科夫HMMLearn库的应用教学

1、Python HMMLearn TutorialEdited By 毛片物语hmmlearnimplements the Hidden Markov Models (HMMs). The HMM is a generative probabilistic model, in which a sequence of observable(mathbfX)variables is generated by a sequence of internal hidden states(mathbfZ). The hidden states are not be observed directly. Th

2、e transitions between hidden states are assumed to have the form of a (first-order) Markov chain. They can be specified by the start probability vector(boldsymbolpi)and a transition probability matrix(mathbfA). The emission probability of an observable can be any distribution with parameters(boldsym

3、boltheta)conditioned on the current hidden state. The HMM is completely determined by(boldsymbolpi),(mathbfA)and(boldsymboltheta).There are three fundamental problems for HMMs: Given the model parameters and observed data, estimate the optimal sequence of hidden states. Given the model parameters an

4、d observed data, calculate the likelihood of the data. Given just the observed data, estimate the model parameters.The first and the second problem can be solved by the dynamic programming algorithms known as the Viterbi algorithm and the Forward-Backward algorithm, respectively. The last one can be

5、 solved by an iterative Expectation-Maximization (EM) algorithm, known as the Baum-Welch algorithm.References:Rabiner89Lawrence R. Rabiner “A tutorial on hidden Markov models and selected applications in speech recognition”, Proceedings of the IEEE 77.2, pp. 257-286, 1989.Bilmes98Jeff A. Bilmes, “A

6、gentle tutorial of the EM algorithm and its application to parameter estimation for Gaussian mixture and hidden Markov models.”, 1998.Available modelshmm.GaussianHMMHidden Markov Model with Gaussian emissions.hmm.GMMHMMHidden Markov Model with Gaussian mixture emissions.hmm.MultinomialHMMHidden Mark

7、ov Model with multinomial (discrete) emissionsRead onfor details on how to implement an HMM with a custom emission probability.Building HMM and generating samplesYou can build an HMM instance by passing the parameters described above to the constructor. Then, you can generate samples from the HMM by

8、 callingsample. import numpy as np from hmmlearn import hmm np.random.seed(42) model = hmm.GaussianHMM(n_components=3, covariance_type=full) model.startprob_ = np.array(0.6, 0.3, 0.1) model.transmat_ = np.array(0.7, 0.2, 0.1,. 0.3, 0.5, 0.2,. 0.3, 0.3, 0.4) model.means_ = np.array(0.0, 0.0, 3.0, -3.

9、0, 5.0, 10.0) model.covars_ = np.tile(np.identity(2), (3, 1, 1) X, Z = model.sample(100)The transition probability matrix need not to be ergodic. For instance, a left-right HMM can be defined as follows: lr = hmm.GaussianHMM(n_components=3, covariance_type=diag,. init_params=cm, params=cmt) lr.start

10、prob_ = np.array(1.0, 0.0, 0.0) lr.transmat_ = np.array(0.5, 0.5, 0.0,. 0.0, 0.5, 0.5,. 0.0, 0.0, 1.0)If any of the required parameters are missing,samplewill raise an exception: hmm.GaussianHMM(n_components=3) X, Z = model.sample(100)Traceback (most recent call last): .sklearn.utils.validation.NotF

11、ittedError: This GaussianHMM instance is not fitted yet. Call fit with appropriate arguments before using this method.Fixing parametersEach HMM parameter has a character code which can be used to customize its initialization and estimation. EM algorithm needs a starting point to proceed, thus prior

12、to training each parameter is assigned a value either random or computed from the data. It is possible to hook into this process and provide a starting point explicitly. To do so1. ensure that the character code for the parameter is missing frominit_paramsand then2. set the parameter to the desired

13、value.For example, consider an HMM with explicitly initialized transition probability matrix model = hmm.GaussianHMM(n_components=3, n_iter=100, init_params=mcs) model.transmat_ = np.array(0.7, 0.2, 0.1,. 0.3, 0.5, 0.2,. 0.3, 0.3, 0.4)A similar trick applies to parameter estimation. If you want to f

14、ix some parameter at a specific value, remove the corresponding character fromparamsand set the parameter value before training.Examples: Sampling from HMMTraining HMM parameters and inferring the hidden statesYou can train an HMM by calling thefitmethod. The input is a matrix of concatenated sequen

15、ces of observations (akasamples) along with the lengths of the sequences (seeWorking with multiple sequences).Note, since the EM algorithm is a gradient-based optimization method, it will generally get stuck in local optima. You should in general try to runfitwith various initializations and select

16、the highest scored model.The score of the model can be calculated by thescoremethod.The inferred optimal hidden states can be obtained by callingpredictmethod. Thepredictmethod can be specified with decoder algorithm. Currently the Viterbi algorithm (viterbi), and maximum a posteriori estimation (ma

17、p) are supported.This time, the input is a single sequence of observed values. Note, the states inremodelwill have a different order than those in the generating model. remodel = hmm.GaussianHMM(n_components=3, covariance_type=full, n_iter=100) remodel.fit(X) GaussianHMM(algorithm=viterbi,. Z2 = rem

18、odel.predict(X)Monitoring convergenceThe number of EM algorithm iteration is upper bounded by then_iterparameter. The training proceeds untiln_itersteps were performed or the change in score is lower than the specified thresholdtol. Note, that depending on the data EM algorithm may or may not achiev

19、e convergence in the given number of steps.You can use themonitor_attribute to diagnose convergence: remodel.monitor_ ConvergenceMonitor(history=., iter=12, n_iter=100, tol=0.01, verbose=False) remodel.monitor_.convergedTrueWorking with multiple sequencesAll of the examples so far were using a singl

20、e sequence of observations. The input format in the case of multiple sequences is a bit involved and is best understood by example.Consider two 1D sequences: X1 = 0.5, 1.0, -1.0, 0.42, 0.24 X2 = 2.4, 4.2, 0.5, -0.24To pass both sequences tofitorpredictfirst concatenate them into a single array and t

21、hen compute an array of sequence lengths: X = np.concatenate(X1, X2) lengths = len(X1), len(X2)Finally just call the desired method withXandlengths: hmm.GaussianHMM(n_components=3).fit(X, lengths) GaussianHMM(algorithm=viterbi, .Examples: Gaussian HMM of stock dataSaving and loading HMMAfter traning

22、 an HMM can be easily persisted for future use with the standardpicklemodule or its more efficient replacement in thejoblibpackage: from sklearn.externals import joblib joblib.dump(remodel, filename.pkl)filename.pkl joblib.load(filename.pkl) GaussianHMM(algorithm=viterbi,.Implementing HMMs with cust

23、om emission probabilitiesIf you want to implement other emission probability (e.g. Poisson), you have to subclass_BaseHMMand override the following methodsbase._BaseHMM._init(X,lengths)Initializes model parameters prior to fitting.base._BaseHMM._check()Validates model parameters prior to fitting.bas

24、e._BaseHMM._generate_sample_from_state(state)Generates a random sample from a given component.base._BaseHMM._compute_log_likelihood(X)Computes per-component log probability under the model.base._BaseHMM._initialize_sufficient_statistics()Initializes sufficient statistics required for M-step.base._Ba

25、seHMM._accumulate_sufficient_statistics(.)Updates sufficient statistics from a given sample.base._BaseHMM._do_mstep(stats)Performs the M-step of EM algorithm.Sampling from HMMThis script shows how to sample points from a Hiden Markov Model (HMM): we use a 4-components with specified mean and covaria

26、nce.The plot show the sequence of observations generated with the transitions between them. We can see that, as specified by our transition matrix, there are no transition between component 1 and 3.print(_doc_)import numpy as npimport matplotlib.pyplot as pltfrom hmmlearn import hmmPrepare parameter

27、s for a 4-components HMM Initial population probabilitystartprob = np.array(0.6, 0.3, 0.1, 0.0)# The transition matrix, note that there are no transitions possible# between component 1 and 3transmat = np.array(0.7, 0.2, 0.0, 0.1, 0.3, 0.5, 0.2, 0.0, 0.0, 0.3, 0.5, 0.2, 0.2, 0.0, 0.2, 0.6)# The means

28、 of each componentmeans = np.array(0.0, 0.0, 0.0, 11.0, 9.0, 10.0, 11.0, -1.0)# The covariance of each componentcovars = .5 * np.tile(np.identity(2), (4, 1, 1)# Build an HMM instance and set parametersmodel = hmm.GaussianHMM(n_components=4, covariance_type=full)# Instead of fitting it from the data,

29、 we directly set the estimated# parameters, the means and covariance of the componentsmodel.startprob_ = startprobmodel.transmat_ = transmatmodel.means_ = meansmodel.covars_ = covars# Generate samplesX, Z = model.sample(500)# Plot the sampled dataplt.plot(X:, 0, X:, 1, .-, label=observations, ms=6,

30、mfc=orange, alpha=0.7)# Indicate the component numbersfor i, m in enumerate(means): plt.text(m0, m1, Component %i % (i + 1), size=17, horizontalalignment=center, bbox=dict(alpha=.7, facecolor=w)plt.legend(loc=best)plt.show()Total running time of the script:(0 minutes 0.528 seconds)Gaussian HMM of st

31、ock dataThis script shows how to use Gaussian HMM on stock price data from Yahoo! finance. For more information on how to visualize stock prices with matplotlib, please refer todate_demo1.pyof matplotlib.from _future_ import print_functionimport datetimeimport numpy as npfrom matplotlib import cm, p

32、yplot as pltfrom matplotlib.dates import YearLocator, MonthLocatortry: from matplotlib.finance import quotes_historical_yahoo_ochlexcept ImportError: # For Matplotlib prior to 1.5. from matplotlib.finance import ( quotes_historical_yahoo as quotes_historical_yahoo_ochl )from hmmlearn.hmm import Gaus

33、sianHMMprint(_doc_)Get quotes from Yahoo! financequotes = quotes_historical_yahoo_ochl( INTC, datetime.date(1995, 1, 1), datetime.date(2012, 1, 6)# Unpack quotesdates = np.array(q0 for q in quotes, dtype=int)close_v = np.array(q2 for q in quotes)volume = np.array(q5 for q in quotes)1:# Take diff of

34、close value. Note that this makes# len(diff) = len(close_t) - 1, therefore, other quantities also# need to be shifted by 1.diff = np.diff(close_v)dates = dates1:close_v = close_v1:# Pack diff and volume for training.X = np.column_stack(diff, volume)Run Gaussian HMMprint(fitting to HMM and decoding ., end=)# Make an HMM instance and execute fitmodel = GaussianHMM(n_components=4, covariance_type=diag, n_iter=1000).fit(X)# Predict the optimal sequence of internal hidden statehidden_states = model.predict(X)print(done)Out:fitting to HMM and decoding .donePrint trained parameters a