How do I use fminsearch to optimize the likelihood function for a Kalman Filter?

Hi,

I have written a Kalman Filter which works and I would now like to find the parameters which optimize the likelihood, using fminsearch. Can anyone help me with this?

This is my Kalman Filter Code and below is the function I have so far to maximize the likelihood.

    function [loglik,errory,FC,FS,mse_store]= kf_loglik_adapted(data, C, xt_init, R, A, Q, sigma)
    %
    N=max(size(data,1));
    xt_upd=xt_init;
    xt_pred =[];

    sigma_upd=sigma;
    sigma_pred=[];
    mse_store=[];
    mse_store=sigma;

    pred_error=[];
    mse=[];

    %
    FC=zeros(N,1);
    FS=zeros(N,6);
    errory = [];
    loglik = 0;
    %Run over the sample

    for i=1:N

         %Prediction
         xt_pred =  A.*xt_upd;  % (2x6)=(2x6).*(2x6)
         sigma_pred = A'*A*sigma_upd + Q; %(6x6)=(6x2)*(2x6)*(6x6)+(6x6)

         %Updating
         pred_error = data(i,:)-diag(C*xt_pred)'; %(1x6) = (1x6)- (1x6)

         mse = C*C' * sigma_pred + R; %(6x6)=(6x2)*(2x6)*(6x6)+(6x6)
         inv_mse = inv(mse); %(6x6)

         %kalg  2x6
         kalg = C'*sigma_pred *inv_mse; %(2x6)=(2x6)*(6x6)*(6x6)

         xt_upd = xt_pred + kalg .* repmat(pred_error,2,1); %(2x6)=(2x6)+ (2x6).*(2x6)
         sigma_upd = (eye(6)-C*kalg)*sigma_pred; %(6x6)=((6x6)-(6x2)*(2x6))*(6x6)

         % store factors
            factors=xt_upd./A;
            FC(i,1) = factors(1,1);
            FS(i,:) = factors(2,:);

         % Store Errors
            errory(i,:) = data(i,:)-diag(C*xt_upd)';

         % store sigma
            if i==1
                mse_store(1:6,1:6) = mse;
            else
                mse_store(6*i+1:6*i+6,1:6) = mse;
            end
    end

    %Log-Likelihood function 
    loglik=loglik-0.5*log(abs(det(mse)))-0.5* errory(end,:)*inv(mse)*errory(end,:)';

function result = maxlik(sigma_store,obs,errory)

for i=1:obs

maxlik = -0.5*log(abs(sigma_store(i*6+1:i*6,1:6)))...
    -0.5* errory(i,:)*inv(sigma_store(i*6+1:i*6,1:6))*errory(i,:)'; 

end

Ideally I need parameters A, C. When I use the optimtool it always tells me I don't have enough input parameters. the help doesn't seem to get me much further

thanks

Answer 1

[Alan Weiss]

As explained in the documentation, to write an objective function, you need to put all your variables into one vector, typically called x. So your objective function should look something like the following:

 function [loglik,errory,FC,FS,mse_store]= kf_loglik_adapted(x)
 A = x(1)
 B = x(2)
 C = x(3)
 % Calculations go here

As explained in the documentation, if you need to pass extra parameters such as data and xt_init, write an anonymous function:

 function [loglik,errory,FC,FS,mse_store]= kf_loglik_adapted(x,data,xt_init)
 ...

Then call the function with the objective

 obj = @(x)kf_loglik_adapted(x,data,xt_init)
 xanswer = fminsearch(obj,x0)

Alan Weiss

MATLAB mathematical toolbox documentation

[Steven]

Hi Alan, thanks for your answer. I've actually changed the code to include the sum of the log-likelihoods and it works.

I need the first row of my vector params to take on the same value. is there any way I can do that with options?

I now use: params(1:2,:) = A ; params(3:4,:) = C';

results = fminsearch('kalmanf',params,[],data,xt_init,sigma,Q,R); 

function [sumll,loglik,errory,FC,FS,mse_store]= kalmanf(params, data, xt_init, sigma,Q,R)

% get params
A = params(1:2,:);
C = params(3:4,:)';

% other initializations
N=max(size(data,1));
xt_upd=xt_init;
xt_pred =[];

sigma_upd=sigma;
sigma_pred=[];
mse_store=[];
mse_store=sigma;

pred_error=[];
mse=[];

%
FC=zeros(N,1);
FS=zeros(N,6);
errory = [];
loglik = [];
%Run over the sample

for i=1:N

       %Prediction
       xt_pred =  A.*xt_upd;  % (2x6)=(2x6).*(2x6)
       sigma_pred = A'*A*sigma_upd + Q; %(6x6)=(6x2)*(2x6)*(6x6)+(6x6)

       %Updating
       pred_error = data(i,:)-diag(C*xt_pred)'; %(1x6) = (1x6)- (1x6)

       mse = C*C' * sigma_pred + R; %(6x6)=(6x2)*(2x6)*(6x6)+(6x6)
       inv_mse = inv(mse); %(6x6)

       %kalg  2x6
       kalg = C'*sigma_pred *inv_mse; %(2x6)=(2x6)*(6x6)*(6x6)

       xt_upd = xt_pred + kalg .* repmat(pred_error,2,1); %(2x6)=(2x6)+ (2x6).*(2x6)
       sigma_upd = (eye(6)-C*kalg)*sigma_pred; %(6x6)=((6x6)-(6x2)*(2x6))*(6x6)

       % store factors
          factors = xt_upd./A; % does this have to be divided by A -- xt_pred is defined as A*xt_upd
          FC(i,1) = factors(1,1);
          FS(i,:) = factors(2,:);

       % Store Errors
          errory(i,:) = data(i,:)-diag(C*xt_upd)';

       % store sigma
          if i==1
              mse_store(1:6,1:6) = mse;
          else
              mse_store(6*i+1:6*i+6,1:6) = mse;
          end

       % store ll
          loglik(i,1)=0-0.5*log(abs(det(mse)))-0.5* errory(i,:)*inv(mse)*errory(i,:)'; %constant set to zero for now

end

%Log-Likelihood function 
sumll = -sum(loglik);

thanks

[Steven]

I suppose I will need to use fmincon, but I don't know how to specify the constraint for this

Answer 2

[Alan Weiss]

You need the first two rows of your vector parameter to be the same? Then why not define

 A = param(1)
 B = param(1)

or some such thing. I mean, just reduce the number of variables that fminsearch uses.

Also, if you have Optimization Toolbox, you will almost certainly get better performance using fminunc than fminsearch.

If, for some reason, you really don't want to rewrite your code to reduce the number of parameters, you can use fmincon with a linear equality constraint of the form

 Aeq*x = beq

where Aeq is a matrix with rows of the form [+1,-1,0,0] and beq = 0. This means

 x(1) - x(2) = 0

which is what you want, I think. But I also think it would be more efficient to reduce the number of variables rather than use a linear equality constraint.

Alan Weiss

MATLAB mathematical toolbox documentation

[Steven]

No sorry I wasn't very clear I think, the goal is that the first row of my params matrix (a 4x6 matrix) takes on the same values. So I need something like this: params(1,1) = params(1,2); params(1,2) = params(1,3); params(1,3) = params(1,4); params(1,4) = params(1,5); params(1,5) = params(1,6);

Can I still use the Aeq*x = beq you proposed? How would I defined the x in this case to only be the first row of the params matrix?

thanks

[Alan Weiss]

• I still think that the best way to approach this is to change your problem to have fewer parameters. You have only four rows, so four parameters. If x is your vector of four unknowns, you could set

 params = [x(:),x(:),x(:),x(:),x(:),x(:)];

• However, if you really don't want to do things that way, then of course you can use linear equalities. You must convert the params matrix to a vector params(:), and use Aeq and beq entries that operate on that form. There are 24 entries in params, so Aeq has 24 columns. Each row of Aeq has a +1 and a -1, with zeros as the other entries.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation