% Demo: The distribution of a squared normal (or sum of many).
% Matlab code.
%
% nfits=0: just show the histogram
% nfits=1: try to fit normal
% nfits=2: try to fit exponential
% nfits=3: show the real chi^2 distribution
%
function chidemo(k,nfits)

if nargin<1, k=1; end
if nargin<2, nfits=0; end

n = 500000;      % Do it this many times.

% Sample from standard normal distribution.
x = normrnd(0,1,k,n);

% Square them
xsq = x.^2;

% Sum each column of k numbers.
y = sum(xsq, 1);

% Descriptive statistics
fprintf('Y = sum of %d squared standard normals\n', k);
fprintf('Sample mean of Y    = %.3f\n', mean(y));
fprintf('Sample stddev. of Y = %.3f\n', std(y));

% Empirical histogram
clf
histogram(y, 'normalization', 'pdf', 'facecolor', 'c');
hold on

% Try to fit some distributions, and plot them on top
yy = linspace(0.01, max(y), 1000);

if nfits>=1
    [nmuhat,nshat] = normfit(y);
    fnorm = normpdf(yy, nmuhat, nshat);
    plot(yy, fnorm, 'linewidth', 2);
end

if nfits>=2
    expmuhat = expfit(y);
    fexp = exppdf(yy, expmuhat);
    plot(yy, fexp, 'linewidth', 2);
end

if nfits>=3
    fchi2 = chi2pdf(yy, k);
    plot(yy, fchi2, 'linewidth', 2);
end

title(sprintf('Distribution of Y = sum of %d squared normals', k));

L = {'empirical'};

if nfits>=1, L{end+1} = sprintf('fit N(%.1f, %.1f^2)', nmuhat, nshat); end
if nfits>=2, L{end+1} = sprintf('fit Exp(%.1f)', expmuhat); end
if nfits>=3, L{end+1} = sprintf('chi2(%d)', k); end
legend(L);