clear all;
close all hidden;


%%%%%%%%%%%%%%%
%% DGP
%%%%%%%%%%%%%%%
n = 1000;
t0 = 0;
t1 = 500;
sigma = 0.15;
time = linspace(t0,t1,n);
ydet = 0.7*exp(-time/10)+1.1*exp(-time/100)+1.4*exp(-time/1000);
y = ydet+sigma*randn(1,n);
plot(time,y,'.b')
drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stima modello con 1 componente
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
optimset('lsqnonlin')
opts=optimset('Display','Iter','TolFun',10^(-6),...
    'MaxFunEvals',10000,'MaxIter',10000,'TolX',10^(-20));
model1 = @(parms,yobs,t) yobs-(parms(1)*exp(-t/parms(2)));
x0 = [1 300];
lb = zeros(2,1);
[X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT]=lsqnonlin(@(x) model1(x,y,time),x0,lb,[],opts);
% RESNORM = squared 2-norm of the residual at X

figure(1)
hold on;
plot(time,y,'.b')
fplot(@(x) X(1)*exp(-x/X(2)),[t0,t1],'r')
title(['A=',num2str(X(1)),' B=',num2str(X(2))]);
hold off;
drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stima modello con 2 componenti
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model2 = @(parms,yobs,t) yobs-(parms(1)*exp(-t/parms(2))+...
    parms(3)*exp(-t/parms(4)));
x0 = [1 8 1 500];
lb = zeros(4,1);
[Xa,RESNORM,RESIDUAL,EXITFLAG,OUTPUT]=lsqnonlin(@(x) model2(x,y,time),x0,lb,[],opts);
figure(2)
hold on;
plot(time,y,'.b')
fplot(@(x) Xa(1)*exp(-x/Xa(2))+Xa(3)*exp(-x/Xa(4)),[t0,t1],'r')
title(['A=',num2str(Xa(1)),' B=',num2str(Xa(2)),...
    ' C=',num2str(Xa(3)),' D=',num2str(Xa(4))]);
hold off;
drawnow;

% LSQCURVEFIT
model2 = @(parms,t) parms(1)*exp(-t/parms(2))+parms(3)*exp(-t/parms(4));
[Xb,RESNORM,RESIDUAL,EXITFLAG,OUTPUT]=lsqcurvefit(model2,x0,time,y,lb,[],opts);
[Xa; Xb]

% NLINFIT - STATISTICS TOOLBOX
model2 = @(parms,t) parms(1)*exp(-t/parms(2))+parms(3)*exp(-t/parms(4));
x0 = [1 8 1 500];
[par,resid,Jac,CovMtx,mse] = nlinfit(time,y,model2,x0);
par_ci = nlparci(par,resid,'covar',CovMtx);
[par_ci(:,1) par' par_ci(:,2)]
sqrt(mse)
nlintool(time,y,model2,x0)
drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Stima modello con 3 componenti
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model3 = @(parms,t) parms(1)*exp(-t/parms(2))+...
    parms(3)*exp(-t/parms(4))+parms(5)*exp(-t/parms(6));
x0 = [1 8 1 500 2 100];
lb = zeros(6,1);
[X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT]=lsqcurvefit(model3,x0,time,y,lb,[],opts);
figure(3)
hold on;
plot(time,y,'.b')
fplot(@(x) X(1)*exp(-x/X(2))+X(3)*exp(-x/X(4))+X(5)*exp(-x/X(6)),[t0,t1],'r')
title(['A=',num2str(X(1)),' B=',num2str(X(2)),' C=',num2str(X(3)),...
    'D=',num2str(X(4)),' E=',num2str(X(5)),' F=',num2str(X(6))]);
hold off;
drawnow;

% NLINFIT - STATISTICS TOOLBOX
model3 = @(parms,t) parms(1)*exp(-t/parms(2))+...
    parms(3)*exp(-t/parms(4))+parms(5)*exp(-t/parms(6));
x0 = [1 8 1 500 2 100];
[par,resid,Jac,CovMtx,mse] = nlinfit(time,y,model3,x0);
par_ci = nlparci(par,resid,'covar',CovMtx);
[par_ci(:,1) par' par_ci(:,2)]
sqrt(mse)
nlintool(time,y,model3,x0)