%--------------------------------------------------------------------------------
% Quinto  laboratorio 17 Dicembre 2002
% Argomento: intervalli di confidenza
%--------------------------------------------------------------------------------

%--------------------------------------------------------------------------------
% 1. INTERVALLO DI CONFIDENZA SIMMETRICO PER LA MEDIA DI UNA POPOLAZIONE CON VARIANZA NOTA
%
%            mu+- = X_n +- z_beta * sigma / sqrt(n)
%--------------------------------------------------------------------------------

% generiamo un campione di 10 osservazioni con legge normale standard
% e dimentichiamo di averlo generato ma supponiamo che sia frutto dell'osservazione
% di un determinato fenomeno che IPOTIZZIAMO abbia media <mu> = 0 e che sappiamo con CERTEZZA
% avere varianza <sigma> = 1

mu	 = 0;
n         	 = 10;
sigma	 = 1;
x        = normrnd (mu, sigma, n, 1);

% ed ora cerchiamo l'ampiezza di un intervallo SIMMETRICO rispetto a xmedio
% che contenga la media con probabilita' del 95% (chiariremo in seguito il significato di questa espressione)

alpha    	 = 0.95;
xmedio	 = mean (x);
beta	 = (1 + alpha)/2;
zbeta        = norminv (beta, 0, 1);
errore   	 = zbeta*sigma/sqrt(n)
intervallo 	 = [xmedio-errore xmedio+errore]

% proviamo con diversi valori di <n>: come varia l'ampiezza dell'intervallo
% al crescere di n?

for n = 1 : 1000
	x = normrnd (mu, sigma, n, 1);
	xmedio = mean (x);
	zbeta = norminv (beta, 0, 1);
	errore(n) = zbeta*sigma/sqrt(n);
end

plot(errore)


% NOTA: con questo procedimento non siamo arrivati a conoscere il valore vero della media
% (ovvero il valore atteso della v.a. che modellizza il fenomeno)
% ma ora sappiamo che, ripetendo un numero grande di volte (<nprove>) l'esperimento,
% nel 95% dei casi il valore vero della media sarà contenuto nell'intervallo centrato in <xmedio> e di ampiezza
% pari a <errore>

% generiamo 200 intervalli di confidenza al 95%. (Ancora una volta dobbiamo "fingere" che i dati così generati
% siano invece frutto dell'osservazione ripetuta di un fenomeno)
% Quanti di questi contengono il valore vero della media?

n      = 50;
nprove = 200;
x      = normrnd (mu, sigma, n, nprove);
xmedio = mean (x);
zbeta       = norminv (beta, 0, 1);
mumeno = xmedio-zbeta*sigma/sqrt (n);
mupiu  = xmedio+zbeta*sigma/sqrt (n);
nok    = sum (mumeno<mu & mu<mupiu);
frazione = nok/nprove % <frazione> dovrebbe essere vicino ad <alpha>

% vediamo la variabilita' degli intervalli:
plot(mumeno);
hold on
plot(mupiu);
hold off


%------------------------------------------------------------------------------------
% 2. INTERVALLO DI CONFIDENZA SIMMETRICO PER LA FREQUENZA DI POPOLAZIONE BERNOULLIANA
%
%             La formula è ancora la stessa ma questa volta approssimiamo sigma:
%
%                mu+- = X_n +- z_beta * sqrt{X_n * (1 - X_n) / n)
%
%                !!!!!! assicurarsi che np>5 e n(1-p)>5 !!!!!!
%------------------------------------------------------------------------------------

p      = 0.2;
mu     = p;
alpha  = 0.95;
beta	 = (1 + alpha)/2;
n      = 100;
nprove = 200;
x      = binornd (1, p, n, nprove);
xmedio = mean (x);
sigma  = sqrt (xmedio.*(1-xmedio)); % attenzione al punto!!!
zbeta  = norminv (beta, 0, 1);
mumeno = xmedio-zbeta*sigma/sqrt (n);
mupiu  = xmedio+zbeta*sigma/sqrt (n);
nok    = sum (mumeno<mu & mu<mupiu);
frazione = nok/nprove
% la variabilita' di mumeno e mupiu e' grande:
min(mumeno), max(mumeno), min(mupiu), max(mupiu)


%----------------------------------------------------------------------------------------------
% 3. INTERVALLO DI CONFIDENZA SIMMETRICO PER LA MEDIA DI UNA POPOLAZIONE CON VARIANZA INCOGNITA
%
%                mu+- = X_n +- t_beta(n-1) S_n / sqrt(n)
%----------------------------------------------------------------------------------------------

% La formula è formalmente identica alla precedente a patto di sostituire la varianza con la
% varianza campionaria e la statistica di riferimento (t di Student al posto della Normale)

% Il procedimento è il medesimo di prima

mu       = 0;
n        = 10;
sigma    = 1;
x        = normrnd (mu, sigma, n, 1);

alpha      	 = 0.95;
beta	 = (1+alpha)/2;
xmedio    	 = mean (x);
S       	 = std(x);
tbeta   	 = tinv (beta, n-1);
errore   	 = tbeta*S/sqrt(n)
intervallo	 = [xmedio-errore xmedio+errore]

% proviamo con diversi valori di <n>: come varia l'ampiezza dell'intervallo
% al crescere di n?

for n = 1 : 1000
	x = normrnd (mu, sigma, n, 1);
	xmedio = mean (x);
	tbeta = tinv (beta, n-1);
	errore(n) = tbeta*S/sqrt(n);
end

plot(errore)

% confrontiamo gli intervalli di confidenza per la media
% con varianza nota e quelli con varianza incognita: 
% si puo' dire a priori quali sono piu' ampi? [no]
zbeta = norminv (beta, 0, 1);
errore2   = zbeta*sigma/sqrt(n)
intervallo2 = [xmedio-errore2 xmedio+errore2]
bar([errore errore2])

% confrontiamo ora gli intervalli di confidenza per la media
% con varianza nota, usando come varianza la varianza campionaria calcolata
% invece di sigma^2, e quelli con varianza incognita: 
% si puo' dire a priori quali sono piu' ampi? [si, errore>errore3]
zbeta = norminv (beta, 0, 1);
errore3   = sbeta*S/sqrt(n)
intervallo3 = [xmedio-errore3 xmedio+errore3]
bar([errore errore3])


% generiamo 200 intervalli di confidenza al 95%.
% Quanti di questi contengono il valore vero della media?

n      = 50;
nprove = 200;
x      = normrnd (mu, sigma, n, nprove);
xmedio = mean (x);
tbeta  = tinv (beta,n-1);
S      = std(x);
mumeno = xmedio-tbeta*S/sqrt (n);
mupiu  = xmedio+tbeta*S/sqrt (n);
nok    = sum (mumeno<mu & mu<mupiu);
frazione = nok/nprove

% vediamo la variabilita' degli intervalli:
plot(mumeno);
hold on
plot(mupiu);
hold off


% confrontiamo il valore ottenuto con la t di student, il
% valore ottenuto approssimando la t di Student con la gaussiana,
% il valore ottenuto considerando la varianza nota.
n=10;
i=1;
while n<1e3,
   xarr(i)    = n;
   x          = normrnd (0, sigma, n, 1);
   xmedio     = mean (x);
   S          = std(x);
   tbeta      = tinv (beta, n-1);
   errore(i)  = tbeta*S/sqrt(n);
   zbeta      = norminv (beta, 0, 1);
   erroren(i) = zbeta*S/sqrt(n);
   errore2(i) = zbeta*sigma/sqrt(n);
   i          = i+1;
   n          = n*2;
end
loglog(xarr, errore, xarr, erroren, xarr, errore2); % gli errori decrescono come n^-1/2
loglog(xarr, abs(erroren-errore)); 


%--------------------------------------------------------------------------------
% 4. INTERVALLO DI CONFIDENZA PER LA VARIANZA DI UNA POPOLAZIONE CON MEDIA NOTA
%
%             sigma+- = n Tn^2/chi2_{(1 pm alpha)/2}(n)
%         formula per l'intervallo unilatero della varianza:
%                sigma^2 < n Tn^2/chi2_{1-alpha}(n)
%--------------------------------------------------------------------------------



mu       = 0;
n        = 10;
sigma    = 1;
alpha    = 0.95;
beta1		= (1+alpha)/2;
beta2		= (1-alpha)/2;
beta		= 1-alpha;
x        = normrnd (mu, sigma, n, 1);
Tn       = sqrt (sum ((x-mu).^2)/n);
qcp      = chi2inv(beta1, n);
qcm      = chi2inv(beta2, n);
intervallo = [n*Tn^2/qcp n*Tn^2/qcm]
ampiezza = n*Tn^2/qcm-n*Tn^2/qcp
qc       = chi2inv(beta, n);
varmax   = n*Tn^2/qc

% provare con n=100, n=1000: come diminuisce l'ampiezza dell'intervallo
% al variare di n?

% generate 100 intervalli di confidenza al 95%. Quanti di questi
% contengono il valore vero di sigma^2?
n      = 50;
nprove = 200;
x      = normrnd (0, sigma, n, nprove);
xmedio = mean (x);
Tn     = sqrt (sum ((x-mu).^2)/n);
qcp    = chi2inv((1+alpha)/2, n);
qcm    = chi2inv((1-alpha)/2, n);
varmeno= n*Tn.^2/qcp;
varpiu = n*Tn.^2/qcm;
qc     = chi2inv(1-alpha, n);
varmax = n*Tn.^2/qc;
nok1   = sum (varmeno<sigma^2 & sigma^2<varpiu);
frazione1 = nok1/nprove % frazione dovrebbe essere vicino ad alpha
nok2   = sum (sigma^2 < varmax);
frazione2 = nok2/nprove % frazione dovrebbe essere vicino ad alpha
plot(varmeno)
hold on
plot(varpiu,'g')
plot(varmax,'r');
hold off



%--------------------------------------------------------------------------------------
% 5. INTERVALLO DI CONFIDENZA PER LA VARIANZA DI UNA POPOLAZIONE CON MEDIA INCOGNITA
%
%                sigma+- = (n-1) Sn^2/chi2_{(1 pm alpha)/2}(n-1)
%               formula per l'intervallo unilatero della varianza:
%                 sigma^2 < (n-1) Sn^2/chi2_{1-alpha}(n-1)
%--------------------------------------------------------------------------------------



mu       = 0;
n        = 10;
sigma    = 1;
alpha    = 0.95;
beta1		= (1+alpha)/2;
beta2		= (1-alpha)/2;
beta		= 1-alpha;
x        = normrnd (mu, sigma, n, 1);
Sn       = std(x);
qcp      = chi2inv(beta1, n-1);
qcm      = chi2inv(beta2, n-1);
intervallo = [n*S^2/qcp n*S^2/qcm]
ampiezza = n*S^2/qcm-n*S^2/qcp
qc       = chi2inv(beta, n-1);
varmax   = n*S^2/qc

% provare con n=100, n=1000: come diminuisce l'ampiezza dell'intervallo
% al variare di n?

% confrontiamo gli intervalli di confidenza per la varianza
% con media nota e quelli con media incognita.

Tn       = sqrt (sum ((x-mu).^2)/n);
qcp      = chi2inv(beta1, n);
qcm      = chi2inv(beta2, n);
intervallo2 = [n*Tn^2/qcp n*Tn^2/qcm]
ampiezza2 = n*Tn^2/qcm-n*Tn^2/qcp
qc       = chi2inv(beta, n);
varmax2   = n*Tn^2/qc
%----------------------------------------------------------------------------
%----------------------------------------------------------------------------
%----------------------------------------------------------------------------