How do I use fminsearch to optimize the likelihood function for a Kalman Filter?
3 ビュー (過去 30 日間)
古いコメントを表示
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
0 件のコメント
回答 (2 件)
Alan Weiss
2013 年 5 月 7 日
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
Alan Weiss
2013 年 5 月 7 日
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
2 件のコメント
Alan Weiss
2013 年 5 月 7 日
- 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
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!