%% Central Limit Theorem Example

%Make up a non-Normal PDF.  

PDF=inline('1/3*((x>-3).*(x<=-2).*(2*(x+3))+(x>-2).*(x<=-1).*(-2*(x+1))+(x>0).*(x<=1).*(x)+(x>1).*(x<=2).*(-x+2))','x');

dx=0.1;
x=(-3:dx:2);

figure
plot(x,PDF(x),'linewidth',2)
xlabel('x','fontsize',14)
ylabel('PDF(x)','fontsize',14)

%% Use Accept-Reject Method to select from the PDF

N=10; %Max number of random variable to add together. 
M=1000; %number selected from the distribution  
R=zeros(M,N);

cnt=1;
while cnt<=(N*M)
    xtest=rand*5-3;
    ytest=rand;
    if ytest<=PDF(xtest)
        R(cnt)=xtest;
        cnt=cnt+1;
    end
end

%Check that our Accept-Reject Method is working  
[counts bin]=hist(R,(-3:dx:2))
pdfcurve=PDF(x);

%normalize to match our normalized PDF
normcounts=counts/sum(counts)/dx;
%check normalization
sum(normcounts*dx)

figure
bar(bin,normcounts,'r')
hold on
plot(bin,pdfcurve,'k','linewidth', 2)
xlabel('x','Fontsize',16)
ylabel('Normalized Counts','Fontsize', 16)
legend('Generated Data','PDF(x)')
 
%% Show the disribution of N random variables added together

n=10;    %number of distributions to combine
u=-1; %PDF mean
s2=13/6; 

Sn=sum(R(:,1:n),2);
Zn=(Sn-n*u)/sqrt(n*s2);

dx=0.1;
x=(-3:dx:3)
[counts bin]=hist(Zn,x)

%normalize to match our normalized PDF
normcounts=counts/sum(counts)/dx;
%check normalization
sum(normcounts*dx)

%compare with 'standard' Normal distribution (u=0, s2=1);
PDF_SN=inline('1/sqrt(2*pi)*exp( -x.^2/2)','x');

figure
bar(bin,normcounts,'r')
hold on
plot(bin,PDF_SN(bin),'b','linewidth', 2)
xlabel('x','Fontsize',16)
ylabel('Normalized Counts','Fontsize', 16)
legend('Generated Data','Standard Normal')
s=sprintf('Zn, n=%d',n)
title(s)



