%% Testo
% 1) Mostrare la funzione di densità della v.c $\chi^2_g$  con g che varia
% da 1 a 20 nell'intervallo x [0 30]. Inserire il grafico per ciascun g in
% un pannello separato. Inserire nel titolo di ogni grafico l'indicazione
% dei rispettivi gradi di libertà. Impostare il minimo e il massimo
% dell'asse x pari a [0 30]
% 
% 2) Calcolare in una v.c. $\chi^2_g$ con g=12 
% 
% a) la probabilità di ottenere un valore inferiore a 11.3
% 
% b) la probabilità di ottenere un valore compreso nell'intervallo [10 11]
% 
% c) la probabilità di ottenere un valore superiore a 20
% 
% d) il quantile 0.9 (ossia il valore che lascia alla sua destra una probabilità 
% del 10%)
% 
% 3) Estrarre 1000000 di numeri casuali dalla $\chi^2_g$ con g=12. 
% 
% Calcolare le frequenze relative empiriche 
% 
% a) di ottenere un valore inferiore a 11.3
% 
% b) di ottenere un valore compreso nell'intervallo [10 11]
% 
% c) di ottenere un valore superiore a 20
% 
% d) calcolare il quantile 0.9 (ossia il valore che lascia alla sua destra 
% 10000 casi)
% 
% Confrontare i risultati empirici con quelli teorici ottenuti al punto 2 

%% Soluzione

%% Studio della densità della v.c. Chi2 al variare dei gradi di libertà
% 1) Mostrare la funzione di densità della v.c $\chi^2_g$  con g che varia
% da 1 a 20 nell'intervallo x [0 30]. Inserire il grafico per ciascun g in
% un pannello separato. Inserire nel titolo di ogni grafico l'indicazione
% dei rispettivi gradi di libertà. Impostare il minimo e il massimo
% dell'asse x pari a [0 30]

x=0:0.1:30;
for j=1:20
    
    subplot(4,5,j)

    denschi2=chi2pdf(x,j);
    
    plot(x,denschi2)
    title(['$\chi^2_{' num2str(j) '}$'],'Interpreter','latex')
    xlim([0 30])
end

%% Funzione di densità ripartizione e quantili di una \chi^2
% 2) Calcolare in una v.c. $\chi^2_g$ con g=12
% a) la probabilità di ottenere un valore inferiore a 11.3
% b) la probabilità di ottenere un valore compreso nell'intervallo [10 11]
% c) la probabilità di ottenere un valore superiore a 20
% d) il quantile 0.9 (ossia il valore che lascia alla sua destra una probabilità del 10%)
% 

g=12;
% Prob X^2<b
b=11.3;
probMinore11e3=chi2cdf(b,g);
disp(['Prob ( \chi^2_g(12) <' num2str(b) ' = ' num2str(probMinore11e3) ')'])
% Prob (a< X^2<b)
a=10;
b=11;
disp(['(' num2str(a)  '< Prob \chi^2_g(12) <' num2str(b) ') = ' num2str(chi2cdf(b,g)-chi2cdf(a,g))])

% Prob X^2>b
b=20;
disp(['Prob ( \chi^2_g(12) >' num2str(b) ' = ' num2str(1-chi2cdf(b,g)) ')'])

% quantile 0.9
disp(['Quantile 0.9 in una \chi^2_g(12) =' num2str(chi2inv(0.9,g))])

%% Estrazione numeri casuali da una $\chi^2_g$
% 3) Estrarre 1000000 di numeri casuali dalla $\chi^2_g$ con g=12. 
% Calcolare le frequenze relative empiriche 
% a) di ottenere un valore inferiore a 11.3
% b) di ottenere un valore compreso nell'intervallo [10 11]
% c) di ottenre un valore superiore a 20
% d) calcolare il quantile 0.9 (ossia il valore che lascia alla sua destra 100000 casi)
% con le probabilità ottenute al punto sopra

% Genero numeri casuali dalla distribuzione chi2 con
% Osservazione: i risultati che seguono non sono riproducibili a meno di
% fissare un seed dei numeri casuali

% Per riprodurre i risultati occorre decommentare la riga che segue
% Fisso uno starting point casuale nel generatore di numeri casuali
% rng(1000)

n=1000000;
y=chi2rnd(12,n,1);

% Frequenza relativa empirica di valori inferiori a 11.3
% MODO 1: estrarre i numeri minori di 11.3 e vedere la lunghezza del vettore
freqrelMinore11e3=length(y(y<11.3))/n;
disp(['Frequenza relativa empirica di valori inferiori a 11.3=' num2str(freqrelMinore11e3)])
% MODO 2: fare direttamente la somma dei valori TRUE nel vettore Booleano y<11.3
freqrelMinore11e3CHK=sum(y<11.3)/n;

% controllo che il modo 1 sia esattamente uguale al modo 2
disp('Controllo sui modi di estrarre le frequenze')
isequal(freqrelMinore11e3,freqrelMinore11e3CHK)


% Frequenza relativa empirica di valori compresi tra 10 e 11
freqrel=length(y(y>10 & y<11))/n;
disp(['Frequenza relativa empirica di valori compresi tra 10 e 11=' num2str(freqrel)])

% Frequenza relativa empirica di valori superiori a 20
freqrel=length(y(y>20))/n;
disp(['Frequenza relativa empirica di valori superiori a 20=' num2str(freqrel)])

% Calcolo del quantile empirico 0.9
quanrequested=prctile(y,90);
disp(['Quantile emprico 0.9 =' num2str(quanrequested)])