% Central limit theorem example using Bernoulli RVs

% Input parameters

p = 0.3;
N = 20;
Ntrials = 100;

X = randB(N,Ntrials,p);  % Generate the RVs
Xav = mean(X);
dxbin = 1/N;
xbin = 0:dxbin:1;
freq = hist(Xav,xbin);

subplot(2,2,1)
bar(xbin,freq)
xlim([0 1])
xlabel(['Xav, ' int2str(N) ' Bernoulli RVs with p =' num2str(p)])
ylabel(['Frequency (out of ' int2str(Ntrials) ' trials)'])
hold on
mu = p;
sigma = sqrt(p*(1-p)/N);
PDF = exp(-(xbin-mu).^2/(2*sigma^2))/(sigma*sqrt(2*pi));
plot(xbin,PDF*dxbin*Ntrials,'r-')
legend('Xav','Normal')
title('Frequency distribution')
hold off

subplot(2,2,2)
cumfrac = cumsum(freq)/Ntrials;
bar(xbin,cumfrac,1)  % Use bars with width=1 to get continuous fill
xlim([0 1])
ylabel('Cumulative fraction Xav <= x')
xlabel('x')
hold on
CDF = 0.5*(1+erf((xbin - mu)/(sigma*sqrt(2))));
plot(xbin,CDF,'r-')
legend('Avg. of RVs','Normal','Location','SouthEast')
title('CDF')
