%% This demonstrates the variance of parameter meets the CRLB


% This code is taken and modified slightly from FitExample.m


%% Parameter setup

A=3.5;
B=7;
N=100;
sigma=10; %this is the sigma in normally distributed noise model

x=(1:N)';
y=A+B*x; %A and B are switched comparted to FitExample to match my notes

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


%% Show and example data set

y_noise=y+sigma*randn(size(y));
figure
plot(x,y,'r','linewidth',2)
hold on
plot(x,y_noise,'.','linewidth',2)
legend('Expectation Model','Noisy data')
xlabel('X (m)','fontsize',16)
ylabel('Y (m)','fontsize',16)
set(gca,'fontsize',14)



%% 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_noise=y+sigma*randn(size(y));
    X0=[1 1]; %these are the starting guesses
    [out]=fminsearch(func,X0,options,x,y_noise);
    A_found(ii)=out(2);
    B_found(ii)=out(1); %I switched these to match my notes
end

%% Compare our results to the CRLB

FI=1/sigma^2*[N sum(x); sum(x) sum(x.^2)]

invFI=inv(FI);

%The CRLB variances are the diagonal elements of invFI
varA_crlb=invFI(1,1);
varB_crlb=invFI(2,2);

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

fprintf('\n')
fprintf('CRLB predicted variance of A: %g\n', varA_crlb)
fprintf('found variance of A:          %g\n\n', varA_observed)

fprintf('CRLB predicted variance of B: %g\n', varB_crlb)
fprintf('found variance of B:          %g\n\n', varB_observed)




