Daniel@0: % Example of fixed lag smoothing
Daniel@0: 
Daniel@0: rand('state', 1);
Daniel@0: S = 2;
Daniel@0: O = 2;
Daniel@0: T = 7;
Daniel@0: data = sample_discrete([0.5 0.5], 1, T);
Daniel@0: transmat = mk_stochastic(rand(S,S));
Daniel@0: obsmat = mk_stochastic(rand(S,O));
Daniel@0: obslik = multinomial_prob(data, obsmat);
Daniel@0: prior = [0.5 0.5]';
Daniel@0: 
Daniel@0: 
Daniel@0: [alpha0, beta0, gamma0, ll0, xi0] = fwdback(prior, transmat, obslik);
Daniel@0: 
Daniel@0: w = 3;
Daniel@0: alpha1 = zeros(S, T);
Daniel@0: gamma1 = zeros(S, T);
Daniel@0: xi1 = zeros(S, S, T-1);
Daniel@0: t = 1;
Daniel@0: b = obsmat(:, data(t));
Daniel@0: olik_win = b; % window of conditional observation likelihoods
Daniel@0: alpha_win = normalise(prior .* b);
Daniel@0: alpha1(:,t) = alpha_win;
Daniel@0: for t=2:T
Daniel@0:   [alpha_win, olik_win, gamma_win, xi_win] = ...
Daniel@0:       fixed_lag_smoother(w, alpha_win, olik_win, obsmat(:, data(t)), transmat);
Daniel@0:   alpha1(:,max(1,t-w+1):t) = alpha_win;
Daniel@0:   gamma1(:,max(1,t-w+1):t) = gamma_win;
Daniel@0:   xi1(:,:,max(1,t-w+1):t-1) = xi_win;
Daniel@0: end
Daniel@0: 
Daniel@0: e = 1e-1;
Daniel@0: %assert(approxeq(alpha0, alpha1, e));
Daniel@0: assert(approxeq(gamma0(:, T-w+1:end), gamma1(:, T-w+1:end), e));
Daniel@0: %assert(approxeq(xi0(:,:,T-w+1:end), xi1(:,:,T-w+1:end), e));
Daniel@0: 
Daniel@0: