function [x, V, VV, loglik] = kalman_filter(y, A, C, Q, R, init_x, init_V, model)

% Kalman filter.

% [x, V, VV, loglik] = kalman_filter(y, A, C, Q, R, init_x, init_V, model)

%

% Inputs:

% y(:,t)   - the observation at time t

% A(:,:,m) - the system matrix for model m

% C(:,:,m) - the observation matrix for model m

% Q(:,:,m) - the system covariance for model m

% R(:,:,m) - the observation covariance for model m

% init_x(:,m) - the initial state for model m

% init_V(:,:,m) - the initial covariance for model m

% model(t) - which model to use at time t (defaults to model 1 if not specified)

%

% Outputs:

% x(:,t) = E[X_t | t]

% V(:,:,t) = Cov[X_t | t]

% VV(:,:,t) = Cov[X_t, X_t-1 | t] t >= 2

% loglik = sum_t log P(Y_t)



[os T] = size(y);

ss = size(A,1);



if nargin<8, model = ones(1, T); end



x = zeros(ss, T);

V = zeros(ss, ss, T);

VV = zeros(ss, ss, T);



loglik = 0;

for t=1:T

  m = model(t);

  if t==1

    prevx = init_x(:,m);

    prevV = init_V(:,:,m);

    initial = 1;

  else

    prevx = x(:,t-1);

    prevV = V(:,:,t-1);

    initial = 0;

  end

  [x(:,t), V(:,:,t), LL, VV(:,:,t)] = ...

      kalman_update(A(:,:,m), C(:,:,m), Q(:,:,m), R(:,:,m), y(:,t), prevx, prevV, initial);

  loglik = loglik + LL;

end