clear all
close all
lambda = 1.5;
alpha = 1;
media = alpha/lambda;
varianza = alpha/lambda^2;
T = 100;
% numero di simulazioni
N = 5000;
counter = zeros(1,N);
for k = 1:N
  t = 0;
  t = t+randg(alpha)/lambda;
  while (t < T)
    counter(k) = counter(k)+1;
    t = t+randg(alpha)/lambda;
  end
  t = T;
end
E = sum(counter)/N;
abs(E/t-1/media)
countermax = max(counter);
mu = T/media;
sigma = sqrt(T*varianza/media^3);
nrange = [0:countermax];
for i = nrange
  f(i+1) = length(find(counter == i));
end
plot(nrange,f/N,'*',...
nrange,1/sqrt(2*pi*sigma^2)*exp(-((nrange-mu).^2)./(2*sigma^2)))