% load data ..
load chl_spring08

% y contains chlorophyll fluorescence [micro g/l]
% x contains time in days since 2001-1-1 [UTC]

% QC: clean up NaNs in y
bad = find(isnan(y));
y(bad) = [];
x(bad) = [];

% modify x so that it's relative to beginning of the data set
x = x-x(1);

% plot (need to get an idea of values for initial parameter guess)
figure
subplot(1,2,1)
plot(x,y,'g*-')
subplot(1,2,2)
semilogy(x,y,'g*-')

% first guess: baseline 0.2, ini pop 0.2, specific growth rate 0.2 
% (doubling in five days) 
ain = [0.2 0.2 0.2];

% plot initial guess model on top of your data
figure
subplot(1,2,1)
plot(x,y,'g*-')
hold on
plot(x,modfunc1(x,ain),'k')
%
subplot(1,2,2)
semilogy(x,y,'g*-')
hold on
semilogy(x,modfunc1(x,ain),'k')

% NOW FIT to model function P = P1 + P2*exp(mu*t) 
% or a1 + a2*exp(a3*t)
% Note: First we must create a matlab function for this model
% in our case called modfunc1.m
% then use nlleasqr
% look at >> help nlleasqr  to learn about input/output parameters
% but don't get intimidated
% for now, call as below
[f,a,kvg,iter,corp,covp,covr,stdresid,Z,r2] = nlleasqr(x,y,ain,'modfunc1');
% we'll only look at a few of the output parameters:
% a is your optimal parameter vector

% let's see how much better it a does:
figure
subplot(1,2,1)
plot(x,y,'g*-')
hold on
plot(x,modfunc1(x,ain),'k')
plot(x,modfunc1(x,a),'r')
%
subplot(1,2,2)
semilogy(x,y,'g*-')
hold on
semilogy(x,modfunc1(x,ain),'k')
semilogy(x,modfunc1(x,a),'r')
legend('data','first guess','optimized model')

% kvg is a flag to say whether convergence was achieved before the routine
% gave up
% iter is the number of iterations that was used
% covp is the covariance matrix of the coefficients (the square roots of
% the diagonal elements are the uncertainties (the sigmas) of the coefficients) 
% r2 is the overall correlation coefficient

coefficients = a'
sigmas = sqrt(diag(covp))'
R2 = r2

% coefficients = 0.6025    0.0733    0.3995
% sigmas =    0.1506    0.0255    0.0293
% R2 =  0.7862

% so mu = 0.40 +- 0.03 per day
% 69% o the variability in the data is explained by our model
% what about P1 and P2?
% do we really need P1?

% Redo the fit with modified function without P1
