wolffd@0: function [alpha, obslik, gamma, xi] = fixed_lag_smoother(d, alpha, obslik, obsvec, transmat, act)
wolffd@0: % FIXED_LAG_SMOOTHER Computed smoothed posterior estimates within a window given previous filtered window.
wolffd@0: % [alpha, obslik, gamma, xi] = fixed_lag_smoother(d, alpha, obslik, obsvec, transmat, act)
wolffd@0: %
wolffd@0: % d >= 2 is the desired window width.
wolffd@0: % Actually, we use d=min(d, t0), where t0 is the current time.
wolffd@0: %
wolffd@0: % alpha(:, t0-d:t0-1) - length d window, excluding t0 (Columns indexed 1..d)
wolffd@0: % obslik(:, t0-d:t0-1) - length d window
wolffd@0: % obsvec              - likelihood vector for current observation
wolffd@0: % transmat            - transition matrix
wolffd@0: % If we specify the optional 'act' argument, transmat{a} should be a cell array, and
wolffd@0: % act(t0-d:t0)          - length d window, last column = current action
wolffd@0: %
wolffd@0: % Output:
wolffd@0: % alpha(:, t0-d+1:t0)     - last column = new filtered estimate
wolffd@0: % obslik(:, t0-d+1:t0)    - last column = obsvec
wolffd@0: % xi(:, :, t0-d+1:t0-1)   - 2 slice smoothed window
wolffd@0: % gamma(:, t0-d+1:t0)     - smoothed window
wolffd@0: %
wolffd@0: % As usual, we define (using T=d)
wolffd@0: % alpha(i,t) = Pr(Q(t)=i | Y(1:t))
wolffd@0: % gamma(i,t) = Pr(Q(t)=i | Y(1:T))
wolffd@0: % xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | Y(1:T))
wolffd@0: %
wolffd@0: % obslik(i,t) = Pr(Y(t) | Q(t)=i) 
wolffd@0: % transmat{a}(i,j) = Pr(Q(t)=j | Q(t-1)=i, A(t)=a)
wolffd@0: 
wolffd@0: [S n] = size(alpha); 
wolffd@0: d = min(d, n+1); 
wolffd@0: if d < 2
wolffd@0:   error('must keep a window of length at least 2');
wolffd@0: end
wolffd@0: 
wolffd@0: if ~exist('act')
wolffd@0:   act = ones(1, n+1);
wolffd@0:   transmat = { transmat };
wolffd@0: end
wolffd@0: 
wolffd@0: % pluck out last d-1 components from the history
wolffd@0: alpha = alpha(:, n-d+2:n); 
wolffd@0: obslik = obslik(:, n-d+2:n); 
wolffd@0: 
wolffd@0: % Extend window by 1
wolffd@0: t = d;
wolffd@0: obslik(:,t) = obsvec;
wolffd@0: xi = zeros(S, S, d-1);
wolffd@0: xi(:,:,t-1) = normalise((alpha(:,t-1) * obslik(:,t)') .* transmat{act(t)});
wolffd@0: alpha(:,t) = sum(xi(:,:,t-1), 1)';
wolffd@0: 
wolffd@0: % Now smooth backwards inside the window
wolffd@0: beta = ones(S, d);
wolffd@0: T = d;
wolffd@0: %fprintf('smooth from %d to 1, i.e., %d to %d\n', d, t0, t0-d+1);
wolffd@0: gamma(:,T) = alpha(:,T);
wolffd@0: for t=T-1:-1:1
wolffd@0:   b = beta(:,t+1) .* obslik(:,t+1); 
wolffd@0:   beta(:,t) = normalise(transmat{act(t)} * b);
wolffd@0:   gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
wolffd@0:   xi(:,:,t) = normalise((transmat{act(t)} .* (alpha(:,t) * b')));
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: