#-----------------------------------------------------# # # # LABORATORIO DI METODI STATISTICI # # PER LA FINANZA # # PROF F. BARTOLUCCI # # # # LAB-04 # #-----------------------------------------------------# ########################################################################## # Lezione 4 ##########################################################################1 # # Indice: # # * test di ipotesi sulla media di una popolazione # * test di ipotesi per il confronto tra le medie di due popolazioni # * test di ipotesi per la varianza di una popolazione # * test di ipotesi per il confronto tra le varianze di due popolazioni # * funzione di verosimiglianza # ########################################################################## # --------------------------------------------------- # TEST SU UNA SINGOLA MEDIA # --------------------------------------------------- # Caricare dati # Il file retribuzioni.txt contiene informazioni sulle retribuzioni # orarie di un insieme di 7986 individui dati = read.table("retribuzioni.txt", header = TRUE) head(dati) retribuzione = dati$retribuzione # verifica di hp sulla media # ----------------------------------- # 1) Sottoporre a verifica il seguente sistema di hp: H0: mu = 16, H1: mu > 16 # (assumendo varianza nota e pari a 75, alpha = 0.05) # Ampiezza campionaria n = length(retribuzione) mu0 = 16 sigma0 = sqrt(75) # Media campionaria xm = mean(retribuzione) # Costruiamo la statistica test z = (xm-mu0)/ (sigma0 *sqrt(1/n)) z # Calcoliamo il quantile della distribuzione Normale standard # (ipotesi alternativa unidirezionale) che lascia a destra una prob di 0.05 zal = qnorm(1 - 0.05) zal # Risultato del test # 1) Regione di rifiuto z > zal # Conclusione: non R H0 # 2) P-value pv = 1 - pnorm(z) pv < 0.05 # Conclusione: non R H0 # 2) Sottoporre a verifica il seguente sistema di hp: H0: mu = 16, H1: mu > 16 # (assumendo varianza non nota, alpha = 0.05) # Varianza campionaria s2 = var(retribuzione) # Calcoliamo la statistica test t = (xm-mu0)/sqrt(s2/n) # Calcoliamo il quantile della t di Student con n-1 gld # (ipotesi alternativa unidirezionale) che lascia a destra una prob di 0.05 gdl = n-1 tal = qt(1-0.05,gdl) # Risultato del test t > tal # 2) P-value pv = 1 - pt(t,gdl) pv < 0.05 # Conclusione: non R H0 -- non c'è sufficiente evidenza statistica sfavorevole ad H0 # ---------------------------------------------------- # Test su due medie # ---------------------------------------------------- # 1) Sottoporre a verifica il seguente sistema di ipotesi -- # H0: mu_laurea = mu_nonLaurea, H1: mu_laurea != mu_nonLaurea # (assumendo varianze non note ma uguali, alpha = 0.05) # Selezioniamo i sottocampioni retribuzione1 = dati$retribuzione[dati$laurea=="SI"] retribuzione2 = dati$retribuzione[dati$laurea=="NO"] # Calcoliamo le ampiezze campionarie n1 = length(retribuzione1) n2 = length(retribuzione2) # Calcoliamo le media campionarie xm = mean(retribuzione1) ym = mean(retribuzione2) # calcoliamo le varianze campionarie s21 = var(retribuzione1) s22 = var(retribuzione2) # Calcoliamo la varianza pooled s2p = ((n1-1)*s21 + (n2-1)*s22)/ (n1+n2-2) # Costruiamo la statistica test t = (xm-ym) / sqrt(s2p * (1/n1 + 1/n2)) # Calcoliamo il quantile della distribuzione t di student con n1+n2-2 gld (ipotesi alternativa bidirezionale) gdl = n1 + n2 -2 tal = qt(1-0.05/2, gdl) # Risultato del test t tal abs(t) > tal pv = 2*(1 - pt(abs(t), gdl)) pv pv < 0.05 # Conclusione: R H0 # le retribuzioni medie nei due gruppi sono significativamente diverse # Ripetiamo il test utilizzando funzione t.test di R t.test(retribuzione1, retribuzione2, mu = 0, alternative="two.sided", var.equal=TRUE) # verifica di hp sulla varianza # ----------------------------------- # 1) Verificare l'ipotesi che la varianza variabile retribuzione oraria sia # pari a 75 contro ipotesi sia maggiore (alpha = 0.05) # Calcoliamo l'ampiezza campionaria n = length(retribuzione) # valore del parametro sotto H0 sigma0 = 75 # Calcoliamo la varianza campionaria s2 = var(retribuzione) # Calcoliamo la statistica test v = (n-1) * s2 / sigma0 # sotto H0 v ha una distribuzione Chi-quadro con n-1 gdl # Calcoliamo il quantile val = qchisq(1 - 0.05, n-1) v > val # P-value pv = 1 - pchisq(v, n-1) pv < 0.05 # Conclusione: non R H0 # Test su due varianze # --------------------------------------------------- # 2) Verificare l'ipotesi che la varianza della retribuzione oraria # nel gruppo dei laureati e dei non laureati sia uguale, contro l'ipotesi alternativa # che sia diversa (alpha = 0.05) # Calcoliamo le varianze campionarie nei due gruppi s21 = var(retribuzione1) s22 = var(retribuzione2) # Cacloliamo la statistica test f = s21 / s22 # Calcoliamo i quantili della distribuzione F di Fisher con n1-1 ed n2-1 gdl # che lasciano a sinistra e a destra una probabilità di 0.05/2 f1al = qf(0.05/2, n1-1, n2-1) f2al = qf(1-(0.05/2), n1-1, n2-1) # Risultato del test f < f1al | f > f2al # Conclusione: R H0 # Ripetiamo il test utilizzando funzione var.test di R help(var.test) var.test(retribuzione1, retribuzione2, alternative="two.sided") # funzione di verosimiglianza # ----------------------------------------- # variabile con distribuzione di Bernoulli # ----------------------------------------- # estraiamo un campione di 1000 osservazioni da una distribuzione di Bernoulli(0.3) y = rbinom(1000, 1, 0.5) # Rappresentiamo likelihood # griglia di valori per p pv = seq(0,1,0.001) # numero di punti considerato np = length(pv) Lv = rep(0,np) # calcoliamo il valore di L al variare di pv for(i in 1:np) Lv[i] = prod(dbinom(y,1,pv[i])) plot(pv,Lv,type="l") # Rappresentiamo log-likelihood lv = log(Lv) plot(pv,lv,type="l") # massimo vicino a 0.5? phat = sum(y)/1000 phat # variabile con distribuzione Normale # ----------------------------------------- # Estraiamo 1000 osservazioni da una Normale standard y = rnorm(100) # Rappresentiamo la likelihood rispetto a mu # griglia di valori per mu muv = seq(-1,1,0.01) # dimensione della griglia np = length(muv) # valore di L al variare di mu Lv = rep(0,np) for(i in 1:np) Lv[i] = prod(dnorm(y,muv[i])) # Rappresentiamo la likelihood plot(muv,Lv,type="l") # Rappresentiamo log-likelihood lv = log(Lv) plot(muv,lv,type="l") # Massimo vicino a 0? muhat = sum(y)/100 muhat # salviamo ed usciamo da R savehistory()