Example of time-varying HMM transition matrix

Multinomial_Logit_HMM.py

import numpy as np

#A simple observation series
obs = np.array([0, 0, 1, 0, 0])

#initial state occupation probabilities
#(flat prior over state occupancies when t = 0)
initProbs = np.ones(4)/4

#These are the time-constant transition *rates* (with zeros
#representing non-constant transitions or true zero 
#rate transitions. The probability of remaining in the
#current state is 1-(p(all other transitions)), so this can 
#be left empty as well

#transition rates are specified as logits, with
#current state as the reference category
epsilon = -0.5
gamma = -0.5
tmat = np.array([[0.0, 0.0, 0.0, 0.0], 
  			[0.0, 0.0, epsilon, 0.0],
				[0.0, 0.0, 0.0, gamma],
				[0.0, 0.0, 0.0, 0.0]])

#Expand the transition matrix with a time-varying component,
#in this case on the (S -> E) transition

#To start, we'll just sample some random values
inf_rates = np.random.uniform(-1,1., len(obs)-1)

#and put them into the (0,1) spot on the  expanded transition 
#matrices
tv_tmat = np.array([tmat]*(len(obs)-1))
tv_tmat[:,0,1] = inf_rates

#now convert these to probabilities
pr_tm = []
for tm in tv_tmat:
	ptm = np.copy(tm)
	for i,row in enumerate(ptm):
		#If there are no transitions in the row
		#then it is an absorbing state; move on
		if np.sum(row) == 0:
			ptm[i,i] = 1.0
		else:
			#Exponentiate probabilities for all
			#off-diagonal transition rates > 0
			nz = row != 0
			all_pr = np.exp(row[nz])
			total_pr = 1 + np.sum(all_pr)
			ptm[i,nz] = all_pr / total_pr
			ptm[i,i] = 1. / total_pr
	pr_tm.append(ptm)

pr_tm = np.array(pr_tm)