%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All rights reserved by Krishna Pillai, http://www.dsplog.com 
% The file may not be re-distributed without explicit authorization
% from Krishna Pillai.
% Checked for proper operation with Octave Version 3.0.0
% Author	: Krishna Pillai
% Email		: krishna@dsplog.com
% Version	: 1.0
% Date		: 17 July 2008
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Matlab/Octave script for plotting the probability density
% function of chi square random variable

close all
clear all
N = 10^6;

z1 = randn(1,N).^2;
z2 = randn(1,N).^2 +  randn(1,N).^2;
z3 = randn(1,N).^2 +  randn(1,N).^2 +  randn(1,N).^2;
z4 = randn(1,N).^2 +  randn(1,N).^2 +  randn(1,N).^2 +  randn(1,N).^2;
z5 = randn(1,N).^2 +  randn(1,N).^2 +  randn(1,N).^2 +  randn(1,N).^2 + randn(1,N).^2;

% probability density function of abs(z)
zBin = [0:0.01:14];
sigma2 = 1;
pzTheory = (1/(sigma2*2)).*exp(-(zBin)/(2*sigma2)); % theory for m=2
[nzSim1 zBinSim1] = hist(z1,zBin); % simulation
[nzSim2 zBinSim2] = hist(z2,zBin); % simulation
[nzSim3 zBinSim3] = hist(z3,zBin); % simulation
[nzSim4 zBinSim4] = hist(z4,zBin); % simulation
[nzSim5 zBinSim5] = hist(z5,zBin); % simulation


figure
plot(zBinSim1,nzSim1/(N*0.01),'b','LineWidth',2);
hold on
plot(zBinSim2,nzSim2/(N*0.01),'y','LineWidth',2);
plot(zBinSim3,nzSim3/(N*0.01),'c','LineWidth',2);
plot(zBinSim4,nzSim4/(N*0.01),'r','LineWidth',2);
plot(zBinSim5,nzSim5/(N*0.01),'g','LineWidth',2);
plot(zBin,pzTheory,'k.-')
xlabel('z');
ylabel('probability density, p(z)');
legend('m=1','m=2','m=3','m=4','m=5','theory (m=2)');
title('Probability density function');
axis([0 14 0 0.7]);
grid on