www.pudn.com > HMM1.zip > fixed_lag_smoother.m
function [alpha, obslik, gamma, xi] = fixed_lag_smoother(d, alpha, obslik, obsvec, transmat, act)
> FIXED_LAG_SMOOTHER Computed smoothed posterior estimates within a window given previous filtered window.
> [alpha, obslik, gamma, xi] = fixed_lag_smoother(d, alpha, obslik, obsvec, transmat, act)
>
> d >= 2 is the desired window width.
> Actually, we use d=min(d, t0), where t0 is the current time.
>
> alpha(:, t0-d:t0-1) - length d window, excluding t0 (Columns indexed 1..d)
> obslik(:, t0-d:t0-1) - length d window
> obsvec - 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+1:t0) - last column = new filtered estimate
> obslik(:, t0-d+1:t0) - last column = obsvec
> xi(:, :, t0-d+1:t0-1) - 2 slice smoothed window
> gamma(:, t0-d+1: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))
>
> obslik(i,t) = Pr(Y(t) | Q(t)=i)
> transmat{a}(i,j) = Pr(Q(t)=j | Q(t-1)=i, A(t)=a)
[S n] = size(alpha);
d = min(d, n+1);
if d < 2
error('must keep a window of length at least 2');
end
if ~exist('act')
act = ones(1, n+1);
transmat = { transmat };
end
> pluck out last d-1 components from the history
alpha = alpha(:, n-d+2:n);
obslik = obslik(:, n-d+2:n);
> Extend window by 1
t = d;
obslik(:,t) = obsvec;
xi = zeros(S, S, d-1);
xi(:,:,t-1) = normalise((alpha(:,t-1) * obslik(:,t)') .* transmat{act(t)});
alpha(:,t) = sum(xi(:,:,t-1), 1)';
> Now smooth backwards inside the window
beta = ones(S, d);
T = d;
>fprintf('smooth from >d to 1, i.e., >d to >d\n', d, t0, t0-d+1);
gamma(:,T) = alpha(:,T);
for t=T-1:-1:1
b = beta(:,t+1) .* obslik(:,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')));
end