% teorema del limite centrale
% sommiamo nv variabili indipendenti con distribuzione di fischer 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;
nu1=6;
nu2=10;
media = nu2/(nu2-2);
varianza = (2*nu2^2*(nu1+nu2-2))/(nu1*(nu2-2)^2*(nu2-4));
nbins=201;
ub=8; %sqrt(varianza) + media;
lb=0; %-sqrt(varianza) + media;
x=linspace(lb,ub,nbins);
dx=(ub-lb)/(nbins-1);
Fpdf = fpdf(x,nu1,nu2);
figure;
hold on;
plot(x,Fpdf,'r')
title('Fischer')

nv=1;
y= frnd(nu1,nu2,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= frnd(nu1,nu2,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= frnd(nu1,nu2,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= frnd(nu1,nu2,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= frnd(nu1,nu2,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= frnd(nu1,nu2,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= frnd(nu1,nu2,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')