% LynxHare.m
%
% This m-file does parameter estimation for the Hudson Bay Company's
% Lynx-Hare data set using the predator-prey model
%
%  dX1/dt =  b(1)*X1 - b(2)*X1*X2
%  dX2/dt = -b(4)*X2 + b(3)*X1*X2
%  X1(0) = X1_0, X2(0) = X2_0 
clear all, close all
%--------------------------------------------------------------------------
% Input data and other problem dependent functions and parameters.
%--------------------------------------------------------------------------
H = [30 47.2 70.2 77.4 36.3 20.6 18.1 21.4 22 25.4 27.1 ...
        40.3 57 76.6 52.3 19.5 11.2 7.6 14.6 16.2 24.7]'; %Hare
L = [4 6.1 9.8 35.2 59.4 41.7 19 13 8.3 9.1 7.4 ...
        8 12.3 19.5 45.7 51.1 29.7 15.8 9.7 10.1 8.6]'; %Linx
yvec=[H;L]; 
tvec=zeros(size(yvec)); tvec(1:21)=1:21;
nlin_fn = @lynxhare_fn;
b_0 = [.47; .024; .023; .76; 30; 4]; 
nparams = length(b_0);
npts = length(yvec);

%--------------------------------------------------------------------------
% Compute optimal parameter values using nlinfit and plot.
%--------------------------------------------------------------------------
[b_opt,r,J,C,mse] = nlinfit(tvec,yvec,nlin_fn,b_0);
model_opt = feval(nlin_fn,b_opt,tvec);
figure(1)
  subplot(211)
    plot(tvec(1:npts/2),yvec(1:npts/2),'o',tvec(1:npts/2),model_opt(1:npts/2),'k')
    title('Data and Model Fit: Hare')
  subplot(212)
    plot(tvec(1:npts/2),yvec(npts/2+1:npts),'o'), 
    hold on, plot(tvec(1:21),model_opt(22:42),'k'), hold off
    title('Data and Model Fit: Linx')

%--------------------------------------------------------------------------
% Sample from beta_hat using Monte Carlo Method.
%--------------------------------------------------------------------------
ndraws = 1000; 
b_est = zeros(nparams,ndraws);
for i=1:ndraws
  h = waitbar(i/ndraws);
  e = sqrt(mse)*randn(npts,1);
  y = model_opt+e;
  b_est(:,i) = nlinfit(tvec,y,nlin_fn,b_0);
end
close(h)
    
%--------------------------------------------------------------------------
% Plot samples and output confidence intervals
%--------------------------------------------------------------------------
sample_plot(b_est,b_opt,2)
ci = nlparci(b_opt,r,'jacobian',J);
for i=1:nparams
  fprintf('b_%d normal theory: 0.025 quantile = %2.5f, mean = %2.5f 0.975 quantile = %2.5f\n',i,ci(i,1),mean([ci(i,1),ci(i,2)]),ci(i,2))
end