close all
clear all
display('john horel 2015')
display('demonstrate central limit theorem')
display(' ')
display(' ')
sides = input('enter the number of sides on the die (6 or 10 or...): ')
dice = input('enter the sample size (number of dice to roll at the same time): ')
rolls = input('enter the number of rolls: ')

%get random numbers between 0 and 1 
rd = rand(rolls,dice);
%compute probability that any one number comes up
%convert 0 - 1/sides = 1; 1/sides - 2/sides = 2, etc.
rdd = ceil(sides*rd);
% 

%define the population as all the possible rolls
pop = rdd(:);
%compute the sum and mean of the sample
if(dice>1)
ss = sum(rdd,2);
sm = mean(rdd,2);
else
ss = rdd;
sm = rdd;
end

%now plot
figure(1)
%plot population
subplot(4,1,1)
bar(pop)
axis('tight')
title('Population time series: each die roll is independent')
%plot population
subplot(4,1,2)
vals = [1:1:sides];
hist(pop,vals)
axis('tight')
title('Population histogram')
subplot(4,1,3)
%plot sample sum
valsd = [dice:1:sides*dice];
hist(ss,valsd)
axis('tight')
title('Sum of sample die rolls')
subplot(4,1,4)
%plot sample mean
%hist(sm,vals)
hist(sm)
axis('tight')
title('Mean of sample die rolls')

%compute the std of the sample means 
sstd = std(sm,0);
%compute std of population
pstd = std(pop,0);
%does the central limit theorem work? is the following close to sstd?
check = pstd/sqrt(dice);


figure(2)
%plot cdf of the population
subplot(3,1,1)
%plot cdf of sample sum
subplot(3,1,2)
%plot cdf of sample mean
subplot(3,1,3)

%compute the mean of the sample means
%compute the mean of the population


%determine whether the null hypothesis
% can be rejected for each sample at the 5% level

for i=1:rolls
    %use the t test for anb appropriate test for flood episodes
end
%plot
figure(3)
subplot(3,1,1)
%plot time series of sample means
subplot(3,1,2)
%plot time series of p values from the t test
subplot(3,1,3)
%plot h values from the t test
plot(h)

%compute the percentage fraction for which the null
%hypothesis can be rejected




