function [alpha, gamma, xi] = fixed_lag_smoother(D, alpha, obsmat, transmat, act)% FIXED_LAG_SMOOTHER% [alpha, gamma, xi] = fixed_lag_smoother(D, alpha, obsmat, transmat, act) %% D is the delay/lag, so D=0 means online filtering.% d=min(D+1, t0) is the size of the window we need to keep as history,% where t0 is the current time.%% alpha(:, t0-d-1:t0-1) - length d window, excluding t0 (Columns indexed 1..d)% obsmat(:, t0-d:t0)    - length d window, last column = likelihood vector for current observation% transmat            - transition matrix% If we specify the optional 'act' argument, transmat{a} should be a cell array, and% act(t0-d:t0)          - length d window, last column = current action%% Output:% alpha(:, t0-d:t0)     - last column = new filtered estimate% xi(:, :, t0-d:t0-1)   - 2 slice smoothed window% gamma(:, t0-d:t0)     - smoothed window%% As usual, we define (using T=d)% alpha(i,t) = Pr(Q(t)=i | Y(1:t))% gamma(i,t) = Pr(Q(t)=i | Y(1:T))% xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | Y(1:T))%% obsmat(i,t) = Pr(Y(t) | Q(t)=i) - use mk_dhmm_obs_mat to create this% transmat{a}(i,j) = Pr(Q(t)=j | Q(t-1)=i, A(t)=a)[S n] = size(alpha); if nargin < 5  act = ones(1, n+1);  transmat = { transmat };endd = D+1;d = min(d, n+1); % we can only grow the window to 1 + its current sized = max(d, 2); % we must always use at least 2 slicesalpha = alpha(:, n-d+2:n); % pluck out last d-1 components% Extend window by 1t = d;xi = zeros(S, S, d-1);xi(:,:,t-1) = normalise((alpha(:,t-1) * obsmat(:,t)') .* transmat{act(t)});alpha(:,t) = sum(xi(:,:,t-1), 1)';% Now smooth backwards inside the windowif D > 0  beta = ones(S, d);  T = d;  gamma(:,T) = alpha(:,T);  for t=T-1:-1:1    b = beta(:,t+1) .* obsmat(:,t+1);     beta(:,t) = normalise(transmat{act(t)} * b);    gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));    xi(:,:,t) = normalise((transmat{act(t)} .* (alpha(:,t) * b')));  endend