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

## 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.

Subscribe to:
Post Comments (Atom)

## No comments:

## Post a Comment