%% Physics 233, MATLAB Exercises for Hughes and Hase: Chapter 3
% MATLAB activities to complement _Measurements and their Uncertainties_,
% Chapter 3, "Uncertainties as probabilities." Section numbers in the book
% correspond to those listed below.
%
% Author: J. S. Dodge, 2016

%% 3.1 Distributions and probability
% The MATLAB Statistics and Machine Learning Toolbox defines several
% common probability distributions. Consider the normal distribution.

mu = 0;
sigma = 1;
x = linspace(-5,5);
plot(x,normpdf(x,mu,sigma))

%%
% Confirm that it is properly normalized. The syntax |@(x)
% normpdf(x,mu,sigma)| defines |normpdf(x,mu,sigma)| as an anonymous
% function of |x|, with |mu = 0| and |sigma = 1| inherited from the earlier
% assignments.

checkNorm = quad( @(x) normpdf(x,mu,sigma), -Inf, Inf );
disp(checkNorm)

%%
% Confirm that 68% of the area lies in the range $x \in (-\sigma,+\sigma)$.

checkArea = quad( @(x) normpdf(x,mu,sigma), -sigma, sigma );
disp(checkArea)

%%
% Determine the first two moments. Analytically, the first moment is zero
% and the second moment is unity. The numerical result for the first moment
% is nonzero but insignificant.

firstMoment = quad( @(x) x.*normpdf(x,mu,sigma), -Inf, Inf );
disp(firstMoment)

%%
%
secondMoment = quad( @(x) x.^2.*normpdf(x,mu,sigma), -Inf, Inf );
disp(secondMoment)

%% 3.2 The Gaussian probability distribution function
%% 3.2.1 Probability calculations
% Hughes and Hase define the error function differently from MATLAB and
% other standard references. In MATLAB the function defined by (3.8) is
% called _the normal cumulative distribution function_, and is given by the
% function <http://www.mathworks.com/help/stats/normcdf.html |normcdf|>.
% The more common notation is
%
% $$ \Phi(x;\bar{x},\sigma) \equiv \frac{1}{\sigma\sqrt{2\pi}}
% \int_{-\infty}^x \exp\left[-\frac{(x-\bar{x})^2}{2\sigma^2}\right]\,dx$$
%
% As discussed in the MATLAB documentation and elsewhere, the _error
% function_ is defined (usually) as
%
% $$ \mathrm{erf}(x) \equiv \frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\,dt,$$
%
% and is related to $\Phi(x;\bar{x},\sigma)$ (or
% $\mathrm{Erf}(x;\bar{x},\sigma)$ in the notation of Hughes and Hase) by
%
% $$ \Phi(x;\bar{x},\sigma) = \frac{1}{2}
% \left[1 - \mathrm{erf}\left(\frac{-(x-\bar{x})}{\sigma\sqrt{2}}\right)\right]$$
%
% Check:

plot(x,normcdf(x,mu,sigma),'k-',...
    x,(1/2)*(1 - erf(-(x-mu)/(sigma*sqrt(2)))),'r:');

xlabel('Variable, x')
ylabel('Normal CDF')
legend('\Phi(x;\mu,\sigma)',...
    '\{1 - erf[-(x-\mu)/(\sigma\surd{2})]\}/2',...
    'Location','SouthEast')

%% 3.3 Confidence limits and error bars
% Check (3.10) using |normcdf|.

oneSigmaFrac = normcdf(mu+sigma,mu,sigma) - normcdf(mu-sigma,mu,sigma);
disp(oneSigmaFrac)

%%
% Alternatively, use |erf|. The symmetry $\textrm{erf}(-x) =
% -\textrm{erf}(x)$ allows simplification:
%
% $$ P = \Phi(\bar{x}+\sigma;\bar{x},\sigma) - 
% \Phi(\bar{x}-\sigma;\bar{x},\sigma) = 
% \textrm{erf}\left(\frac{1}{\surd{2}}\right)$$

disp(erf(1/sqrt(2)))

%% 3.3.1 Extended ranges
% Check Table 3.1.

rangeInSigma = [1, 1.65, 2, 2.58, 3];
fracPercent = erf(rangeInSigma/sqrt(2))*100;

disp([rangeInSigma;fracPercent])

%% 3.3.2 Rejecting outliers
%% 3.3.3 Experimental example of a Gaussian distribution
%% 3.3.4 Comparing experimental results with an accepted value

%% Poisson probability function for discrete events
% The MATLAB function <http://www.mathworks.com/help/stats/poisspdf.html
% poisspdf> returns the Poisson distribution function.

%%
% Reproduce Fig. 3.5.

N = (0:8);
Nbar = 1.5;
p = poisspdf(N,Nbar);

bar(N,p)
xlim([-0.5 8.5])
xlabel('Number of counts, N')
ylabel('Poisson PDF')
box off
text(0.95,0.95,'(a)',...
    'Units','Normalized')

%%

N = (0:40);
Nbar = 15;
p = poisspdf(N,Nbar);

bar(N,p)
xlim([-0.5 40.5])
xlabel('Number of counts, N')
ylabel('Poisson PDF')
box off
text(0.95,0.95,'(b)',...
    'Units','Normalized')

%% 3.4.1 Worked example - Poisson counts
%% 3.4.2 Error bars and confidence limits for Poisson statistics
%% 3.4.3 Approximations for higher means
%% 3.5 The central limit theorem
%% 3.5.1 Examples of the central limit theorem
% Reproduce the top row of Fig. 3.7. The left-most plot is the PDF for a
% uniform distribution over (0,1).

xVar = linspace(0,1);
f = ones(size(xVar));
plot(xVar,f);
xlabel('Variable, x')
ylabel('Probability Density, f(x)')

%%
% The center plot is a histogram of 1000 trials.

rand('seed',1)
%rng(1);
x = rand(1,1000);
hist(x,10)
xlabel('Variable, x')
ylabel('Occurrence')

%%
% The right-most plot is a histogram of 1000 trials, where the result of
% each trial is the mean of five numbers drawn from the uniform
% distribution. To simplify the computation, note that the |mean| function
% operates independently on the columns of an array.

x = rand(5,1000);
xBar = mean(x);
hist(xBar,20);
xlabel('Mean')
ylabel('Occurrence')

%%
% Let's look at how this evolves from |N=2| to |N=5|, using 10,000 trials
% to get better statistics.

for N = 2:5
    x = rand(N,10000);
    xBar = mean(x);
    figure;
    hist(xBar,20);
    xlabel('Mean')
    ylabel('Occurrence')
end


%%
% $Id: P233HH3.m 151 2016-02-02 18:38:44Z physcmpt $