%------------------------------------------------------------------------------- %-------------------------------------------------------------------------------- % Sesto laboratorio 14 Gennaio 2003 % Argomento: test d'ipotesi %-------------------------------------------------------------------------------- % Test di ipotesi sulla media con varianza nota % Una ditta sta preparando un nuovo shampoo ed e' interessata a conoscere % l'altezza della schiuma che tale shampoo genera. Si suppone che tale altezza % sia normalmente distribuita con deviazione standard nota e pari a 20 mm. % Conoscendo i seguenti 16 dati riguardo all'altezza prodotta: % 194 195 226 224 222 178 163 180 174 205 156 235 186 166 158 185 % la compagnia vuole effettuare il seguente test di ipotesi sulla media: % H0 la media e' uguale a 180 mm. % H1 la media e' diversa da 180 mm. % 1) % Utilizzando soltanto i primi dieci dati e come statistica test la media % campionaria, trovare la regione di accettazione per l'ipotesi nulla ad un % livello di significativita' del 5% e trarre una conclusione sull'ipotesi fatta % 2) % Confrontare i risultati al punto 1) con quelli ottenuti con un livello di % significativita' pari all'1% % 3) % Confrontare i risultati al punto 1) con quelli ottenuti utilizzando tutti i % 16 dati a disposizione al livello di significativita` pari all'1%. % Calcolare, inoltre, il P-Value. Si poteva a priori concludere quale dei due intervalli % fosse piu` ampio? %------------------------------------------------------------------------------- % Impostiamo il problema: sigma = 20; Dati = [194 195 226 224 222 178 163 180 174 205 156 235 186 166 158 185]; mu = 180; % Ipotesi nulla %------------------------------------------------------------------------------- % 1) n = 10; % ampiezza del campione Xbar_n = mean(Dati(1:n)); % Poiche' alpha = P(rigettiamo l'ipotesi nulla anche se vera) cioe' % alpha = P(|Xbar - mu|*sqrt(n)/sigma > q_beta) quindi l'intervallo % e' dato da mu+- = Xbar_n +- q_beta * sigma / sqrt(n) alpha = 0.05; beta = 1 - alpha/2; q_beta = norminv(beta , 0 , 1); %q_beta = mynorminv(beta , 0 , 1); %se non funziona l'altra errore = q_beta * sigma / sqrt(n); mupiu = Xbar_n + errore; mumeno = Xbar_n - errore; intervallo = [mumeno mupiu] ampiezza = 2 * errore % Si noti che la regione di accettazione e` l'intervallo % di confidenza a livello 1-alpha % Verifichiamo se l'ipotesi nulla e' soddisfatta ovvero se mu appartiene all'intervallo: Test = mu > mumeno & mu < mupiu % Test = 1 ipotesi accettata, Test = 0 ipotesi rifiutata % Calcolo del P-Value = 2 * (1 - Phi(|Xbar_n - mu| * sqrt(n) / sigma)) PValue = 2 * (1 - normcdf(abs(Xbar_n - mu) * sqrt(n) / sigma)) %------------------------------------------------------------------------------- % 2) n = 10; % ampiezza del campione Xbar_n = mean(Dati(1:n)); alpha = 0.01; beta = 1 - alpha/2; q_beta = norminv(beta , 0 , 1); %q_beta = mynorminv(beta , 0 , 1); %se non funziona l'altra errore = q_beta * sigma / sqrt(n); mupiu = Xbar_n + errore; mumeno = Xbar_n - errore; intervallo = [mumeno mupiu] ampiezza = 2 * errore % Verifichiamo se l'ipotesi nulla e' soddisfatta ovvero se mu appartiene all'intervallo: Test = mu > mumeno & mu < mupiu % Test = 1 ipotesi accettata, Test = 0 ipotesi rifiutata % Calcolo del P-Value = 2 * (1 - Phi(|Xbar_n - mu| * sqrt(n) / sigma)) % e` inutile in quanto indipendente da alpha e quindi coincidente con quello % calcolato al punto 1 %------------------------------------------------------------------------------- % Vediamo come varia l'ampiezza dell'intervallo al diminuire di alpha. alpha(1) = 0.5 ; ampiezza(1) = 2 * norminv(1 - alpha(1)/2 , 0 , 1) * sigma / sqrt(n); %ampiezza(1) = 2 * mynorminv(1 - alpha(1)/2 , 0 , 1) * sigma / sqrt(n); for i = 2:10 ; alpha(i) = alpha(i-1) / 2 ; ampiezza(i) = 2 * norminv(1 - alpha(i)/2 , 0 , 1) * sigma / sqrt(n); %ampiezza(i) = 2 * mynorminv(1 - alpha(i)/2 , 0 , 1) * sigma / sqrt(n); end plot(alpha, ampiezza) %------------------------------------------------------------------------------- % 3) % Osserviamo che l'ampiezza dell'intervallo di confidenza a livello % beta e` amp(n,beta)=2*sigma^2*norminv((1+beta)/2)/sqrt(n) % e di conseguenza l'ampiezza dell'intervallo di accettazione del test % ad un livello di significativita` alpha (=1-beta) e` % delta(n,beta)=2*sigma^2*norminv(1-alpha)/2)/sqrt(n) % decrescente in entrambe le variabili. Pertanto % delta(10,0.05) mumeno & mu < mupiu % Test = 1 ipotesi accettata, Test = 0 ipotesi rifiutata % Si osservi che se varia n allora non si puo` piu` dire nulla a priori % sull'accettazione o meno dell'ipotesi nulla in quanto, anche se aumentasse % l'ampiezza dell'intervallo (se si diminuisse nel contempo alpha), % in generale cambiera` il punto medio dell'intervallo di accettazione %------------------------------------------------------------------------------- % Vediamo come varia l'ampiezza dell'intervallo al crescere di n. for i = 1:100 ; ampiezza(i) = 2 * norminv(1 - alpha/2 , 0 , 1) * sigma / sqrt(i); %ampiezza(i) = 2 * mynorminv(1 - alpha/2 , 0 , 1) * sigma / sqrt(i); %ampiezza(i) = 2 * MDnorminv(1 - alpha/2 , 0 , 1) * sigma / sqrt(i); end plot(1:100, ampiezza) %------------------------------------------------------------------------------- % Calcolo del P-Value = 2 * (1 - Phi(|Xbar_n - mu| * sqrt(n) / sigma)) PValue = 2 * (1 - normcdf(abs(Xbar_n - mu) * sqrt(n) / sigma)) %------------------------------------------------------------------------------- % Vediamo come varia il P-Value al crescere di n. % Generiamo 1000 dati casuali da una normale N(1,4) n=1000; Dati=normrnd(1,4,n,1); sigma=4; % Plottiamo l'andamento medio della media campionaria for i=1:n; Xbar(i)=mean(Dati(1:i)); end plot(1:n,Xbar) % Prendiamo ora come mu di riferimento del test mu=1 % Testiamo cosi` la bonta` del generatore mu=1; Xbar_n=0; for i = 1:n ; PValue(i) = 2 * (1 - normcdf(abs(Xbar(i) - mu) * sqrt(i) / sigma)); end plot(1:n, PValue) % Guardiamo il risultato su campioni di ampiezza multipla di 100 PValue(100:100:1000) % Prendiamo ora come mu di riferimento del test un'altro valore % (es. quello suggerito dal grafico di Xbar se differente da 1) mu=0.75; Xbar_n=0; for i = 1:n ; PValue(i) = 2 * (1 - normcdf(abs(Xbar(i) - mu) * sqrt(i) / sigma)); end plot(1:n, PValue) % Guardiamo il risultato su campioni di ampiezza multipla di 100 PValue(100:100:1000) %------------------------------------------------------------------------------- %------------------------------------------------------------------------------- % Test chiquadro di adattamento % L'esperimento di ibridazione di G. Mendel prevedeva, a livello teorico, % una distribuzione dei piselli nelle categorie RY RG WY WG pari a % 9/16 3/16 3/16 1/16. I risultati dell'esperimento furono % RY RG WY WG % 315 108 101 32 % La teoria di Mendel risulta confermata al 5%? % Qual e` il p-value del test? %------------------------------------------------------------------------------- % La statistica da utilizzare e` % Q= sum_i (F_i-FT_i)^2/FT_i % L'ipotesi nulla di adattamento si accetta a livello alpha % se e solo se Q < chi2_(1-alpha)(N_c-1-r) % dove N_c e` il numero di classi e r il numero di parametri della % distribuzione teorica stimati dai dati (in questo caso r=0, N_c=4) % Il p-value e` pvalue=1-F_chi2(N_c-1-r)(Q) % Creiamo il vettore delle frequenze assolute misurate F=[315, 108, 101, 32]; % Creiamo il vettore delle frequenze relative teoriche ft=[9/16,3/16,3/16,1/16]; % Creiamo il vettore delle frequenze assolute teoriche n=315+108+101+32; FT=ft*n; Q=sum((FT-F).^2./FT); ts=chi2inv(.95, 3); Test = Q < ts % Test = 1 ipotesi accettata, Test = 0 ipotesi rifiutata % Il p-value pvalue=1-chi2cdf(Q,3)