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

%% 2.1 Analysing distributions: some statistical ideas
% Reproduce Fig. 2.1 using the |randn| and |histogram| functions, discussed
% in the notes for Ch. 1. Judging from Fig. 2.1(d), the distribution has
% mean $\bar{x} \approx 10$ and standard deviation $\sigma_x\approx 1$.

x0 = 10;
sigma = 1;

%rng(1);                         % Initialize the RNG
x = x0 + sigma*randn(1000,1);   % Generate all 1000 points at once
disp(x(1:5))                    % Show the first 5

%%
%
hist(x(1:5),(8:0.25:12))   % The second argument here sets the bins -
                                %  each bin is 0.25 units wide, and they
                                %  span the range 8 to 12.

%%
% Note that the data are different from those shown in Hughes and Hase -
% they're random, after all. Now format the figure so that it is closer to
% what the book shows.

hist(x(1:5),(8:0.25:12))
box off                         % Remove the box
xlabel('Variable, x')           % Label the x-axis
ylabel('Occurrence')            % Label the y-axis
text(0.95,0.95,'(a)',...        % Place the figure label '(a)' in the upper
    'Units','Normalized')       %  right corner; use '...' to continue the
                                %  command onto the next line. Set the text
                                %  units to 'Normalized' so that the label
                                %  position is placed with respect to the
                                %  graph axes.

%%
hist(x(1:50),(8:0.25:12))
box off
xlabel('Variable, x')
ylabel('Occurrence')
text(0.95,0.95,'(b)',...
    'Units','Normalized')

%%
hist(x(1:100),(8:0.25:12))
box off
xlabel('Variable, x')
ylabel('Occurrence')
text(0.95,0.95,'(c)',...
    'Units','Normalized')

%%
hist(x,(8:0.25:12))
box off
xlabel('Variable, x')
ylabel('Occurrence')
text(0.95,0.95,'(d)',...
    'Units','Normalized')

%% 2.2 The mean
% The MATLAB function <https://www.mathworks.com/help/matlab/ref/mean.html
% |mean|> computes the mean. Remember that we set the |x0 = 10| at the
% beginning.

xbar = mean(x);
disp(xbar)

%%
% Compare this with Eq. (2.1) in Hughes and Hase. 
disp( sum(x)/length(x) )

%% 2.3 The width of the distribution: estimating the precision
% Consider the data in Table 2.1.

T = [10.0, 9.4, 9.8, 9.6, 10.5, 9.8, 10.3, 10.2, 10.4, 9.3];

%%
% Plot a histogram of the data. Let MATLAB choose the bins automatically.

hist(T)
xlabel('Period (s)')
ylabel('Occurrence')

%%
% Try using six bins.

hist(T,6)
xlabel('Period (s)')
ylabel('Occurrence')

%%
% Display the mean.
Tmean = mean(T);
disp(Tmean)

%% 2.3.1 Rough-and-ready approach to estimating the width
% To find the full range of the data, use the functions |max| and |min|.

Tmax = max(T);
Tmin = min(T);
dmax = abs(Tmean - [Tmin, Tmax]);

disp(dmax)

%%
% The _total spread_ of the data is (9.9 � 0.6) s, as given in Hughes and
% Hase (note the significant figures). The 'rough-and-ready' method to
% estimating the spread gives:

sigmaRough = (2/3)*max(dmax);
disp(sigmaRough)

%%
% Compare this to the standard deviation, using the MATLAB function |std|.

sigma = std(T);
disp(sigma)

%%
% Check that this agrees with (2.3).

nT = length(T);
sigmaAlt = sqrt(sum((T-mean(T)).^2)/(nT-1));
disp(sigmaAlt)

%%
% The 'rough-and-ready' estimate is only useful for small sample sizes, and
% even then it's probably better just to calculate the standard deviation.
% The estimate gets worse as the sample size increases.

roughWidth = zeros(4,1);    % Initialize arrays
stdDev = zeros(4,1);
nSample = [3; 10; 30; 100; 300; 1000];

for iSample =  1:length(nSample)
    sample = randn(nSample(iSample),1);
    dMax = abs([max(sample),min(sample)] - mean(sample));
    roughWidth(iSample) = (2/3)*max(dMax);
    stdDev(iSample) = std(sample);
    disp(['N = ' num2str(nSample(iSample))])
    disp(['Rough method: ' num2str(roughWidth(iSample))])
    disp(['Standard dev.: ' num2str(stdDev(iSample))])
    disp(' ')
end

%% 2.4 Continuous distributions
% To reproduce Fig. 2.3, use the RNG to make 10,000 normally distributed
% numbers with the same mean and standard deviation as we used to reproduce
% Fig. 2.1.
x0 = 10;
sigma = 1;
nSample = 10000;

x = x0 + sigma*randn(nSample,1);   % Generate all 1000 points at once

%%
% The bins in Fig. 2.3(a) are 1.6 wide; the width decreases by a factor of
% two in each successive histogram. Use a loop to generate the four
% plots, using a character array to produce the different labels. To
% superpose the normal distribution on each plot, use the MATLAB function
% |linspace| to make an array |xVar|, and use (2.4); normalize to match the
% maximum of the histogram.
binStart = 6;
binEnd = 14;

xVar = linspace(binStart,binEnd);   % Default size of LINSPACE is [1,100]
f = (1/(sigma*sqrt(2*pi)))*exp(-((xVar - x0).^2/(2*sigma^2)));

binWidth = [1.6, 0.8, 0.4, 0.2];
labelArray = 'abcd';

for iPlot = 1:4
    
    figure;                 % Make a new figure for each panel
    h = hist(x, (binStart : binWidth(iPlot) : binEnd));
    hMax = max(h.Values);   % h.Values contains the bin values
    box off
    
    hold on                 % Overlay next plot on current plot
                            % Plot the distribution as a black line
    plot(xVar,f*hMax*sqrt(2*pi)*sigma,'k-');
    
    xlabel('Variable, x')
    ylabel('Occurrence')
    label = ['(' labelArray(iPlot) ')'];
    text(0.95,0.95,label,...
        'Units','Normalized')
    
    hold off                % Return plot to default behaviour
end

%% 2.5 The normal distribution
% Reproduce Fig. 2.4. Define an
% <http://www.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html
% _anonymous function_> for (2.4) to simplify the computation of the three
% distributions. Note the LaTeX markup for the legend.

xVar = linspace(0,20);
f = @(x,mu,sigma) (1/(sigma*sqrt(2*pi)))*exp(-((x - mu).^2/(2*sigma^2)));

plot(xVar, f(xVar, 10, 1), 'k-', ...    % sigma = 1, solid black line
    xVar, f(xVar, 10, 2), 'k:', ...     % sigma = 2, dotted black line
    xVar, f(xVar, 10, 3), 'k--');       % sigma = 3, dashed black line

box off
xlabel('Random variable, x')
ylabel('Probability distribution, f(x)')
legend('\sigma = 1','\sigma = 2','\sigma = 3')

%% 2.6 Sample and parent distribution
% The default behaviour for
% <http://www.mathworks.com/help/matlab/ref/std.html |std|> is to compute
% the sample standard deviation using (2.6), but it accepts an optional
% parameter that causes it to use (2.7) instead.

sigmaParent = std(T,1);
disp(sigmaParent)

%%
% Check that this agrees with (2.5).

sigmaParentAlt = sqrt(sum((T-mean(T)).^2)/nT);
disp(sigmaParentAlt)

%%
% The sample standard deviation from Sec. 2.3.2 is a bit bigger.
disp(std(T))

%%
% We can use the RNG to compare (2.5) and (2.6). Make a |5 x 1000| array
% of normally distributed random numbers, and compute both versions of the
% standard deviation for each of the 10000 columns (with $N = 5$).

data = randn(5,1000);
stdSample = std(data);
stdParent = std(data,1);

disp(stdSample(1:5))
disp(stdParent(1:5))

%%
% The parent distribution for the RNG has $\sigma = 1$. Both |stdSample|
% and |stdParent| are functions of five random numbers, so they are also
% random numbers. Comparing their averages over all 10,000 columns, we find
% that (2.5) underestimates the parent standard deviation when $N = 5$; the
% $N-1$ in (2.6) corrects for this.

disp(mean(stdSample))
disp(mean(stdParent))

%%
% For large data sets, though, the difference between |N| and |N-1| is
% negligible: consider |x| from Fig. 2.3.

disp(length(x))
disp(std(x))
disp(std(x,1))

%% 2.7 The standard error
% The sample standard deviation is an estimate of the distribution width.
% The standard _error_ $\alpha$ estimates how well we know its mean value.
% If $N = 1$, then $\alpha = \sigma$, but our knowledge of $\mu$ improves
% with every new measurement, as Fig. 2.1 shows.

%%
% Return to the data in Table 2.1. For these ten measurements, we can
% estimate the mean of the parent distribution to be (9.9 � 0.1) s.

disp(Tmean)

Tstderr = std(T)/sqrt(nT);
disp(Tstderr)

%%
% Confirm (2.7) for the |5 x 1000| random numbers in the array |data|:
% take the mean of each column, then find the standard deviation of the
% mean.

dataMean = mean(data);
dataMeanUnc = std(dataMean);
dataStd = std(data);
disp(dataMeanUnc)
disp(mean(dataStd)/sqrt(5))

%% 2.7.1 The error in the error
% For the data in Table 2.1, the fractional error in the error is 20% (to
% one significant figure).

disp(1/sqrt(2*nT - 2))

%%
% Confirm (2.8) for the |5 x 1000| random numbers in the array |data|:
% take the standard deviation of each column, then find the standard
% deviation of the standard deviation.

dataStdUnc = std(dataStd);
disp(dataStdUnc)
disp(mean(dataStd)/sqrt(2*5 - 2))

%%
% $Id: P233HH2.m 150 2016-02-02 09:07:14Z physcmpt $