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

%  clear workspace & plots
clear;  close all;

%  experimental parameters
Ndata = 1000;

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

done = 0;  choice = 0;
%  menu
while (done==0)
    %  collect Ndata random values 
    data = zeros(1,Ndata);
    
    switch choice
        case{1}  % uniform
            n = input('n = # of possible values = ');
            data = ceil(n*rand(Ndata,1));
            bins = 1:n;
            theory = Ndata * (1/n);
            
        case{2}  % bernoulli
            n = input('n = # of trials = ');
            p = input('p = probability of 1 = ');
            data = sum(double(rand(n,Ndata)<p));
            bins = 0:n;
            choose = gamma(n+1) ./ gamma(bins+1) ./ gamma(n-bins+1);
            theory = Ndata * choose .* (p.^bins) .* ((1-p).^(n-bins)); 
            
        case{3}  % geometric
            p = input('p = probability of 1 = ');
            for jj = 1:Ndata
                count = 0;  val = 1;
                while(val>p)
                    val = rand(1,1);
                    count = count + 1;
                end

                data(jj) = count;
            end
            bins = 1:max(data);
            theory = Ndata * p * (1-p).^(bins-1);
            
        case{4}  % poisson
            lam = input('lambda = ');
            for jj = 1:Ndata
                count = -1;  test = exp(-lam);  val = 1;
                while(val>test)
                    val = val*rand(1,1);
                    count = count + 1;
                end

                data(jj) = count;
            end
            bins = 0:max(data);
            factorials = gamma(bins+1);
            theory = Ndata * test * ((lam.^bins)./factorials);
            
        case{5}  %  random prescribed      
            n = input('# of possible values = ');
            dpf = rand(1,n);  dpf = dpf ./ sum(dpf);
            cdf = cumsum(dpf);
            for jj = 1:Ndata
                die  = ((cdf-rand(1,1))>=0);
                [junk index] = max(die);

                data(jj) = index;
            end
            bins = 1:n;
            theory = Ndata * dpf;
    end

    %  statistics
    switch choice
        case{1,2,3,4,5}
            %  histogram of data
            clf;  hist(data,bins);  hold on
            plot(bins,theory,'r*-')
            
            title(['\bf histogram of outcomes'])
            xlabel('\bf random variable (X)')
            ylabel(['\bf frequency of occurrences in ' num2str(Ndata)])
    
            %  some mean & var
            disp(' ')
            disp(['mean prob = ' num2str(mean(data))])
            disp(['variance  = ' num2str( var(data))])
            disp(' ')
    end
        
    %  selection prompt
    disp('1)  uniform')
    disp('2)  bernoulli trials')
    disp('3)  geometric')
    disp('4)  poisson')
    disp('5)  random prescribed')
    disp('R)  done')
    choice = input('    choice = ');
    if(isempty(choice))
        done = 1;
    end
end