wolffd@0
|
1 function [gamma, loglik, marginals, marginalsT] = bk_ff_fb(prior, transmat, obslik, filter_only, hnodes, ns)
|
wolffd@0
|
2 % BK_FF_FB Fully factored Boyen-Koller version of forwards-backwards
|
wolffd@0
|
3 % [gamma, loglik, marginals, marginalsT] = bk_ff_hmm(prior, transmat, obslik, filter_only, hnodes, ns)
|
wolffd@0
|
4
|
wolffd@0
|
5 ss = length(ns);
|
wolffd@0
|
6 S = length(prior);
|
wolffd@0
|
7 T = size(obslik, 2);
|
wolffd@0
|
8 marginals = cell(ss,T);
|
wolffd@0
|
9 marginalsT = cell(ss,T);
|
wolffd@0
|
10 scale = zeros(1,T);
|
wolffd@0
|
11 alpha = zeros(S, T);
|
wolffd@0
|
12
|
wolffd@0
|
13 transmat2 = transmat';
|
wolffd@0
|
14 for t=1:T
|
wolffd@0
|
15 if t==1
|
wolffd@0
|
16 [alpha(:,t), scale(t)] = normalise(prior(:) .* obslik(:,t));
|
wolffd@0
|
17 else
|
wolffd@0
|
18 [alpha(:,t), scale(t)] = normalise((transmat2 * alpha(:,t-1)) .* obslik(:,t));
|
wolffd@0
|
19 end
|
wolffd@0
|
20 [marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(alpha(:,t), hnodes, ns);
|
wolffd@0
|
21 alpha(:,t) = combine_marginals_into_joint(marginalsT(:,t), hnodes, ns);
|
wolffd@0
|
22 %fprintf('alpha t=%d\n', t);
|
wolffd@0
|
23 %celldisp(marginals(1:8,t))
|
wolffd@0
|
24 end
|
wolffd@0
|
25 loglik = sum(log(scale));
|
wolffd@0
|
26
|
wolffd@0
|
27 if filter_only
|
wolffd@0
|
28 gamma = alpha;
|
wolffd@0
|
29 return;
|
wolffd@0
|
30 end
|
wolffd@0
|
31
|
wolffd@0
|
32 beta = zeros(S,T);
|
wolffd@0
|
33 gamma = zeros(S,T);
|
wolffd@0
|
34 t = T;
|
wolffd@0
|
35 beta(:,t) = ones(S,1);
|
wolffd@0
|
36 gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
|
wolffd@0
|
37 [marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(gamma(:,t), hnodes, ns);
|
wolffd@0
|
38
|
wolffd@0
|
39 for t=T-1:-1:1
|
wolffd@0
|
40 b = beta(:,t+1) .* obslik(:,t+1);
|
wolffd@0
|
41 beta(:,t) = normalise((transmat * b));
|
wolffd@0
|
42 [junk, tempT] = project_joint_onto_marginals(beta(:,t), hnodes, ns);
|
wolffd@0
|
43 beta(:,t) = combine_marginals_into_joint(tempT, hnodes, ns);
|
wolffd@0
|
44 %gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
|
wolffd@0
|
45 %[marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(gamma(:,t), hnodes, ns);
|
wolffd@0
|
46 end
|
wolffd@0
|
47
|
wolffd@0
|
48 gamma2 = zeros(S,T);
|
wolffd@0
|
49 for t=T-1:-1:1
|
wolffd@0
|
50 b = beta(:,t+1) .* obslik(:,t+1);
|
wolffd@0
|
51 xi(:,:,t) = normalise((transmat .* (alpha(:,t) * b')));
|
wolffd@0
|
52 if t==T-1
|
wolffd@0
|
53 gamma2(:,T) = sum(xi(:,:,T-1), 1)';
|
wolffd@0
|
54 end
|
wolffd@0
|
55 gamma2(:,t) = sum(xi(:,:,t), 2);
|
wolffd@0
|
56 [marginals(:,t), marginalsT(:,t)] = project_joint_onto_marginals(gamma2(:,t), hnodes, ns);
|
wolffd@0
|
57 end
|
wolffd@0
|
58
|
wolffd@0
|
59
|