% teorema del limite centrale
% sommiamo nv variabili indipendenti con distribuzione di poisson ottenute
% con np prove e confrontiamo la frequenza relativa della distribuzione
% della loro somma con la previsione teorica, vale a dire con una gaussiana
% di media e varianza note dalla teoria

clear all;
close all;
clc;
rng('default')

np=10000;
lambda = 4.5;
media = lambda;
varianza = lambda;
nbins=201;
ub=20; %sqrt(varianza) + media;
lb=0; %-sqrt(varianza) + media;
x=linspace(lb,ub,nbins);
dx=(ub-lb)/(nbins-1);
ppdf = poisspdf(x,lambda);
figure;
hold on;
plot(x,ppdf,'r')
title('Poisson')

nv=1;
y= poissrnd(lambda,nv,np);
sy=sum(y,1)/nv;
freq=hist(sy,x)/(dx*np);
dist=1/sqrt(2*pi*varianza/nv)*exp(-(x-media).^2/(2*varianza/nv));
figure;
hold on;
bar(x,freq);
plot(x,dist,'r')
title('n=1')

nv=2;
y= poissrnd(lambda,nv,np);
sy=sum(y,1)/nv;
freq=hist(sy,x)/(dx*np);
dist=1/sqrt(2*pi*varianza/nv)*exp(-(x-media).^2/(2*varianza/nv));
figure;
hold on;
bar(x,freq);
plot(x,dist,'r')
title('n=2')

nv=5;
y= poissrnd(lambda,nv,np);
sy=sum(y,1)/nv;
freq=hist(sy,x)/(dx*np);
dist=1/sqrt(2*pi*varianza/nv)*exp(-(x-media).^2/(2*varianza/nv));
figure;
hold on;
bar(x,freq);
plot(x,dist,'r')
title('n=5')

nv=30;
y= poissrnd(lambda,nv,np);
sy=sum(y,1)/nv;
freq=hist(sy,x)/(dx*np);
dist=1/sqrt(2*pi*varianza/nv)*exp(-(x-media).^2/(2*varianza/nv));
figure;
hold on;
bar(x,freq);
plot(x,dist,'r')
title('n=30')

nv=50;
y= poissrnd(lambda,nv,np);
sy=sum(y,1)/nv;
freq=hist(sy,x)/(dx*np);
dist=1/sqrt(2*pi*varianza/nv)*exp(-(x-media).^2/(2*varianza/nv));
figure;
hold on;
bar(x,freq);
plot(x,dist,'r')
title('n=50')

nv=100;
y= poissrnd(lambda,nv,np);
sy=sum(y,1)/nv;
freq=hist(sy,x)/(dx*np);
dist=1/sqrt(2*pi*varianza/nv)*exp(-(x-media).^2/(2*varianza/nv));
figure;
hold on;
bar(x,freq);
plot(x,dist,'r')
title('n=100')

nv=500;
y= poissrnd(lambda,nv,np);
sy=sum(y,1)/nv;
freq=hist(sy,x)/(dx*np);
dist=1/sqrt(2*pi*varianza/nv)*exp(-(x-media).^2/(2*varianza/nv));
figure;
hold on;
bar(x,freq);
plot(x,dist,'r')
title('n=500')