clear all
close all
N = 10000;
n = 20;
alpha = 3;
lambda = 0.5;
alphaN = alpha * ones(1,N);
v = randg(alphaN)./lambda;
minv = min(v);
maxv = max(v);
level = linspace(minv,maxv,n+1);
for i = 1:n
  f(i) = length(find(v > level(i) & v <= level(i+1)));
end
f = n / (maxv - minv) * f / N;
hold on
bar((level(1:n)+level(2:n+1))/2,f)
f = @(x) lambda^alpha*x.^(alpha-1).*exp(-lambda*x)/gamma(alpha);
x = linspace(minv,maxv,100);
plot(x,f(x),'r')
