HMM with multiple sequences

[jcatulo]

Hi!

I am new to bayesian modeling and PyMC3 and I am trying to estimate the parameters of a Hidden Markov Model (HMM) with gaussian emissions for multiple short sequences.

I followed the examples on the examples directory and I was able to correctly implement a HMM with only one sequence of observations. However, in my problem, I have multiple sequences of observations and I want to estimate the parameters of the HMM with 3 hidden states. I already tried several approaches but I was not successful on any of them.

I tried to generate, for each sequence, a DiscreteMarkovChain and aSwitchingProcess random variable in a for loop, that iterated over each sequence of observations, but the code was very slow and I could not obtain any results. I saw that theano.scan allows for more efficient iterations, but I was not able to implement it either.

Here is my code:

N_states = 3
    
with pm.Model() as model:
    Pt = pm.Dirichlet("p_transition", np.ones( (N_states, N_states) ), shape=(N_states, N_states))
    P0 = pm.Dirichlet("p_init", np.ones((N_states,)), shape=(N_states,))

    mu = pm.Normal('mu', mu=-15, sigma=2, shape = (N_states,) )
    sigma = pm.HalfStudentT("sigma", sigma=2, nu=10, shape=(N_states,) )

    comp_dists = [ pm.TruncatedNormal.dist(mu = mu[i], sigma = sigma[i], lower=-30, upper=0) for i in range(0, N_states)]

    Z_rv = []
    X_rv = []

    for idx, i in enumerate(obs):
        z_i = DiscreteMarkovChain("Z_t"+str(idx), tt.shape_padleft(Pt), P0, shape= np.shape(i)[-1] )
        x_i = SwitchingProcess("X_t"+str(idx), comp_dists, z_i, observed= i)

        Z_rv += [z_i]
        X_rv += [x_i]

How can I estimate the parameters of a Hidden Markov Model with multiple sequences?
I am sorry if this is a very basic question, but I am new to this and I am really trying to learn!

Thank you in advance!

[brandonwillard]

It looks like obs is missing from your example. We'll likely need an MWE (with no external dependencies, like data files) and the resulting output in order to properly understand the situation.

Aside from that, if you want to estimate multiple series at once, the usual answer is to expand the dimensions of the DiscreteMarkovChain and SwitchingProcess distributions in your model; however, DiscreteMarkovChain's logp implementation does not support that (see here). PRs that remove this restriction are, of course, welcome.
Ultimately, this means that you will need to create distinct DiscreteMarkovChain instances if you want multiple state sequences, as your example is currently doing.

You might be able to combine those distinct state sequences (e.g. using tt.stack) and use them as a single input to a SwitchingProcess, but it's possible that the FFBSStep step method might not like that. There are some potential computational advantages to using extra dimensions like this, so that's why I mention it.

Regardless, you can always create distinct DiscreteMarkovChain and SwitchingProcess instances for each observed series—like you've done in the example—and get effectively the same result.

I tried to generate, for each sequence, a DiscreteMarkovChain and aSwitchingProcess random variable in a for loop, that iterated over each sequence of observations, but the code was very slow and I could not obtain any results. I saw that theano.scan allows for more efficient iterations, but I was not able to implement it either.

Without an explicit value of obs, I can't say anything about why it would've been slow. For instance, if obs is a vector of N-many observations, then it looks like your code would create N-many separate HMMs, and that probably isn't what you'd want—and could explain the latency if N is large enough.

[jcatulo]

Thank you for your reply.

In my problem I have 8000 sequences. Each of the sequences has different lengths, which can also cause problems. In order to start from a simpler problem, I am only taking 10 of those sequences where each of them has the same length (20).

You might be able to combine those distinct state sequences (e.g. using tt.stack) and use them as a single input to a SwitchingProcess, but it's possible that the FFBSStep step method might not like that. There are some potential computational advantages to using extra dimensions like this, so that's why I mention it.

I tried to combine distinct state sequences and use them as a single input to a SwitchingProcess, however, I cannot implement the FFBSStep step method, because this method only accepts one variable as input (here). Is there any way to handle this?

Here is the updated code:

N_states = 3
with pm.Model() as model:
    Pt = pm.Dirichlet("p_transition", np.ones( (N_states, N_states) ), shape=(N_states, N_states))
    P0 = pm.Dirichlet("p_init", np.ones((N_states,)), shape=(N_states,))

    mu = pm.Normal('mu', mu=-30, sigma=2, shape = (N_states,) )
    sigma = pm.HalfStudentT("sigma", sigma=2, nu=10, shape=(N_states,) )

    comp_dists = [ pm.TruncatedNormal.dist(mu = mu[i], sigma = sigma[i], lower=-30, upper=0) for i in range(0, N_states)]

    Z_rv = []
    X_rv = []

    for idx, sequence in enumerate( obs[0:10] ):
        z_i = DiscreteMarkovChain("Z_t"+str(idx), tt.shape_padleft(Pt), P0, shape= len(sequence) )
        Z_rv += [z_i]

    Z_rv = tt.stack(Z_rv)
    X = SwitchingProcess("Y_t", comp_dists, Z_rv, observed= obs[0:10])

Additionally, is it possible to expand this problem by estimating the parameters of HMM with multiple sequences with different lengths?

[brandonwillard]

Which samplers does PyMC say you're using when you perform those slow runs? My guess is that the FFBS sampler isn't happy with the fact that you have multiple observed distributions, so it's defaulting to the NUTS sampler for everything, and that can easily take forever.

[jcatulo]

I think I managed to estimate the parameters of the HMM with multiple sequences. I removed the step input from the function pm.sample, in order for the step methods to be assigned automatically. I read here that Discrete variables are assigned to Metropolis and Continuous variables to NUTS. In fact, PyMC says that I am using Metropolis for the DiscreteMarkovChain random variables and NUTS for the rest.

However, I would like to expand my problem to sequences of multiple lengths. I created distinct DiscreteMarkovChain instances for each sequence and combined them using np.stack. Then, I tried to use the combined state sequences as single input to a SwitchingProcess, but I get the following error: ValueError: setting an array element with a sequence. What am I doing wrong?

How can I estimate the parameters of the HMM with multiple sequences of different lengths?

Edit: I analyzed the code from the class SwitchingProcess with more attention and I concluded that my approach was wrong. The functions inside the class SwitchingProcess transform both the states and observations sequences into Theano Tensors (here and here), so the sequences must have the same length.
I cannot combine the distinct DiscreteMarkovChain instances for each sequence and use them as single input of a SwitchingProcess. I think that the only solution is to create distinct DiscreteMarkovChain and SwitchingProcess instances for each observed series. Is there any other valid approach?

[brandonwillard]

Yes, you will need to create distinct DiscreteMarkovChain and SwitchingProcess for every series, but that's mostly because the math/model itself likely won't make sense otherwise.