##
## ueb07 - Demonstration Histograms and Parzen Window Estimates
##
## Christoph Dalitz, 22.11.2011
## Hochschule Niederrhein
## Fachbereich Elektrotechnik und Informatik
## Krefeld, Germany
##

## histogram parameters
#n = 1000;       # number of random values
n = 100;       # number of random values
k = 10;        # number of equally spaced histogram bins
h = 0.5;       # window width
h2 = h/2;

## parameters of probability distribution
shape = 1.8;
scale = 7;
mu = 0;
sigma2 = 0.5;

## create random numbers and plot the points and the true PDF
##x = wblrnd(shape,scale,1,n);
x = normrnd(mu,sigma2,1,n);
x = sort(x);
xmin = x(1);
xmax = x(length(x));
hold off
plot(x,zeros(1,length(x)),"*;;");
hold on;
xx=[xmin-h2:0.01:xmax+h2];
##plot(xx,wblpdf(xx,shape,scale),"-2;true PDF;")
plot(xx,normpdf(xx,mu,sigma2),"-2;true PDF;")

## draw equally spaced histogram
## note that normalization must be computed from bin width:
## 1 = sum(binheights)*binwidth
binwidth = (xmax-xmin) / k;
#hist(x,k,1/binwidth);

## draw floating histogram
floathist = 0*xx;  # memory preallocation
i = 0;
for y = xx
  i += 1;
  left = (y-h2)*ones(1,n);
  right = (y+h2)*ones(1,n);
  floathist(i) = sum((left < x) & (x < right)) / (n*h);
endfor
plot(xx,floathist,"-4;floating histogram;");

## do gaussian kernel smoothing
h = h/4; h2 = h/2; # needs smaller window width
gkernel = 0*xx;
norm = 1/sqrt(2*pi);
i = 0;
for y = xx
  i += 1;
  gkernel(i) = norm * sum(exp(-0.5*((y-x)/h).^2)) / (n*h);
endfor
plot(xx,gkernel,"-5;gaussian smoothing;");