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

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

%  experimental parameters (try Nrvs = 1,6,1000)
Nrvs = 500;  mean_rv = 0.5;  db = mean_rv/8;  dx = 0.01;%dx = db/10;

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

%  random data:  uniform on (0,1)
rvs = rand(Nrvs,1);
xt = [-db/2 0 0 1 1 1+db/2];  theory = [0 0 1 1 0 0];

plot(xt,theory,'b');  hold on

choice = 2;
switch choice
    case {1}
        %  histogram points
        xk = -db/2:dx:1+db/2;  myhist = zeros(size(xk));
        
        %  build histogram
        for j=1:Nrvs
            s = (db/2-abs(xk-rvs(j))>0);
            myhist = myhist + s/sum(s)/dx;
        end
        title(['\bf estimated PDF via box histogram'])
    case {2}
        %  histogram points
        xk = -mean_rv:dx:1+mean_rv;  myhist = zeros(size(xk));
        
        %  build histogram
        for j=1:Nrvs
            s = exp(-(xk-rvs(j)).^2 / (0.5*db^2)) .* (abs(xk-rvs(j))<mean_rv);
            myhist = myhist + s/sum(s)/dx;
        end
        title(['\bf estimated PDF via gaussian histogram'])
end

yk = (xk>=0).*(1>=xk);

plot(xk,myhist/Nrvs,'k.')
plot(xk,myhist/Nrvs,'r')
plot(xk,yk,'kx')

error = sqrt(sum((yk-myhist/Nrvs).^2)/size(yk,2))

xlabel(['\bf x-axis'])
text(0,1.3,['\bf (\Delta b = ' num2str(db) ' , \Delta x = ' num2str(dx) ')'])
ylabel(['\bf ePDF (' num2str(Nrvs) ' rvs)'])
%text(0.2,0.2,'\bf note:  x_k''s extend beyond (0,1) - why?')