|
wolffd@0
|
1 function [alpha, obslik, gamma, xi] = fixed_lag_smoother(d, alpha, obslik, obsvec, transmat, act)
|
|
wolffd@0
|
2 % FIXED_LAG_SMOOTHER Computed smoothed posterior estimates within a window given previous filtered window.
|
|
wolffd@0
|
3 % [alpha, obslik, gamma, xi] = fixed_lag_smoother(d, alpha, obslik, obsvec, transmat, act)
|
|
wolffd@0
|
4 %
|
|
wolffd@0
|
5 % d >= 2 is the desired window width.
|
|
wolffd@0
|
6 % Actually, we use d=min(d, t0), where t0 is the current time.
|
|
wolffd@0
|
7 %
|
|
wolffd@0
|
8 % alpha(:, t0-d:t0-1) - length d window, excluding t0 (Columns indexed 1..d)
|
|
wolffd@0
|
9 % obslik(:, t0-d:t0-1) - length d window
|
|
wolffd@0
|
10 % obsvec - likelihood vector for current observation
|
|
wolffd@0
|
11 % transmat - transition matrix
|
|
wolffd@0
|
12 % If we specify the optional 'act' argument, transmat{a} should be a cell array, and
|
|
wolffd@0
|
13 % act(t0-d:t0) - length d window, last column = current action
|
|
wolffd@0
|
14 %
|
|
wolffd@0
|
15 % Output:
|
|
wolffd@0
|
16 % alpha(:, t0-d+1:t0) - last column = new filtered estimate
|
|
wolffd@0
|
17 % obslik(:, t0-d+1:t0) - last column = obsvec
|
|
wolffd@0
|
18 % xi(:, :, t0-d+1:t0-1) - 2 slice smoothed window
|
|
wolffd@0
|
19 % gamma(:, t0-d+1:t0) - smoothed window
|
|
wolffd@0
|
20 %
|
|
wolffd@0
|
21 % As usual, we define (using T=d)
|
|
wolffd@0
|
22 % alpha(i,t) = Pr(Q(t)=i | Y(1:t))
|
|
wolffd@0
|
23 % gamma(i,t) = Pr(Q(t)=i | Y(1:T))
|
|
wolffd@0
|
24 % xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | Y(1:T))
|
|
wolffd@0
|
25 %
|
|
wolffd@0
|
26 % obslik(i,t) = Pr(Y(t) | Q(t)=i)
|
|
wolffd@0
|
27 % transmat{a}(i,j) = Pr(Q(t)=j | Q(t-1)=i, A(t)=a)
|
|
wolffd@0
|
28
|
|
wolffd@0
|
29 [S n] = size(alpha);
|
|
wolffd@0
|
30 d = min(d, n+1);
|
|
wolffd@0
|
31 if d < 2
|
|
wolffd@0
|
32 error('must keep a window of length at least 2');
|
|
wolffd@0
|
33 end
|
|
wolffd@0
|
34
|
|
wolffd@0
|
35 if ~exist('act')
|
|
wolffd@0
|
36 act = ones(1, n+1);
|
|
wolffd@0
|
37 transmat = { transmat };
|
|
wolffd@0
|
38 end
|
|
wolffd@0
|
39
|
|
wolffd@0
|
40 % pluck out last d-1 components from the history
|
|
wolffd@0
|
41 alpha = alpha(:, n-d+2:n);
|
|
wolffd@0
|
42 obslik = obslik(:, n-d+2:n);
|
|
wolffd@0
|
43
|
|
wolffd@0
|
44 % Extend window by 1
|
|
wolffd@0
|
45 t = d;
|
|
wolffd@0
|
46 obslik(:,t) = obsvec;
|
|
wolffd@0
|
47 xi = zeros(S, S, d-1);
|
|
wolffd@0
|
48 xi(:,:,t-1) = normalise((alpha(:,t-1) * obslik(:,t)') .* transmat{act(t)});
|
|
wolffd@0
|
49 alpha(:,t) = sum(xi(:,:,t-1), 1)';
|
|
wolffd@0
|
50
|
|
wolffd@0
|
51 % Now smooth backwards inside the window
|
|
wolffd@0
|
52 beta = ones(S, d);
|
|
wolffd@0
|
53 T = d;
|
|
wolffd@0
|
54 %fprintf('smooth from %d to 1, i.e., %d to %d\n', d, t0, t0-d+1);
|
|
wolffd@0
|
55 gamma(:,T) = alpha(:,T);
|
|
wolffd@0
|
56 for t=T-1:-1:1
|
|
wolffd@0
|
57 b = beta(:,t+1) .* obslik(:,t+1);
|
|
wolffd@0
|
58 beta(:,t) = normalise(transmat{act(t)} * b);
|
|
wolffd@0
|
59 gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
|
|
wolffd@0
|
60 xi(:,:,t) = normalise((transmat{act(t)} .* (alpha(:,t) * b')));
|
|
wolffd@0
|
61 end
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63
|
|
wolffd@0
|
64
|
|
wolffd@0
|
65
|
