annotate _beattracker/viterbi_path.m @ 9:4ea6619cb3f5 tip

removed log files
author matthiasm
date Fri, 11 Apr 2014 15:55:11 +0100
parents b5b38998ef3b
children
rev   line source
matthiasm@8 1 function [path,p] = viterbi_path(prior, transmat, obslik)
matthiasm@8 2 % VITERBI Find the most-probable (Viterbi) path through the HMM state trellis.
matthiasm@8 3 % path = viterbi(prior, transmat, obslik)
matthiasm@8 4 %
matthiasm@8 5 % Inputs:
matthiasm@8 6 % prior(i) = Pr(Q(1) = i)
matthiasm@8 7 % transmat(i,j) = Pr(Q(t+1)=j | Q(t)=i)
matthiasm@8 8 % obslik(i,t) = Pr(y(t) | Q(t)=i)
matthiasm@8 9 %
matthiasm@8 10 % Outputs:
matthiasm@8 11 % path(t) = q(t), where q1 ... qT is the argmax of the above expression.
matthiasm@8 12
matthiasm@8 13
matthiasm@8 14 % delta(j,t) = prob. of the best sequence of length t-1 and then going to state j, and O(1:t)
matthiasm@8 15 % psi(j,t) = the best predecessor state, given that we ended up in state j at t
matthiasm@8 16
matthiasm@8 17 scaled = 1;
matthiasm@8 18
matthiasm@8 19 T = size(obslik, 2);
matthiasm@8 20 prior = prior(:);
matthiasm@8 21 Q = length(prior);
matthiasm@8 22
matthiasm@8 23 delta = zeros(Q,T);
matthiasm@8 24 psi = zeros(Q,T);
matthiasm@8 25 path = zeros(1,T);
matthiasm@8 26 scale = ones(1,T);
matthiasm@8 27
matthiasm@8 28
matthiasm@8 29 t=1;
matthiasm@8 30 delta(:,t) = prior .* obslik(:,t);
matthiasm@8 31 if scaled
matthiasm@8 32 [delta(:,t), n] = normalise(delta(:,t));
matthiasm@8 33 scale(t) = 1/n;
matthiasm@8 34 end
matthiasm@8 35 psi(:,t) = 0; % arbitrary value, since there is no predecessor to t=1
matthiasm@8 36 for t=2:T
matthiasm@8 37 for j=1:Q
matthiasm@8 38 [delta(j,t), psi(j,t)] = max(delta(:,t-1) .* transmat(:,j));
matthiasm@8 39 delta(j,t) = delta(j,t) * obslik(j,t);
matthiasm@8 40 end
matthiasm@8 41 if scaled
matthiasm@8 42 [delta(:,t), n] = normalise(delta(:,t));
matthiasm@8 43 scale(t) = 1/n;
matthiasm@8 44 end
matthiasm@8 45 end
matthiasm@8 46 [p(T), path(T)] = max(delta(:,T));
matthiasm@8 47 for t=T-1:-1:1
matthiasm@8 48 path(t) = psi(path(t+1),t+1);
matthiasm@8 49 end
matthiasm@8 50
matthiasm@8 51 % keyboard
matthiasm@8 52
matthiasm@8 53 % If scaled==0, p = prob_path(best_path)
matthiasm@8 54 % If scaled==1, p = Pr(replace sum with max and proceed as in the scaled forwards algo)
matthiasm@8 55 % Both are different from p(data) as computed using the sum-product (forwards) algorithm
matthiasm@8 56
matthiasm@8 57 if 0
matthiasm@8 58 if scaled
matthiasm@8 59 loglik = -sum(log(scale));
matthiasm@8 60 %loglik = prob_path(prior, transmat, obslik, path);
matthiasm@8 61 else
matthiasm@8 62 loglik = log(p);
matthiasm@8 63 end
matthiasm@8 64 end