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
|