get_kalman_ests.m

%This version of Kalman is latest as of July 2005
%It corrects for missing values
%It estimates the initx
%Written by E. E. Holmes (eli.holmes@noaa.gov)

function kalman=get_kalman_ests(data)
%data is a row vector of raw counts
%here I use Dennis et al. method to get some reasonable init values for mu and sigma2
dennis_ests = get_dennis_ests(data);

muinit = dennis_ests.mu;
sigma2init = 0.5 * dennis_ests.s2p;
sigma2npinit = 0.5 * dennis_ests.s2p;
myopt=optimset('Display','off');

%deal with 0s by adding +1
if(ismember(0,data)) data=data+1; end
data(data==-98) = -99;  %need to adjust for having added +1 to the -99s

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

initvals = [muinit log(max(0.0001,sigma2init)) log(max(0.0001,sigma2npinit)) y(1)];
V1 = .001; %generic small seems to work ok
a=fminsearch('kalman_loglik_initx_free', initvals ,myopt,y,V1,0); %0 is a flag to say don't plot

kalman.mu=a(1);
kalman.s2p=exp(a(2));
kalman.s2s=exp(a(3));
kalman.initx=a(4);

%x is true N estimate
[kalman.negloglik kalman.x]= kalman_loglik_initx_free( a,y, V1, 0 ); %0 is just a flag to say don't plot