%
%  math495/stat490 -- 08 oct 03 -- djm
%
%  w06pdf.m:  

%  clear workspace & plots
clear;  close all;
figure(1);  hold on

%  experimental parameters
Nrvs = 2000;  lambda = 1.0;

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

%  random data
choice = 1;  pdfx = 0:(1/lambda/20):(5/lambda);
switch choice
    case{1}
        rvs = -log(rand(Nrvs,1))/lambda;
        plot(pdfx,exp(-lambda*pdfx)*lambda,'r')
    case{2}
        rvs = sqrt(-log(rand(Nrvs,1)));
        plot(pdfx,2*pdfx.*exp(-pdfx.^2),'r')
end

%  histogram pdf (bins = bin centres)
db = (1/lambda)/5;  bins = db/2:db:(5/lambda)+3*db/2;
rhist = hist(rvs,bins);  plot(bins(1:end-1),rhist(1:end-1)/Nrvs/db,'ko')

Nshift = 4;  
for jj = 1:Nshift-1
    bshift = bins - jj*db/Nshift;
    rhist = hist(rvs,bshift);  
    %  fix 1st bin
    rhist(1)  = rhist(1)*(Nshift/(Nshift-jj));
    bshift(1) = (-jj*db/Nshift + db)/2;
    plot(bshift(1:end-1),rhist(1:end-1)/Nrvs/db,'kx')
end

title(['\bf histogram of outcomes & empirical cdf'])
xlabel('\bf random variable (X)')
ylabel(['\bf % frequency of occurrences in ' num2str(Nrvs)])
axis([-db/2 5/lambda 0 max(1/lambda,1)])

%  empirical cdf
if(1==1)
    rv_sort = sort(rvs);
    cdf0 = (0:Nrvs-1)/Nrvs;  cdf1 = (1:Nrvs)/Nrvs;

    plot(rv_sort,cdf1);  hold on
end