%
%  math495/stat490 -- 10 sept 03 -- djm
%
%  w02pairs.m:  

%  clear workspace & plots
clear;  close all;

%  experimental parameters
Nrolls = 36;  Ndata = 100;

%  initialize random number generator
phone = 2914814;
rand('state',phone);  %  initialize with my office phone number

%  collect pair stats on Ndata sets of Nroll
data = zeros(1,Ndata);
for jj = 1:Ndata
    pairs = 0;

    for kk = 1:Nrolls
        %  choose two uniformly random integers (1-6)
        die1 = ceil(6*rand(1,1));
        die2 = ceil(6*rand(1,1));
        
        % count if a pair
        pairs = pairs + (die1==die2);            
    end

    %  save the observed prob of pairs
    data(jj) = pairs/Nrolls;
end

%  plot a histogram
figure(1);  clf;  hold on
hist(data,21)

title(['\bf histogram of pair probability in ' num2str(Nrolls) ...
       ' rolls (based on ' num2str(Ndata) ' samples)'])
xlabel('\bf observed pair probability')
ylabel(['\bf frequency of occurrences in ' num2str(Nrolls) ' rolls'])

%  some statistics

disp(['mean prob = ' num2str(mean(data))])
disp(['variance  = ' num2str( var(data))])