clear all
close all
% Y_i v.a. gamma
lambda = 1.5;
alpha = 1;
media = alpha/lambda;
varianza = alpha/lambda^2;
T = 100; % dovrebbe andare all'infinito
% numero di simulazioni
N = 5000;
counter = zeros(1,N);
for k = 1:N
  t = 0;
  t = t+randg(alpha)/lambda; % va bene qualunque altra v.a.
  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];
f = zeros(1,countermax+1);
for k = 1:N
  f(counter(k)+1) = f(counter(k)+1)+1;
end
plot(nrange,f/N,'*',...
nrange,1/sqrt(2*pi*sigma^2)*exp(-((nrange-mu).^2)./(2*sigma^2)))