function testStats(n,sig,xbar)
% testStats  Study behavior of distributions created by randn function
%            Generate samples from randn and compare to predictions for
%            normal distribution
%
% Synopsis:  testStats
%            testStats(n)
%            testStats(n,sig)
%            testStats(n,sig,xbar)
%
% Inputs:  n = number of samples
%          sig = target standard deviation of the samples
%          xbar = target mean of the samples

if nargin<1, n=500;      end
if nargin<2, sig=1.5;    end
if nargin<3, xbar=12.2;  end

% -- Generate samples, and compute actual mean and standard deviation
x = xbar + sig*randn(n,1);
xave = mean(x);
sdev = std(x);

n1 = 0;  %  Prepare to count values greater than one stdev from mean

for i = 1:length(x)

  % -- How far is this point from the mean of the sample
  dx = abs(x(i)-xave);     %  abs(...) is important!
  
  % -- Add to count if x(i) is greater than sig from the mean
  if dx > sig
    n1 = n1 + 1;
  end
  
end

% -- Plot a histogram, and save handle to internal data structure
h = histogram(x);

% -- Get maximum of counts in any bin to help draw vertical dashed lines
maxCount = max(h.Values);   %  or max( get(h,'Values') );

% -- Decorate the plot
hold('on')
plot(xave*[1,1],[0 maxCount],'k-','Linewidth',2)
plot(xave-sdev*[1,1],[0 maxCount],'r--',...
     xave+sdev*[1,1],[0 maxCount],'r--','LineWidth',3)
text(xave+1.08*sdev,0.95*maxCount,'\sigma','FontSize',22)
hold('off')

% -- Print a summary
fprintf('\n%d total samples\n',n);
fprintf('\nxbar, mean(x), xbar-mean(x) = %7.4f %7.4f  %12.3e\n',...
         xbar,xave,xbar-xave)
fprintf('sig, std(x), sig-std(x)       = %7.4f %7.4f  %12.3e\n',...
         sig,sdev,sig-sdev)
fprintf('\nn1 = %8d   n1/n = %7.4f\n',n1,n1/n);

end
