% demo sampling distn of N Bernoulli trials. Barnett 5/23/13
clear
N = 1e6;  % trials 1e3 or 1e6
R = 1e3;   % repeats of expt (for N=1e6 make 1e3)
S = zeros(R,1);
for r=1:R, S(r) = sum(rand(1,N)>0.5); end  % sum of Bernoulli samples.
mu = N/2; sig = sqrt(N/4);  % predicted normal
figure; h = ceil(sig/5); b = ceil(mu)+h*(-30:30);   % good bin width, centers
d = hist(S,b); bar(b,d); % get histogram
hold on; x=mu+sig*(-4:.01:4); % where to plot curve
plot(x,R*h*exp(-(x-mu).^2/(2*sig^2))/sqrt(2*pi)/sig,'r-','linewidth',5);
title(sprintf('sampling distribution for sum of N=%d fair trials in {0,1}',N))
