I wonder if anyone can point me to an example of Baum-Welch algorithm in pymc3 that can be used to learn the parameters of a given Hidden Markov Model.
I have trying to use python hmmlearn to do that but doesnt seem to be very good. The problem I am having from the following code is that I do not get the same parameter that I used with my toy example to generate the sequences.
The following is toy example in which first I initiate a hmm, then I generate some sequences using known parameters, finally I initiate a second hmm and need to train (fit) that with my generated sequences to learn the initial parameters. But the issue is never get close enough to my original parameter, even for this simple example, that makes me worried how I am going to use this for my complex data.
I therefore wonder if there is anything similar in pymc3, as I think stan has a fairly good coverage of HMMs.
def initialize_components(n_components, n_emissions): """ Initialize start probabilities, transition and emission matrices, start for a model with a fixed number of components, for Multinomial model with a certain number of dimensions. """ trans_mat = np.abs(np.random.randn(n_components, n_components)) trans_mat = (trans_mat.T / trans_mat.sum( axis=1 )).T emission_mat = np.abs(np.random.randn(n_components, n_emissions)) emission_mat = (emission_mat.T/emission_mat.sum(axis=1)).T # print(emission_mat.shape) start_probs = np.abs( np.random.randn(n_components) ) start_probs /= start_probs.sum() return (start_probs, trans_mat, emission_mat) def generate_sequences(Pi, T, E, n_components=3, n_features=6, N=200, K=20): """ Use the given parameters, initiate a hmm model and generate some sequences """ states =  seqs =  seq_lengths =  model = hmm.MultinomialHMM(n_components=3) model.startprob_ = Pi model.transmat_ = T model.emissionprob_ = E model.n_features = n_features print('*** Original parameters *********') print('p_init') print(model.startprob_) print('p_tans ') print(model.transmat_) print('p_emm') print(model.emissionprob_) for j in range(N): seq_length = K ### generate sequence seq = model.sample(seq_length) z = model.sample(seq_length) # print(z) # print(seq) ### generate sequence seq = model.sample(seq_length) seq_lengths.append(seq_length) states.append(z) # print(seq) if len(seqs) > 0: seqs = np.append(seqs, seq, axis=0) else: seqs = seq model.fit(seqs, seq_lengths) score = model.score(seqs) print('Score: ', score) return(seqs, states, seq_lengths) def train(hmm_model, Seqs, Lengths, nIteration): """ Train the model with given data, and find the best describing parameters (Baub-Welch) """ best_model = None best_trans = None best_emissions = None best_pi = None max_score = -np.inf for i in range(nIteration): pi, A, B = initialize_components(n_states, n_emissions) hmm_model.emissionprob_ = B hmm_model.transmat_ = A hmm_model.startprob_ = pi hmm_model.n_features_ = 6 hmm_model.fit(Seqs, Lengths) this_score = hmm_model.score(Seqs) if this_score > max_score: max_score = this_score best_model = hmm_model best_pi = pi best_trans = A best_emissions = B ### used these start poing if i % 10 == 0: print('^'*50) print('\niteration', i) print('best score: ', max_score) # print(i , '\n ================== After training ==============') # print(model2.score(seqs)) print(max_score) print(best_model.startprob_) print(best_model.transmat_) print(best_model.emissionprob_) return(best_model) ################## S = ["S1", "S2", 'S3'] n_states = len(S) E = [i+1 for i in range(6)] ## sides of a dice trans_prob = np.array( [[0.9, 0.05, 0.05], [0.1, 0.8, 0.1], [0.05, 0.05, 0.9]] ) start_prob = np.array([0.8, 0.1, 0.1]) ### Lets assume that we are throwing 3 6-sided dices: one is fair with equel probablity of each side, and the other two are loaded em_prob = np.array([[1/4, 3/20, 3/20, 3/20, 3/20, 3/20], [2/15, 2/15, 2/15, 2/15, 2/15, 1/3], [1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6]]) X, Z, L = generate_sequences(start_prob, trans_prob, em_prob, n_components=3, n_features=6, N=200, K=20) model2 = hmm.MultinomialHMM(n_components=3, n_iter=200) fitted_hmm = train(model2, X, L, 100)