## 20160708

### Numerically stable Viterbi algorithm in Python for hidden markov model state inference

The Python demonstration of the Viterbi algorithm on Wikipedia is numerically unstable for moderate to large problem sizes. This implementation is phrased in terms of log probabilities, and so has better numerical properties. Please reuse, creative commons.

def hmm_viterbi(Y,logP,logA,logB):
'''
See https://en.wikipedia.org/wiki/Viterbi_algorithm

Parameters
----------
Y : 1D array
Observations (integer states)
logP : array shape = (nStates ,)
1D array of priors for initial state
given in log probability
logA : array (nStates,nStates)
State transition matrix given in log probability
logB : ndarray K x N
conditional probability matrix
log probabilty of each observation given each state
'''
K = len(logP)         # Number of states
T = len(Y)            # Number of observations
N = np.shape(logB)[1] # Number of states
Y = np.int32(Y)

assert np.shape(logA)==(K,K)
assert np.shape(logB)==(K,N)

# The initial guess for the first state is initialized as the
# probability of observing the first observation given said
# state, multiplied by the prior for that state.
logT1 = np.zeros((K,T),'float') # Store probability of most likely path
logT1[:,0] = logP + logB[:,Y[0]]

# Store estimated most likely path
T2 = np.zeros((K,T),'float')

# iterate over all observations from left to right
for i in range(1,T):
# iterate over states 1..K (or 0..K-1 with zero-indexing)
for s in range(K):
# The likelihood of a new state is the likelihood of
# transitioning from either of the previous states.
# We incorporate a multiplication by the prior here
log_filtered_likelihood = logT1[:,i-1] + logA[:,s] + logB[s,Y[i]]
best = np.argmax(log_filtered_likelihood)
logT1[s,i] = log_filtered_likelihood[best]
# We save which state was the most likely
T2[s,i] = best

# At the end, choose the most likely state, then
# Iterate backwards over the data and fill in the state estimate
X     = np.zeros((T,) ,'int'  ) # Store our inferred hidden states
X[-1] = np.argmax(logT1[:,-1])
for i in range(1,T)[::-1]:
X[i-1] = T2[X[i],i]
return X