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

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

%  experimental parameters
Nrvs = 1000;  smooth = (0.1)^2;  

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

%  random data
choice = 1;
switch choice
    case {1}
        rvs = rand(Nrvs,1);
        mean_rv = 0.5;  xmax = 5*mean_rv;
        x = [-mean_rv 0 0 1 1 xmax];  theory = [0 0 1 1 0 0];
    case {2}
        mean_rv = 0.5;  xmax = 5*mean_rv;
        lambda = 1/mean_rv;
        rvs = -log(rand(Nrvs,1))/lambda;
        x = [0:mean_rv/10:xmax];  theory = lambda*exp(-lambda*x);
        x = [-mean_rv 0 x];  theory = [0 0 theory];
    case {3}
        mean_rv = sqrt(pi/2);  xmax = 5*mean_rv;
        rvs = sqrt(-log(rand(Nrvs,1)));
        x = [0:mean_rv/10:xmax];  theory = 2*x.*exp(-x.^2);
        x = [-mean_rv x];  theory = [0 theory];
end
plot(x,theory,'b')

%  build kernel-based pdf
dx = mean_rv/20;  xx = -mean_rv:dx:5*mean_rv;
factor = 1/sqrt(2*pi*smooth);

pdf = zeros(size(xx));
for jj = 1:Nrvs
    in  = (abs(xx-rvs(jj)) < mean_rv);
    ker = factor*in.*exp(-(xx-rvs(jj)).^2 / (2*smooth));
    if(sum(ker)~=0)
        pdf = pdf + ker/sum(ker)/Nrvs/dx;
    end
end

plot(xx,pdf,'r.');
axis([-mean_rv 4*mean_rv 0 1.1*max(theory)])

title(['\bf kernel-based pdf'])
xlabel('\bf x-axis')
ylabel(['\bf computed f(x) using ' num2str(Nrvs) ' points'])