%% Demonstration of Standard Error of the Mean and Confidence Intervals


%A measurement has a statistical error such that measured values are 
%normally distibuted around  the actual value. 

%In general, mean and sigma are unknown and we won't use them after generating data
sigma=5;   
mu=10;

%generate N measured values

n=50;
data=sigma*randn(n,1)+mu;

%show our values

xbins=(mu-3*sigma:sigma/(2*n):mu+3*sigma)
figure
hist(data,xbins)
xlabel('measured value')
ylabel('counts')

%The best guess for the actual value is the mean

est_mu=mean(data)   %estimated actual value

%How good is our estimate?

%The standard error of the mean (SEM) is the standard deviation of the 
%error in the sample mean relative to the true mean.

SE_mu=std(data)/sqrt(n)     %standard error of the mean (uncorrected, biased)

SE_mu=(std(data)*sqrt(n/(n-1)))/sqrt(n);     %standard error of the mean (unbaised)

%A measured parameter is usually given as est_mu +/- SE_mu

%We can also calculate a 'confidence interval'.  
%If the experiment is repeated many times, 95% of the time the true value
%would fall in the 95% CI

%95% confidence interval is 
% est_mu +/- SE_mu*1.96

%Note: This formalism requires 'n' to be greater than ~ 20

%This is not the same as saying that the there is a 95% probably that the
%true value falls withing the CI for a single measurement. 

%Let's check this...

%% Confidence Interval check

sigma=5;
mu=10;

n=50;
Nloops=10000;
ObservationInInterval=0;

for nn=1:Nloops
 
    data=sigma*randn(n,1)+mu;
    est_mu=mean(data);   %estimated actual value
    SE_mu=std(data)/sqrt(n);     %standard error of the mean
    SE_mu=(std(data)*sqrt(n/(n-1)))/sqrt(n);     %standard error of the mean (unbaised)
    
    %95% confidence interval is
    lb=est_mu-SE_mu*1.96;
    ub=est_mu+SE_mu*1.96;
    
    if (mu>lb)&(mu<ub)
        ObservationInInterval=ObservationInInterval+1;
    end
    
end

%percentage of times true value was found in 95% CI
perc=ObservationInInterval/Nloops*100














