%
%  math495/stat490 -- 30 aug 03 -- djm
%
%  w01dice.m:  

%  clear workspace & plots
clear;  close all;

%  experimental parameters
Ndice = 2;  Nsides = 6;  Nrolls = 10*(Nsides^2);

%  initialize random number generator:  help "rand"
rand('state',0);  %  comment out for "real" randomness

%  initialize identical dice:  "help zeros"
prob = zeros(Ndice,Nsides);

for jj=1:Ndice
    vals(jj,:) =  1:Nsides ;
    prob(jj,:) = (1:Nsides)/Nsides;
%    prob(jj,:) = [1 2 3 4 5 6]/6;
end

%  collect stats on Nrolls: "help if", "help max"
data = zeros(1,Nrolls);
for kk = 1:Nrolls
    total = 0;
    for jj = 1:Ndice
        roll = 0;
        while (roll==0)
            %  choose a random number
            roll = rand(1,1);
        end
        %  obtain the die roll -- this is the only tricky part
        die  = ((prob(jj,:)-roll)>=0);
        [junk index] = max(die);
        total = total + vals(jj,index);
    end
    %  save the 2-dice total
    data(kk) = total;
end

%  plot a histogram:  "help hist", "help bar"
lo = Ndice-1;  hi = Nsides*Ndice+1;  vect = lo:1:hi;

figure(1);  clf;  hold on
stats = hist(data,vect)
bar(vect,stats,'w');

%  some titles!!
title(['\bf histogram of ' num2str(Ndice) '-dice totals'])
xlabel('\bf totals')
ylabel(['\bf # occurrences in ' num2str(Nrolls) ' rolls'])

%  a bit of 2-dice theory (s2.1, Ross):  "help plot"
if(Ndice==2)
    theory = [0:1:Nsides,Nsides-1:-1:0]*(Nrolls/Nsides^2)
    plot(vect,theory,'r*-')
    
    %  some useful results (#1.13, Ross)
    disp(['empirical probability of (7,11)   = ' ...
            num2str((stats(7)+stats(11))/Nrolls)])
    disp(['empirical probability of (2,3,12) = ' ...
            num2str((stats(2)+stats(3)+stats(12))/Nrolls)])
end