comparison toolboxes/FullBNT-1.0.7/HMM/fixed_lag_smoother.m @ 0:e9a9cd732c1e tip

first hg version after svn

+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
+
+
+
+