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