%
%  math495/stat490 -- 24 sept 03 -- djm
%
%  w04exp2.m:  
%
%  - displays some of the problems inherent
%    in estimating pdf with numerical 
%    differentiation 
%    (set variable:  choice = 1,2,3 at bottom)

%  exponential random data

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

plottype = 0;

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

if(plottype==1|plottype==0)
    %  plot histogram
    figure(1);  clf;  subplot(2,2,4);
    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 = 2*bins.*exp(-bins.^2);
    plot(bins,density,'r*')
    plot([0 bins],[0 density],'r')

    title(['\bf histogram of x = sqrt(-log(1-u))'])
    ylabel('\bf f(x), density (freq/binsize)')
    xlabel('\bf x-axis')

    text(1.25*x_mean,0.75*top,'\bf exact density = 2 x e^{-x^2}')
end
    
if(plottype==2|plottype==0)
    %  empirical cdf
    subplot(2,1,1)
    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 (' num2str(Nrvs) ' points)'])
    ylabel('\bf F(x), cdf')
    xlabel('\bf x-axis')

    %  exact cdf
    cdf = 1-exp(-x_sort.^2);
    plot(x_sort,cdf,'r')
    axis([0 3*x_mean 0 1])
    
    text(x_mean,1/2,'\bf exact cdf = 1 - e^{-x^2}')

    subplot(2,2,3);
    plot(x_sort,cdf1,'k.');  hold on
    plot(edfx(:),edfy(:),'b')
    plot(x_sort,cdf,'r')

    title(['\bf empirical CDF (much expanded)'])
    ylabel('\bf F(x), cdf')
    xlabel('\bf x-axis')
end
    
disp(['mean = ' num2str(mean(x_rv))])

xdata = [0 x_sort'];
ydata = [0 cdf1];

choice = 1;
switch choice
    case{1}
        figure(3);  subplot(2,1,1);  k = 1;
        dx_dy = (xdata(1+2*k:end)-xdata(1:end-2*k))/(2*k/Nrvs);
    case{2}
        figure(3);  subplot(2,1,2);  k = 25;
        dx_dy = (xdata(1+2*k:end)-xdata(1:end-2*k))/(2*k/Nrvs);
    case{3}
        figure(4); clf 
        k = 25;  pdf = zeros(size(xdata(1+k:end-k)));
        Navg = 8;
        for jj = 0:Navg-1
            pdf = pdf ...
                + (xdata(1+jj+2*(k-jj):end-jj)-xdata(1+jj:end-2*(k-jj)-jj)) ...
                / (2*(k-jj)/Nrvs);
        end
        dx_dy = pdf / Navg;
end

slope = 1./dx_dy;
plot(xdata(1+k:end-k),slope,'k.');  hold on
plot([0 bins],[0 density],'r')

xdata_mag = abs(xdata); 
ln_raw = -0.25*log(xdata_mag); 
dy = exp(ln_raw); 

title(['\bf empirical PDF by numerical differentiation across ' ...
        num2str(2*k) ' steps'])
ylabel('\bf empirical f(x)')
xlabel('\bf empirical f(x)')

axis([-db/2 3*x_mean 0 1.2])