wolffd@0: function [xsmooth, Vsmooth, VVsmooth, loglik] = kalman_smoother(y, A, C, Q, R, init_x, init_V, varargin)
wolffd@0: % Kalman/RTS smoother.
wolffd@0: % [xsmooth, Vsmooth, VVsmooth, loglik] = kalman_smoother(y, A, C, Q, R, init_x, init_V, ...)
wolffd@0: %
wolffd@0: % The inputs are the same as for kalman_filter.
wolffd@0: % The outputs are almost the same, except we condition on y(:, 1:T) (and u(:, 1:T) if specified),
wolffd@0: % instead of on y(:, 1:t).
wolffd@0: 
wolffd@0: [os T] = size(y);
wolffd@0: ss = length(A);
wolffd@0: 
wolffd@0: % set default params
wolffd@0: model = ones(1,T);
wolffd@0: u = [];
wolffd@0: B = [];
wolffd@0: 
wolffd@0: args = varargin;
wolffd@0: nargs = length(args);
wolffd@0: for i=1:2:nargs
wolffd@0:   switch args{i}
wolffd@0:    case 'model', model = args{i+1};
wolffd@0:    case 'u', u = args{i+1};
wolffd@0:    case 'B', B = args{i+1};
wolffd@0:    otherwise, error(['unrecognized argument ' args{i}])
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: xsmooth = zeros(ss, T);
wolffd@0: Vsmooth = zeros(ss, ss, T);
wolffd@0: VVsmooth = zeros(ss, ss, T);
wolffd@0: 
wolffd@0: % Forward pass
wolffd@0: [xfilt, Vfilt, VVfilt, loglik] = kalman_filter(y, A, C, Q, R, init_x, init_V, ...
wolffd@0: 					       'model', model, 'u', u, 'B', B);
wolffd@0: 
wolffd@0: % Backward pass
wolffd@0: xsmooth(:,T) = xfilt(:,T);
wolffd@0: Vsmooth(:,:,T) = Vfilt(:,:,T);
wolffd@0: %VVsmooth(:,:,T) = VVfilt(:,:,T);
wolffd@0: 
wolffd@0: for t=T-1:-1:1
wolffd@0:   m = model(t+1);
wolffd@0:   if isempty(B)
wolffd@0:     [xsmooth(:,t), Vsmooth(:,:,t), VVsmooth(:,:,t+1)] = ...
wolffd@0: 	smooth_update(xsmooth(:,t+1), Vsmooth(:,:,t+1), xfilt(:,t), Vfilt(:,:,t), ...
wolffd@0: 		      Vfilt(:,:,t+1), VVfilt(:,:,t+1), A(:,:,m), Q(:,:,m), [], []);
wolffd@0:   else
wolffd@0:     [xsmooth(:,t), Vsmooth(:,:,t), VVsmooth(:,:,t+1)] = ...
wolffd@0: 	smooth_update(xsmooth(:,t+1), Vsmooth(:,:,t+1), xfilt(:,t), Vfilt(:,:,t), ...
wolffd@0: 		      Vfilt(:,:,t+1), VVfilt(:,:,t+1), A(:,:,m), Q(:,:,m), B(:,:,m), u(:,t+1));
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: VVsmooth(:,:,1) = zeros(ss,ss);
wolffd@0: