% plotAvgHistograms.m
% Plots histograms of averages of random variables, 
% and calculates their means and standard deviations.
% Jose M Soler. Nov.2010

% Set values and probabilities
clear all
x = [0,1];        % head/tail of a coin
p = [0.5,0.5];    % probabilities
% x = (1:6);        % dice values
% p = ones(1,6)/6;  % probabilities

% Find mean and standard deviation of variable x
x0 = meanValue(x,p);
dx1 = standDev(x,p);

% Set a mesh of values between xmin and xmax
xMesh = linspace(min(x),max(x),100);

% Set sum parameters
nVals = [16 64 256 1024];    % number of values to be tested
nv = numel(nVals);           % size of vector nVals
nAvgs = 1000;                % number of averages

% Find random averages
xAvg = zeros(nv,nAvgs);      % average of nSums, for each value of nTerms
for iv = 1:nv
    xAvg(iv,:) = randomAvgs(x,p,nVals(iv),nAvgs);
end

% Plot histograms of averages
for iv = 1:nv
    dx = dx1/sqrt(nVals(iv));      % expected stand. dev. of average of nVals
    y = exp(-(xMesh-x0).^2/(2*dx^2))/sqrt(2*pi*dx^2);  % expected normal distr.
    subplot(2,2,iv)
    [nh,xh] = hist(xAvg(iv,:));    % actual histogram of averages
    dxh = xh(2)-xh(1);             % interval of histogram divisions
    bar(xh,nh), hold on            % bar plot
    plot( xMesh, nAvgs*dxh*y, 'r' ), hold off   % compare with normal distr
    set(gca,'FontSize',15)
    axis( [0,1, 0,nAvgs/3] )
    xlabel('<x>')
    ylabel('Number of averages')
    title([num2str(nVals(iv)),' values'])
end

% Find averages of averages, and their standard deviations
for iv = 1:nv
    meanX(iv) = meanValue(xAvg(iv,:));
    deltaX(iv) = standDev(xAvg(iv,:));
    theoryDeltaX(iv) = standDev(x,p) / sqrt(nVals(iv));
end

% Print averages and standard deviations
nVals
meanX
deltaX
theoryDeltaX
