I have one HMM model for multiple sequences. Each sequence has different length. so, I built pm in the following way:

```
all_seq # the observed data sequences
with pm.Model() as model:
P = pm.Dirichlet('P_transition', a=np.ones((N_states,N_states)), shape=(N_states,N_states))
Pe = pm.Dirichlet('P_emission', a=np.ones((N_states,N_labels)), shape=(N_states,N_labels))
....
pos_states_seqs = []
pos_emission_seqs = []
for i in range(sequence_length):
all_seq_i = all_seq[i]
states_seq_i = HMMStates(...)
emission_seq_i = HMMMultinomialEmissions(..., observed=all_seq_i)
pos_states_seqs += [states_seq_i]
pos_emission_seqs += [states_seq_i]
print("find MAP")
start = pm.find_MAP()
print("run EM steps")
step1 = pm.Metropolis(vars=[P, Pe]+pos_emission_seqs)
step2 = pm.CategoricalGibbsMetropolis(vars=pos_states_seqs)
trace = pm.sample(10000, start=start, step=[step1,step2])
```

My questions are:

- how observed data`s log-likelihood are computed in this looping way? add each sequence together?
- if you know any better inferences for this HMM problem?