How to infer the number of states in a Hidden Markov Model with Gaussian mixture emissions

I have a time series made up of an unknown number of hidden states. Each state contains a set of values unique to that state. I am trying to use a GMM HMM (as implemented in Python's hmmlearn package) to identify these hidden states (so I'm effectively clustering a time series). This seems to work reasonably well when I know the number of hidden states (K) to look for, but the whole point is that I don't know the number of states and need to find this out before applying the model.

I've been trying out a range of values for K and getting the log-likelihood for each. Then, I've been using np.polyfit to fit a curve to this data and find a maximum.

This all feels a bit clumsy though and I'm wondering if there is a better way to get the best K value (especially as I have to guess a sensible range to test K for which could lead to a lot of models being built)? I'm very new to HMMs in general so if anyone can point me in the direction of some good introductory guides that could help me understand this problem better then that would also be much appreciated.

• This is the classical clustering problem, if you dont know how many clusters you have or want you should look into stats.stackexchange.com/questions/241381/… – knk Dec 28 '18 at 19:30
• I don't think my problem here is understanding what the clustering problem is, but rather how to implement something with HMMs. For example I can see plenty of references to Hierarchical HMM clustering, but no information on how to implement this - what do you use as your linkage criteria? Having read further, I can see that iHMMs and their implementation in pomegranate might help me, but I'm not 100% sure yet... – big-o Dec 28 '18 at 23:21
• As an update to my last comment - the pomegranate python package does NOT implement iHMMs. They use the word "infinite" in their literature, but this is to describe HMMs with an infinitely long chain, rather than a true Infinite HMM, which is something different. – big-o Sep 25 '19 at 17:20

Have you tried using AIC (Akaike Information criterion) or BIC (Bayesian Information Criterion)?

It's possible to implement AIC or BIC to work with hmmlearn. Here is my implementation for GaussianHMM for covariance_type='diag'. If the covariance_type changes then the number of parameters will have to be adjusted for covars_. You can extend it to GMMHMM if you know the number for components of the GMM.

import hmmlearn as hmm
import numpy as np

def bic_general(likelihood_fn, k, X):
"""likelihood_fn: Function. Should take as input X and give out   the log likelihood
of the data under the fitted model.
k - int. Number of parameters in the model. The parameter that we are trying to optimize.
For HMM it is number of states.
For GMM the number of components.
X - array. Data that been fitted upon.
"""
bic = np.log(len(X))*k - 2*likelihood_fn(X)
return bic

def bic_hmmlearn(X):
lowest_bic = np.infty
bic = []
n_states_range = range(1,7)
for n_components in n_states_range:
hmm_curr = hmm.GaussianHMM(n_components=n_components, covariance_type='diag')
hmm_curr.fit(X)

# Calculate number of free parameters
# free_parameters = for_means + for_covars + for_transmat + for_startprob
# for_means & for_covars = n_features*n_components
n_features = hmm_curr.n_features
free_parameters = 2*(n_components*n_features) + n_components*(n_components-1) + (n_components-1)

bic_curr = bic_general(hmm_curr.score, free_parameters, X)
bic.append(bic_curr)
if bic_curr < lowest_bic:
lowest_bic = bic_curr
best_hmm = hmm_curr

return (best_hmm, bic)

best_hmm, bic = bic_hmmlearn(X)


Similar question. Number of parameters in Markov model

• Thanks - I think this is about as good as you can get with a classic HMM but unfortunately it doesn't get around the problem of computational expense from fitting multiple models. It's essentially the same approach as what I was doing before, but with a different metric. I think the only way to avoid this problem might be using a different type of model. – big-o Jun 14 '20 at 16:41

The AIC/BIC are good for this, as mentioned by another answer. Another alternative would be to treat number of hidden startes as a hyperparameter and do a hyperparameter search, optimizing your desired metric on a validation set.