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