%  ### EXgaussian1.m ###           9.24.08
% Create histogram plot of the distribution of N averages of M (uniformly-distributed)
% random #s to demonstrate the emergence of the Gaussian (then performs a
% nonlinear regression to fit a Gaussian function to such; see EXregression5.m)
clear
% -----------
M= 1000;        % # of (uniformly distributed) random #s to average
N=1000;         % # of repeats (i.e., how many averages to compute) for histogram
binN= 20;      % # of bins for histogram
NLregress= 1;   % determine best fitting Gaussian function [0-no,1-yes]
% -----------
figure(1); clf; hold on; grid on;
% +++
% loop thru to compute the N averages (each loop deals with the M random #s)
for nn=1:N
    xR= rand(M,1);  % determine array of M random #s
    mu(nn)= mean(xR);   % compute/store mean value
end
% +++
[jj,kk]=hist(mu,binN);   % detrmine histogram distribution
bar(kk,jj);             % plot the histogram (as a bar plot)
xlabel('Mean value'); ylabel('# of occurences'); 
% +++
% perform a nonlinear regression to fit Gaussian function (as in
% EXregression5.m, this requires several external functions)
if NLregress==1
    freeList = {'AA','BB','CC'};    % fit three params.
    initP.AA = max(jj);   % provide some rough starting point for each parameters
    initP.BB = 0.01; 
    initP.CC = 0.5;
    figure(2); clf
    plot_handle = plot(kk,jj,'b-','LineWidth',2);    hold on
    plot(kk,jj,'r.');
    text_handle = text(mean(kk),max(jj)+1,'','HorizontalAlignment','center','FontSize',14);
    xlabel('X');ylabel('Y');
    [bestP,err] = fit('gaussianFunc',initP,freeList,kk,jj,plot_handle,text_handle);
    disp(['Best fit: amplitude= ',num2str(bestP.AA),', STD= ',num2str(bestP.BB),', mean value= ',num2str(bestP.CC)])
end