%
%  math495/stat490 -- 24 sept 03 -- djm
%
%  w04exp.m:  

%  exponential random data

%  initialize settings
rand('state',7654321)
Nrvs = 5000;  lambda = 1;  x_mean = 1/lambda;

plottype = 0;

%  mapping random data from uniform
u_rv = rand(Nrvs,1);
x_rv = -log(1-u_rv)/lambda;

if(plottype==1|plottype==0)
    %  plot histogram
    figure(1);  clf
    db = x_mean/8;
    bins = db/2:db:ceil(max(x_rv)/db)*db+(db/2);
    xhist = hist(x_rv,bins);  bar(bins,xhist/(db*Nrvs),'w');  hold on
    top = 1.1*max(xhist/(db*Nrvs));
    axis([-db/2 3*x_mean 0 top])

    %  exact probability density function
    density = lambda*exp(-lambda*bins);
    plot(bins,density,'r*')
    plot([0 bins],[lambda density],'r')

    title(['\bf histogram of mapped random variables (\lambda = ' ...
            num2str(lambda) ')'])
    ylabel('\bf f(x), density (frequency/binsize)')
    xlabel('\bf x-axis')

    text(x_mean,top/2,'\bf exact density = \lambda e^{-\lambda x}')
end
    
if(plottype==2|plottype==0)
    %  empirical cdf
    figure(2);  clf
    x_sort = sort(x_rv);
    cdf0 = (0:Nrvs-1)/Nrvs;  cdf1 = (1:Nrvs)/Nrvs;

    plot(x_sort,cdf1,'k.');  hold on

    %  connect the dots
    edfx = [x_sort' ; x_sort'];
    edfy = [cdf0 ; cdf1];
    plot(edfx(:),edfy(:),'b')
    
    title(['\bf empirical CDF'])
    ylabel('\bf F(x), cdf')
    xlabel('\bf x-axis')

    %  exact cdf
    cdf = 1-exp(-lambda*x_sort);
    plot(x_sort,cdf,'r')
    axis([0 4*x_mean 0 1])
    
    text(x_mean,1/2,'\bf exact cdf = 1 - e^{-\lambda x}')
end
    
disp(['mean = ' num2str(mean(x_rv))])