%
%  math495/stat490 -- 06 oct 03 -- djm
%
%  w06rbino.m:  

%  clear workspace & plots
clear;  %close all;

%  experimental parameters
Ndata = 2000;  data = zeros(Ndata,1);
n = input('n = # of trials = ');
p = 0.5;

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

% bernoulli trials loop
for jj = 1:Ndata            
    p        = rand(1,1);
    data(jj) = sum(double(rand(n,1)<p));
end

bins = 0:n;
choose = gamma(n+1) ./ gamma(bins+1) ./ gamma(n-bins+1);
theory = ones(size(bins))/(n+1);
            
%  histogram of data
figure(1)
clf;  hdata = hist(data,bins);  bar(bins,hdata/Ndata,'w');  hold on
plot(bins,theory,'r*-')
            
title(['\bf histogram of outcomes'])
xlabel('\bf random variable (X)')
ylabel(['\bf % frequency of occurrences in ' num2str(Ndata)])