
% Fitting the Houston pop data using algebraic, log and nlsq


lw=2;

data=load('houston.dat');X=data(:,1);Y=data(:,2);
t0=min(X); 
X=X-t0;
t0=min(X); t1=max(X);
T=t0:(t1-t0)/1000:t1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plotting data
figure(1)
plot(X,Y,'o','LineWidth',lw);
title('Data')
set(gca,'LineWidth',1.5);


pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% algebraic fit
figure(2)
R1_3=log(Y(3)/Y(1))/(X(3)-X(1))
R1_5=log(Y(5)/Y(1))/(X(5)-X(1))
R1_9=log(Y(9)/Y(1))/(X(9)-X(1))
plot(X,Y,'o',T,Y(1)*exp(R1_3*(T-X(1))), ...
             T,Y(1)*exp(R1_5*(T-X(1))), ...
             T,Y(1)*exp(R1_9*(T-X(1))), ...
             'LineWidth',lw)
title('Algebraic fit')
legend('data','R1-3','R1-5','R1-9',2)
axis([t0 t1 0 max(Y)])
set(gca,'LineWidth',1.5);

pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NLSQ fit using R0
figure(3)
% optional arguments for the optimization routines
%options=optimset('TolFun',1e-10,'TolX',1e-10,'MaxIter',1e10);
options=optimset;

% define the exponential growth with N0 fixed
% p(1) = R0 = growth rate
ExpGrowN0 = inline('18.63*exp(p(1)*x)','p','x');

% define inital guess for the paramaters
R0=0.04; param0=[R0];

% apply a nonlinear square fit and store paramaters in pyeast
paramN0=lsqcurvefit(ExpGrowN0,param0,X,Y,[],[],options)

plot(X,Y,'o',T,ExpGrowN0(param0,T),T,ExpGrowN0(paramN0,T),'LineWidth',lw)
title('NLSQ fit using R0')
axis([t0 t1 0 max(Y)])
legend('data','guess','fit',2)

pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NLSQ fit using R0 and N0
figure(4)
% optional arguments for the optimization routines
%options=optimset('TolFun',1e-10,'TolX',1e-10,'MaxIter',1e10);
options=optimset;

% define the exponential growth
% p(1) = R0 = growth rate
% p(2) = N0 = inital population
ExpGrow = inline('p(2)*exp(p(1)*x)','p','x');

% define inital guess for the paramaters
R0=0.04; N0=10; param0=[R0,N0];

% apply a nonlinear square fit and store paramaters in pyeast
param=lsqcurvefit(ExpGrow,param0,X,Y,[],[],options)

plot(X,Y,'o',T,ExpGrow(param0,T),T,ExpGrow(param,T),'LineWidth',lw)
title('NLSQ fit using R0 and N0')
axis([t0 t1 0 max(Y)])
legend('data','guess','fit',2)

pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NLSQ fit to log(data) using R0 and N0
figure(5)
% optional arguments for the optimization routines
%options=optimset('TolFun',1e-10,'TolX',1e-10,'MaxIter',1e10);
options=optimset;

% define the exponential growth
% p(1) = R0 = growth rate
% p(2) = L0 = log(N0) = log(inital population)
ExpGrowLog = inline('p(2)+p(1)*x','p','x');

% define inital guess for the paramaters
R0=0.04; L0=2; param0=[R0,L0];

% apply a nonlinear square fit and store paramaters in pyeast
paramLog=lsqcurvefit(ExpGrowLog,param0,X,log(Y),[],[],options)

subplot(2,1,1)
plot(X,log(Y),'o',T,ExpGrowLog(param0,T),...
                  T,ExpGrowLog(paramLog,T),'LineWidth',lw)
title('NLSQ fit to log(data) using R0 and N0')
axis tight
legend('log(data)','guess','fit',2)

subplot(2,1,2)
plot(X,Y,'o',T,exp(ExpGrowLog(param0,T)),...
             T,exp(ExpGrowLog(paramLog,T)),'LineWidth',lw)
axis([t0 t1 0 max(Y)])
legend('data','guess','fit',2)

pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot all fits
figure(6)
lw=1;
plot(X,Y,'o',T,exp(ExpGrowLog(paramLog,T)),...
             T,ExpGrowN0(paramN0,T),...
             T,ExpGrow(param,T),...
             'LineWidth',lw)
title('Comparing all the fittings')
axis([t0 t1 0 max(Y)])
legend('data','log(data)','data(R0)','data(R0,N0)',2)


%print -djpeg  figure.jpg
