Input Output Hidden Markov Model Implementation in Python

Viewed 863

I am trying to implement Hidden Markov Models with Input Output Architecture but I could not find any good python implementation for the same.

Can anybody share the Python package the would consider the following implementation for HMM.

Allow continuous emissions. Allow functionality of covariates(i.e. Independent Variables in I/O HMM).

At this moment, I am struggling to find the python implementation for the same.

I could not find the relevant examples in hmmlearn.

Here are few of the libraries that I have tested out:

hmmlearn: hmmlearn allows to pass multiple features to emissions/observations but does not provide the support to include co-variates(i.e. Independent Variables).

hmms: Does not support the functionality to add continuous emissions as well as does not support addition of Independent Variables.

IOHMM: I was able to train the HMM Model using this library, but could not find the documentation to make predictions after training the model.

Therefore, I am looking for the package that serves the purpose.

from IOHMM import UnSupervisedIOHMM
from IOHMM import OLS, DiscreteMNL, CrossEntropyMNL, forward_backward

SHMM = UnSupervisedIOHMM(num_states=3, max_EM_iter=200, EM_tol=1e-6)

SHMM.set_models(model_emissions = [OLS(est_stderr=True)], 
                model_transition=CrossEntropyMNL(solver='lbfgs'),
                model_initial=CrossEntropyMNL(solver='lbfgs'))

SHMM.set_inputs(covariates_initial = [], covariates_transition = [], covariates_emissions = [['Insulin']])


SHMM.set_outputs([['Glucose']])

SHMM.set_data([data])

SHMM.train() 

I could not figure out how to get emission probabilities and sequence of hidden states after the above training.

1 Answers

Referring to "https://web.ece.ucsb.edu/Faculty/Rabiner/ece259/Reprints/tutorial%20on%20hmm%20and%20applications.pdf" and the library "https://hmmlearn.readthedocs.io/en/latest/" I have found this solution:

1- Through log_gamma (posterior distribution):

state_sequences = []
for i in range(100):
    for j in range(lengths[i]):
        state_sequences.append(np.argmax(np.exp(SHMM.log_gammas[i])[j]))
pred_state_seq = [state_sequences[df[df['unit'] == i].index[0]:df[df['unit'] == i].index[-1] + 1] for i in
                  range(1, df_A['unit'].max() + 1)]

2- Viterbi Algorithm:

from hmmlearn import _hmmc

transmat = np.empty((num_states, num_states))
for i in range(num_states):
    transmat = np.concatenate((transmat, np.exp(SHMM.model_transition[i].predict_log_proba(np.array([[]])))))
transmat = transmat[num_states:]

startprob = np.exp(SHMM.model_initial.predict_log_proba(np.array([[]]))).squeeze()


def log_mask_zero(a):
    """
    Compute the log of input probabilities masking divide by zero in log.

    Notes
    -----
    During the M-step of EM-algorithm, very small intermediate start
    or transition probabilities could be normalized to zero, causing a
    *RuntimeWarning: divide by zero encountered in log*.

    This function masks this unharmful warning.
    """
    a = np.asarray(a)
    with np.errstate(divide="ignore"):
        return np.log(a)


def _do_viterbi_pass(framelogprob):
    n_samples, n_components = framelogprob.shape
    state_sequence, logprob = _hmmc._viterbi(n_samples, n_components, log_mask_zero(startprob),
                                             log_mask_zero(transmat), framelogprob)
    return logprob, state_sequence


def _decode_viterbi(X):
    framelogprob = SHMM.log_Eys[X]
    return _do_viterbi_pass(framelogprob)


def decode():
    decoder = {"viterbi": _decode_viterbi}["viterbi"]
    logprob = 0
    sub_state_sequences = []
    for sub_X in range(100):
        # XXX decoder works on a single sample at a time!
        sub_logprob, sub_state_sequence = decoder(sub_X)
        logprob += sub_logprob
        sub_state_sequences.append(sub_state_sequence)
    return logprob, np.concatenate(sub_state_sequences)


def predict():
    """
    Find most likely state sequence corresponding to ``X``.

    Parameters
    ----------
    X : array-like, shape (n_samples, n_features)
        Feature matrix of individual samples.
    lengths : array-like of integers, shape (n_sequences, ), optional
        Lengths of the individual sequences in ``X``. The sum of
        these should be ``n_samples``.

    Returns
    -------
    state_sequence : array, shape (n_samples, )
        Labels for each sample from ``X``.
    """
    logprob, state_sequence = decode()
    return logprob, state_sequence


_, state_seq = predict()

pred_state_seq = [state_seq[df[df['unit'] == i].index[0]:df[df['unit'] == i].index[-1] + 1] for i in
                  range(1, df_A['unit'].max() + 1)]
Related