%% This demonstrates the error on fit parameters.
% The likelihood is modeled by a multivariate gaussian to extract sigma


%% Parameter setup

A=20;
B=5;

x=(.5:.5:50)';
u=A+B*x;

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

%% An example noisy data set

sigma_noise=5;
y=u+sigma_noise*randn(size(x));

figure
plot(x,u_noise,'r.','linewidth',2)
hold on
plot(x,u,'linewidth',2)
xlabel('time (\mus)')
ylabel('measured value')
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 B_found

for ii=1:N_experiments    
    y=u+sigma_noise*randn(size(x));
    X0=[y(1) (y(end)-y(1))/(x(end)-x(1))]; %these are the starting guesses
    [out]=fminsearch(func,X0,options,x,y);
    A_found(ii)=out(1);
    B_found(ii)=out(2); 
end


%% Compare our results to the calculated standard error
N=size(x,1);

d2udA2=-N/sigma_noise^2;
d2udB2=-1/sigma_noise^2*sum(x.^2);
d2udAdB=-1/sigma_noise^2*sum(x);

clear SIGMAinv

SIGMAinv(1,1)=-d2udA2;
SIGMAinv(1,2)=-d2udAdB;
SIGMAinv(2,1)=-d2udAdB;
SIGMAinv(2,2)=-d2udB2;

SIGMA=inv(SIGMAinv);

varA=SIGMA(1,1);
varB=SIGMA(2,2);

stderrA=sqrt(varA);
stderrB=sqrt(varB);

stderrA_observed=std(A_found);
stderrB_observed=std(B_found);

fprintf('\n----------------------------------------------\n')
fprintf('Reported Standard Error of A: \t\t\t\t\t\t\t\t%g\n', stderrA)
fprintf('Fit (MLE) Standard Deviation of A from Repeat Trials:\t\t%g\n\n', stderrA_observed)

fprintf('Reported Standard Error of B: \t\t\t\t\t\t\t\t%g\n', stderrB)
fprintf('Fit (MLE) Standard Deviation of B from Repeat Trials:\t\t%g\n', stderrB_observed)
fprintf('----------------------------------------------\n')


