%% This demonstrates the variance of parameter meets the CRLB


%% Parameter setup

A=100;
tau=2;
B=10;

t=(.01:.01:10);
u=A*exp(-t/tau)+B;

options=optimset;
func=@LL_ExpDecay; %this points to the function that will be minimized

%% An example noisy data set

u_noise=GenPoisson(u);
figure
plot(t,u_noise,'r','linewidth',2)
hold on
plot(t,u,'linewidth',2)
xlabel('time (\mus)')
ylabel('counts')
legend('noisy observation','model')


%% Create many data sets (observations) and estimate the parameters

N_experiments=1000; %(N_experiments would have to be infinite to reach the CRLB)

clear A_found
clear tau_found
clear B_found

for ii=1:N_experiments    
    u_noise=GenPoisson(u);
    X0=[u_noise(1) 5 u_noise(end)]; %these are the starting guesses
    [out]=fminsearch(func,X0,options,t,u_noise);
    A_found(ii)=out(1);
    tau_found(ii)=out(2);
    B_found(ii)=out(3); 
end


%% Compare our results to the CRLB
clear I

dudA=exp(-t/tau);
dudt=A*t/tau^2.*exp(-t/tau);
dudB=1;

I(1,1)=sum(dudA.*dudA./u);
I(1,2)=sum(dudA.*dudt./u);
I(1,3)=sum(dudA.*dudB./u);
I(2,1)=I(1,2);
I(2,2)=sum(dudt.*dudt./u);
I(2,3)=sum(dudt.*dudB./u);
I(3,1)=I(1,3);
I(3,2)=I(2,3);
I(3,3)=sum(dudB.*dudB./u);

CRLB=inv(I);
varA_crlb=CRLB(1,1);
vartau_crlb=CRLB(2,2);
varB_crlb=CRLB(3,3);

varA_observed=var(A_found);
vartau_observed=var(tau_found);
varB_observed=var(B_found);

fprintf('\n----------------------------------------------\n')
fprintf('CRLB predicted variance of A: %g\n', varA_crlb)
fprintf('found variance of A:          %g\n', varA_observed)
fprintf('precision of A:    %g\n', sqrt(varA_crlb))
fprintf('CRLB/Var: %g\n\n', varA_crlb/varA_observed)

fprintf('CRLB predicted variance of tau: %g\n', vartau_crlb)
fprintf('found variance of tau:          %g\n', vartau_observed)
fprintf('precision of tau:    %g\n', sqrt(vartau_crlb))
fprintf('CRLB/Var: %g\n\n', vartau_crlb/vartau_observed)

fprintf('CRLB predicted variance of B: %g\n', varB_crlb)
fprintf('found variance of B:          %g\n', varB_observed)
fprintf('precision of B:    %g\n', sqrt(varB_crlb))
fprintf('CRLB/Var: %g\n', varB_crlb/varB_observed)
fprintf('----------------------------------------------\n')


