Lab1.m

%Kalman Lab 1 demonstrates the kalman filter step by step for a simple CDA
%User inputs parameter values for the CDA
%Written by E. E. Holmes (eli.holmes@noaa.gov)

function kalman_lab1()

%get some data
[data, dataname] = getdata;

h=figure(1); %this will be a 4x1 plot
clf %clear figure
set(h,'Position',[676  131   668   768]);

%next for lines are to properly deal with missing data so that kalman filter knows they are missing
tmp = data;
tmp(data==-99) = 1;
y = log(tmp);
y(data==-99) = -99;

t = 1:length(data);

for (j = 1:3) %make plots with user input params 3x
   subplot(4,1,j)
   
   %get some user inputs
   disp(' ');
   disp(['Panel ' num2str(j)]);
   mu = input('input mu: ');
   s2 = input('input process error (not 0): ');
   s2np = input('input non-process error (not 0): ');
   initx = y(1); %for the purpose of this lab, init x set at y(1), but normally one minimizes for this
   V1 = .01;
   
   %These are the params set by the user; this defines a CDA
   userparams = [mu log(s2) log(s2np) y(1)];
   
   %now the kalman filter is used to fit this CDA to the data; gives the ML estimate of true values
   %kalman_loglik gets passed 
   %userparams = CDA params; y = data; V1 = init var initx; 2 = yes plot
   negloglik = kalman_loglik(userparams, y, V1,2 );
   
   title(['CDA params = ' num2str([mu s2 s2np])  '; -log L = ' num2str(negloglik)]);
   legend('ML fit');
   if(j==2) ylabel('log data'); end
end

%for the last plot show the fit with the parameters that minimize the fit
subplot(4,1,4)

%set up reasonable starting values for the search
Nratios = data(2:end)./data(1:(end-1));
Nratios = Nratios( Nratios > 0 );
Nratios = log(Nratios);

muinit = mean(Nratios);
sigma2init = 0.5 * var(Nratios);
sigma2npinit = 0.5 * var(Nratios);
myopt=optimset('Display','off');

%fminsearch finds the parameters, a, that minimize the negative loglikelihood
initvals = [muinit log(max(0.0001,sigma2init)) log(max(0.0001,sigma2npinit)) y(1)];  %where to start the search
a=fminsearch('kalman_loglik', initvals ,myopt,y, V1, 0 );

%get the negloglik for the ML parameters; plot it
negloglik = kalman_loglik( a,y,V1, 2);

title(['MLE CDA params = ' num2str([a(1) exp(a(2)) exp(a(3))])  '; -log L = ' num2str(negloglik)]);
legend('Fit with ML parameters');

h=xlabel(dataname);
set(h,'FontSize',16);