Lab2.m

%Kalman Lab 2 demonstrates the kalman filter likelihood profiles for different data sets
%Written by E. E. Holmes (eli.holmes@noaa.gov)

function  Lab2()

%getdata is a generic function that gives data used for the workshop
[data, dataname] = getdata;

%the next few lines set up some reasonable initial values
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);

%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;

disp('Searching the likelihood surface for the maximum');
%kalman_loglik is a function that get the max likelihood fits of CDA to data
%kalman_loglik( a,y,V1, 2); a = parameters mu, s2, s2np, xinit; y = data; V1 = init var initx; 2 means plot
%kalman_loglik returns the negative log likelihood
%fminsearch finds the parameters, a, that minimize the negative loglikelihood
myopt=optimset('Display','off');
initvals = [muinit log(max(0.0001,sigma2init)) log(max(0.0001,sigma2npinit)) y(1)];  %where to start the search
V1 = .001;  %something small seems to work fine; could solve for min, but I don't here
a=fminsearch('kalman_loglik', initvals, myopt, y, V1, 0 );

%now plot it by setting demo flag to 2
h=figure(2);
set(h,'Position',[722   180   668   768]);
subplot(4,1,1)
negloglik = kalman_loglik( a,y,V1, 2);
hold on;
t = 1:length(data);
plot(t(y ~= -99),y(y ~= -99),'.r');
hold off;
h=title(dataname);
set(h,'FontSize',16);
legend('Fit with ML parameters');
ylabel('log data');

%Make the likelihood profiles
disp('Making profiles; this takes awhile');
for(j = 1:3) % do each param
   subplot(4,1,j+1) %make a 3x1 plot
   [profpoints, likeprofile] = kalman_likeprofile(a, y, V1, j); %j=1 is mu; j=2 is s2; j=3 is s2np
   plot(profpoints,likeprofile-min(likeprofile),'-b');
   hold on;
   ylim([0 10]);
   xlim([profpoints(1) profpoints(end)]);
   plot([profpoints(1) profpoints(end)],[1.92/2 1.92/2],'-r');
   hold off;
   if(j==1) title(['mu (ML = ' num2str(a(1)) ')']); end
   if(j==2) title(['process error (ML = ' num2str(exp(a(2))) ')']); end
   if(j==3) title(['non-process error (ML = ' num2str(exp(a(3))) ')']); end
   ylabel('normed -log L');
   legend('L profile','95% estimated level');
end