% Gaussian Demo

% Brian D. Jeffs
% 9/18/05

echo on
X = randn(100000,1);

plot(X(1:100),'+');
title('100 Gaussian samples');
pause

plot(X(1:1000),'+');
title('1000 Gaussian samples');
pause

sound(X,22000);

pause

Nbins = 10
Nsamps = 100
[pdf_hat, x] = hist(X(1:Nsamps),Nbins);
binWidth = x(3)-x(2);
pdf_hat = pdf_hat/Nsamps/binWidth;
xSamps = [-5:.1:5];
trueGaussian = normpdf(xSamps);
bar(x,pdf_hat);
hold on
plot(xSamps,trueGaussian,'r');
legend(['Histogram for ',num2str(Nsamps),'samples'],'True Gaussian pdf');
title('Comparison of histogram pdf estimate to true pdf');
hold off
pause

Nbins = 10
Nsamps = 1000
[pdf_hat, x] = hist(X(1:Nsamps),Nbins);
binWidth = x(3)-x(2);
pdf_hat = pdf_hat/Nsamps/binWidth;
xSamps = [-5:.1:5];
trueGaussian = normpdf(xSamps);
bar(x,pdf_hat);
hold on
plot(xSamps,trueGaussian,'r');
legend(['Histogram for ',num2str(Nsamps),'samples'],'True Gaussian pdf');
title('Comparison of histogram pdf estimate to true pdf');
hold off
pause

Nbins = 10
Nsamps = 100000
[pdf_hat, x] = hist(X(1:Nsamps),Nbins);
binWidth = x(3)-x(2);
pdf_hat = pdf_hat/Nsamps/binWidth;
xSamps = [-5:.1:5];
trueGaussian = normpdf(xSamps);
bar(x,pdf_hat);
hold on
plot(xSamps,trueGaussian,'r');
legend(['Histogram for ',num2str(Nsamps),'samples'],'True Gaussian pdf');
title('Comparison of histogram pdf estimate to true pdf');
hold off
pause

Nbins = 100
Nsamps = 100000
[pdf_hat, x] = hist(X(1:Nsamps),Nbins);
binWidth = x(3)-x(2);
pdf_hat = pdf_hat/Nsamps/binWidth;
xSamps = [-5:.1:5];
trueGaussian = normpdf(xSamps);
bar(x,pdf_hat);
hold on
plot(xSamps,trueGaussian,'r');
legend(['Histogram for ',num2str(Nsamps),'samples'],'True Gaussian pdf');
title('Comparison of histogram pdf estimate to true pdf');
hold off
pause

X = 8*(rand(100000,1)-.5);

plot(X(1:100),'+');
title('100 Uniform samples, a = -4, b = 4');
pause

plot(X(1:1000),'+');
title('1000 Uniform samples, a = -4, b =  4');
pause

sound(X,22000);

pause

Nbins = 10
Nsamps = 100
[pdf_hat, x] = hist(X(1:Nsamps),Nbins);
binWidth = x(3)-x(2);
pdf_hat = pdf_hat/Nsamps/binWidth;
xSamps = [-5:.1:5];
trueUniform = 1/8*(abs(xSamps)<=4);
bar(x,pdf_hat);
hold on
plot(xSamps,trueUniform,'r');
legend(['Histogram for ',num2str(Nsamps),' samples'],'True Uniform pdf');
title('Comparison of histogram pdf estimate to true pdf');
hold off
pause

Nbins = 20
Nsamps = 1000
[pdf_hat, x] = hist(X(1:Nsamps),Nbins);
binWidth = x(3)-x(2);
pdf_hat = pdf_hat/Nsamps/binWidth;
xSamps = [-5:.1:5];
trueUniform = 1/8*(abs(xSamps)<=4);
bar(x,pdf_hat);
hold on
plot(xSamps,trueUniform,'r');
legend(['Histogram for ',num2str(Nsamps),' samples'],'True Uniform pdf');
title('Comparison of histogram pdf estimate to true pdf');
hold off
pause

Nbins = 100
Nsamps = 100000
[pdf_hat, x] = hist(X(1:Nsamps),Nbins);
binWidth = x(3)-x(2);
pdf_hat = pdf_hat/Nsamps/binWidth;
xSamps = [-5:.1:5];
trueUniform = 1/8*(abs(xSamps)<=4);
bar(x,pdf_hat);
hold on
plot(xSamps,trueUniform,'r');
legend(['Histogram for ',num2str(Nsamps),' samples'],'True Uniform pdf');
title('Comparison of histogram pdf estimate to true pdf');
hold off
pause


X = 8*(rand(100000,10)-.5);
X = sum(X')';

plot(X(1:100),'+');
title('100 samples of sum of 10 uniform RVs, a = -4, b = 4');
pause

plot(X(1:1000),'+');
title('1000 samples of sum of 10 uniform RVs, a = -4, b =  4');
pause


Nbins = 100
Nsamps = 100000
[pdf_hat, x] = hist(X(1:Nsamps),Nbins);
binWidth = x(3)-x(2);
pdf_hat = pdf_hat/Nsamps/binWidth;
xSamps = [-40:.5:40];
% Note, variance of rand, i.e. U(0,1), is 1/12, std = 1/sqrt(12).  We scaled std by 8,
% in line 126 above, producing std(X)=4/sqrt(3), then summed 10 random 
% variates to get std(X)=sqrt(10)*4/sqrt(3) = 7.303
trueGaussian = normpdf(xSamps,0,7.303);
bar(x,pdf_hat);
hold on
plot(xSamps,trueGaussian,'r');
legend(['Histogram for ',num2str(Nsamps),' samples'],'True Gaussian pdf');
title('Comparison of histogram pdf estimate to true pdf for sum of 10 uniform RVs');
hold off
pause

% Sum of binomial samples
p = .3;
N = 10;
K = 10; %number to sum
X = binornd(N*ones(100000,K),p);
X = (X-N*p)./sqrt(N*p*(1-p));  % normalize to zero mean, unit var.
X = sum(X')'./sqrt(K);  % sum K binomials, normalize by sqrt(K).

plot(X(1:100),'+');
title(['100 samples of sum of 10 normalized binomial RVs, p = ',num2str(p)]);
pause

plot(X(1:1000),'+');
title(['1000 samples of sum of 10 normalized binomial RVs, p = ',num2str(p)]);
pause


Nbins = 200
Nsamps = 100000
[pdf_hat, x] = hist(X(1:Nsamps),Nbins);
binWidth = x(3)-x(2);
pdf_hat = pdf_hat/Nsamps/binWidth;
xSamps = [-5:.02:5];
trueGaussian = normpdf(xSamps,0,1);
plot(x,pdf_hat);
hold on
plot(xSamps,trueGaussian,'r');
legend(['Histogram for ',num2str(Nsamps),' samples'],'True Gaussian pdf');
title('Comparison of histogram pdf estimate to true pdf for sum of 10 binomial RVs');
hold off
pause

plot(x,cumsum(pdf_hat)./sum(pdf_hat),xSamps,cumsum(trueGaussian)./sum(trueGaussian),'r');
legend(['Sample PDF for ',num2str(Nsamps),' samples'],'True Gaussian PDF');
title('Comparison of sample PDF distribution estimate to true PDF for sum of 10 binomial RVs');
pause

% Sum of binomial samples, more terms in sum
p = .3;
N = 10;
K = 50; %number to sum
X = binornd(N*ones(100000,K),p);
X = (X-N*p)./sqrt(N*p*(1-p));  % normalize to zero mean, unit var.
X = sum(X')'./sqrt(K);  % sum K binomials, normalize by sqrt(K).

plot(X(1:100),'+');
title(['100 samples of sum of 10 normalized binomial RVs, p = ',num2str(p)]);
pause

plot(X(1:1000),'+');
title(['1000 samples of sum of 10 normalized binomial RVs, p = ',num2str(p)]);
pause


Nbins = 200
Nsamps = 100000
[pdf_hat, x] = hist(X(1:Nsamps),Nbins);
binWidth = x(3)-x(2);
pdf_hat = pdf_hat/Nsamps/binWidth;
xSamps = [-5:.02:5];
trueGaussian = normpdf(xSamps,0,1);
plot(x,pdf_hat);
hold on
plot(xSamps,trueGaussian,'r');
legend(['Histogram for ',num2str(Nsamps),' samples'],'True Gaussian pdf');
title('Comparison of histogram pdf estimate to true pdf for sum of 50 binomial RVs');
hold off
pause

plot(x,cumsum(pdf_hat)./sum(pdf_hat),xSamps,cumsum(trueGaussian)./sum(trueGaussian),'r');
legend(['Sample PDF for ',num2str(Nsamps),' samples'],'True Gaussian PDF');
title('Comparison of sample PDF distribution estimate to true PDF for sum of 10 binomial RVs');
pause

echo off