% test of the central limit theorem

N = 500; % number of samples to draw

K = 5000;  % number of trials
J = 3; %Number integers to sample
I = [1:J];
x = sum(randi(J,N,K))/N;  % draw integers 1 to J and find the average 

mn = sum(I)/J;  %Calculate the theoretical mean
vr = sum(I.^2)/J-mn^2;  % calculate the theoretical variance
xm = max(abs(x-mn));  % gives a range for the histogram
dx = 2*xm/30;  % number of grid points for the histogram
xv = [mn-xm:dx:mn+xm];  % range for the histogram
h = hist(x,xv)/(K*dx);  % normalized histogram
figure(1)
bar(xv,h)
hold on
y = [1:.01:J];
f = sqrt(N)* exp(-N*(y-mn).^2/(2*vr))/(sqrt(2*pi*vr));  % normal distribution
 
plot(y,f,'--','linewidth',2)

hold off

 