matthiasm@8: function [path,p] = viterbi_path(prior, transmat, obslik)
matthiasm@8: % VITERBI Find the most-probable (Viterbi) path through the HMM state trellis.
matthiasm@8: % path = viterbi(prior, transmat, obslik)
matthiasm@8: %
matthiasm@8: % Inputs:
matthiasm@8: % prior(i) = Pr(Q(1) = i)
matthiasm@8: % transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i)
matthiasm@8: % obslik(i,t) = Pr(y(t) | Q(t)=i)
matthiasm@8: %
matthiasm@8: % Outputs:
matthiasm@8: % path(t) = q(t), where q1 ... qT is the argmax of the above expression.
matthiasm@8: 
matthiasm@8: 
matthiasm@8: % delta(j,t) = prob. of the best sequence of length t-1 and then going to state j, and O(1:t)
matthiasm@8: % psi(j,t) = the best predecessor state, given that we ended up in state j at t
matthiasm@8: 
matthiasm@8: scaled = 1;
matthiasm@8: 
matthiasm@8: T = size(obslik, 2);
matthiasm@8: prior = prior(:);
matthiasm@8: Q = length(prior);
matthiasm@8: 
matthiasm@8: delta = zeros(Q,T);
matthiasm@8: psi = zeros(Q,T);
matthiasm@8: path = zeros(1,T);
matthiasm@8: scale = ones(1,T);
matthiasm@8: 
matthiasm@8: 
matthiasm@8: t=1;
matthiasm@8: delta(:,t) = prior .* obslik(:,t);
matthiasm@8: if scaled
matthiasm@8:   [delta(:,t), n] = normalise(delta(:,t));
matthiasm@8:   scale(t) = 1/n;
matthiasm@8: end
matthiasm@8: psi(:,t) = 0; % arbitrary value, since there is no predecessor to t=1
matthiasm@8: for t=2:T
matthiasm@8:   for j=1:Q
matthiasm@8:     [delta(j,t), psi(j,t)] = max(delta(:,t-1) .* transmat(:,j));
matthiasm@8:     delta(j,t) = delta(j,t) * obslik(j,t);
matthiasm@8:   end
matthiasm@8:   if scaled
matthiasm@8:     [delta(:,t), n] = normalise(delta(:,t));
matthiasm@8:     scale(t) = 1/n;
matthiasm@8:   end
matthiasm@8: end
matthiasm@8: [p(T), path(T)] = max(delta(:,T));
matthiasm@8: for t=T-1:-1:1
matthiasm@8:   path(t) = psi(path(t+1),t+1);
matthiasm@8: end
matthiasm@8: 
matthiasm@8: %  keyboard
matthiasm@8: 
matthiasm@8: % If scaled==0, p = prob_path(best_path)
matthiasm@8: % If scaled==1, p = Pr(replace sum with max and proceed as in the scaled forwards algo)
matthiasm@8: % Both are different from p(data) as computed using the sum-product (forwards) algorithm
matthiasm@8: 
matthiasm@8: if 0
matthiasm@8: if scaled
matthiasm@8:   loglik = -sum(log(scale));
matthiasm@8:   %loglik = prob_path(prior, transmat, obslik, path);
matthiasm@8: else
matthiasm@8:   loglik = log(p);
matthiasm@8: end
matthiasm@8: end