mercoledì 5 agosto 2009

I test di verifica di normalità

Numerose volte nei posts precedenti abbiamo proceduto con l'analisi assumendo che i dati seguissero una distribuzione normale (applicando quindi la statistica parametrica). Vediamo adesso come verificare che i propri dati seguano una distribuzione gaussiana. I test a disposizione sono molto numerosi, e la maggior parte di questi sono disponibili nei vari pacchetti di R. Vediamo qui di seguito l'applicazione di alcuni di essi, ed in particolare:


  • D'agostino's K-squared test

  • Jarque-Bera test

  • Anderson-Darling test for normality

  • Cramer-von Mises normality test

  • Lilliefors (Kolmogorov-Smirnov) test for normality

  • Shapiro-Francia normality test

  • Pearson chi-square normality test

  • Shapiro-Wilk normality test

  • Kolmogorov-Smirnov test



Per tutti gli esempi che analizzeremo, userò un set di dati ottenuto applicando la funzione rnorm(), la quale ci fornisce 50 numeri che seguono una distribuzione normale, nota la media (10) e la deviazione standard (2). Ci aspettiamo quindi che tutti i test che vedremo affermino che i dati seguono una distribuzione normale.

x <- rnorm(50, 10, 2)




D'agostino's K-squared test (D'Agostino-Pearson normality test)

Che io sappia questo test non è stato implementato in nessun pacchetto di R. Tuttavia, conoscendo le formule (disponibili su WikiPedia), è facilmente ricostruibile una funzione. Qui di seguito riporto il codice:


dagostino.pearson.test <- function(x) {
# from Zar (1999), implemented by Doug Scofield, scofield at bio.indiana.edu
DNAME <- deparse(substitute(x))
n <- length(x)
x2 <- x * x
x3 <- x * x2
x4 <- x * x3
# compute Z_g1
k3 <- ((n*sum(x3)) - (3*sum(x)*sum(x2)) + (2*(sum(x)^3)/n)) /
((n-1)*(n-2))
g1 <- k3 / sqrt(var(x)^3)
sqrtb1 <- ((n - 2)*g1) / sqrt(n*(n - 1))
A <- sqrtb1 * sqrt(((n + 1)*(n + 3)) / (6*(n - 2)))
B <- (3*(n*n + 27*n - 70)*(n+1)*(n+3)) / ((n-2)*(n+5)*(n+7)*(n+9))
C <- sqrt(2*(B - 1)) - 1
D <- sqrt(C)
E <- 1 / sqrt(log(D))
F <- A / sqrt(2/(C - 1))
Zg1 <- E * log(F + sqrt(F*F + 1))
# compute Z_g2
G <- (24*n*(n-2)*(n-3)) / (((n+1)^2)*(n+3)*(n+5))
k4 <- (((n*n*n + n*n)*sum(x4)) - (4*(n*n + n)*sum(x3)*sum(x)) -
(3*(n*n - n)*sum(x2)^2) + (12*n*sum(x2)*sum(x)^2) -
(6*sum(x)^4)) /(n*(n-1)*(n-2)*(n-3))
g2 <- k4 / var(x)^2
H <- ((n-2)*(n-3)*abs(g2)) / ((n+1)*(n-1)*sqrt(G))
J <- ((6*(n*n - 5*n + 2)) / ((n+7)*(n+9))) * sqrt((6*(n+3)*(n+5)) /(n*(n-2)*(n-3)))
K <- 6 + (8/J)*(2/J + sqrt(1 + 4/(J*J)))
L <- (1 - 2/K) / (1 + H*sqrt(2/(K-4)))
Zg2 <- (1 - 2/(9*K) - (L^(1/3))) / (sqrt(2/(9*K)))
K2 <- Zg1*Zg1 + Zg2*Zg2
pk2 <- pchisq(K2, 2, lower.tail=FALSE)
RVAL <- list(statistic = c(K2 = K2), p.value = pk2, method =
"D'Agostino-Pearson normality test\n\nK2 is distributed as Chi-squared
with df=2", alternative = "distribution is not normal", data.name =
DNAME)
class(RVAL) <- "htest"
return(RVAL)
}


L'ipotesi nulla è che i nostri dati seguono una distribuzione normale. Applichiamo la funzione ai nostri dati:


dagostino.pearson.test(x)

D'Agostino-Pearson normality test

K2 is distributed as Chi-squared with df=2

data: x
K2 = 1.2788, p-value = 0.5276
alternative hypothesis: distribution is not normal


Essendo p-value > 0.05, accettiamo l'ipotesi nulla: siamo di fronte a una distribuzione normale.




Jarque-Bera normality test

E' disponibile nel pacchetto tseries, che deve essere scaricato, ed installato.


library(tseries)
jarque.bera.test(x)

Jarque Bera Test

data: x
X-squared = 1.1181, df = 2, p-value = 0.5717


Anche per questo test, l'ipotesi nulla è la normalità del vettore x. Essendo p-value > 0.05, accettiamo l'ipotesi nulla.




Anderson-Darling test for normality

Test disponibile nella libreria nortest (libreria che raccoglie numerosi test di normalità).


library(nortest)
ad.test(x)

Anderson-Darling normality test

data: x
A = 0.1931, p-value = 0.8898


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0





Cramer-von Mises normality test


library(nortest)
cvm.test(x)

Cramer-von Mises normality test

data: x
W = 0.0289, p-value = 0.8576


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0




Lilliefors (Kolmogorov-Smirnov) test for normality


library(nortest)
lillie.test(x)

Lilliefors (Kolmogorov-Smirnov) normality test

data: x
D = 0.0808, p-value = 0.5715


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0




Shapiro-Francia normality test


library(nortest)
sf.test(x)

Shapiro-Francia normality test

data: x
W = 0.9918, p-value = 0.9454


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0




Pearson chi-square test for normality


library(nortest)
pearson.test(x)

Pearson chi-square normality test

data: x
P = 1.2, p-value = 0.991


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0




Shapiro-Wilk normality test

Questo è uno dei test più utilizzati, ed è già disponibile nel pacchetto di base di R (stats), quindi non necessita di alcuna installazione.


shapiro.test(x)

Shapiro-Wilk normality test

data: x
W = 0.9865, p-value = 0.834


Ipotesi nulla: i dati seguono una distribuzione normale.
p-value > 0.05: accetto H0




Kolmogorov–Smirnov test

Anche questo test è disponibile nel pacchetto di base. Può essere utilizzato per effettuare due diverse analisi. Possiamo ad esempio chiederci "i dati X e Y provengono dalla stessa distribuzione?". In questo caso la verifica si effettua in questo modo:


x <- rnorm(50)
y <- rnorm(30)

ks.test(x, y)

Two-sample Kolmogorov-Smirnov test

data: x and y
D = 0.1533, p-value = 0.7195
alternative hypothesis: two-sided


Essendo p-value > 0.05, accettiamo l'ipotesi nulla: i due vettori provengono dalla stessa distribuzione.

Un'altra verifica che possiamo effettuare è la seguente: "il vettore X proviene da una distribuzione normale, o da una distribuzione chi-quadro a 1 grado di libertà?". In R si procede così:


ks.test(x, "pnorm")

One-sample Kolmogorov-Smirnov test

data: x
D = 0.0931, p-value = 0.7436
alternative hypothesis: two-sided

ks.test(x, "pchisq", 1)

One-sample Kolmogorov-Smirnov test

data: x
D = 0.54, p-value = 3.564e-14
alternative hypothesis: two-sided


In questo caso, oltre a specificare il vettore da analizzare, si deve pure specificare la funzione di confronto, messa tra virgolette. Per testare altre distribuzioni, potete trovare la sintassi precisa digitando ??distribution: così leggete l'elenco delle distribuzioni implementate in R; quindi, richiamando l'help relativo alla distribuzione di vostro interesse (ad esempio help(Normal)), potete vedere la corretta sintassi da utilizzare (generalmente è il nome della funzione, preceduto dalla lettera p: pnorm, pgamma, pchisq).
Come potevamo prevedere, il test di normalità dà p-value > 0.05, mentre il confronto con la distribuzione chi-quadro dà p-value < 0.05. Quindi il vettore x segue una distribuzione normale.

6 commenti:

  1. Nel ringraziare per l'utilissimo blog e per il notevole lavoro a questo correlato, segnalo come il D'Agostino-Pearson normality test sia disponibile nel package fBasics. Per "x <- rnorm(50, 10, 2)", la funzione "dagoTest" di fBasics restituisce il seguente output:

    Title:
    D'Agostino Normality Test

    Test Results:
    STATISTIC:
    Chi2 | Omnibus: 2.341
    Z3 | Skewness: -1.504
    Z4 | Kurtosis: 0.281
    P VALUE:
    Omnibus Test: 0.3102
    Skewness Test: 0.1326
    Kurtosis Test: 0.7787

    RispondiElimina
  2. Sembra che il test di Pace risulti spesso più sensibile degli altri qui elencati-
    L. Pace,(1983): Un nuovo test per la verifica dell'ipotesi funzionale complessa, con particolare riferimento al controllo della normalità. Statistica Applicata, 16,299-308.

    RispondiElimina
  3. In internet circola un esempio di applicazione del test di K-S, quello sulla distribuzione del passaggio di una varietà di uccelli nelle diverse ore.
    La distribuzione osservata è:
    a=c(0,1,1,9,4).
    La distribuzione attesa potrebbe essere:
    b=c(3,3,3,3,3).
    Dopo aver trovate le distribuzioni cumulate si calcola, manualmente, D = 0,4667... (risultato che concorda con quanto presente su internet).
    Se invece si applica il test di Kolmogorov-Smirnov con R, si trova un valore di D diverso:
    > a=c(0,1,1,9,4)
    > b=c(3,3,3,3,3)
    > ks.test(a,b)

    Two-sample Kolmogorov-Smirnov test

    data: a and b
    D = 0.6, p-value = 0.3291
    alternative hypothesis: two-sided .......

    Se addirittura si usa il test di ks per confrontare la distribuzione osservata con una distribuzione uniforme (come penso sia lecito) la differenza è amaggiore.
    Non mi riesce di capire perchè ci sia questa differenza. C'è qualcuno che può aiutarmi???
    masgiorgi@gmail.com

    RispondiElimina
  4. domanda da principiante...

    come faccio ad estendere il comando del test di normalità a diverse distribuzioni di frequenza contemporaneamente? esempio (inventato di sana pianta): i miei dati riportano per anni diversi le distribuzioni di frequenza di una certa lunghezza (sempre la stessa, misurata con lo stesso sistema).

    rangelunghezza(m) 2000 2001 2002 2003
    0-1 1 0 6 3
    1-2 1 0 5 8
    2-3 2 0 4 3
    3-4 0 0 3 0
    4-5 0 0 2 0
    5-6 0 0 1 0

    come posso fare il test ks della normalità esteso a tutte gli anni in un solo comando ottenendo diversi D e p-values?

    RispondiElimina
  5. Well done! It is so well written and interactive. Keep writing such brilliant piece of work. Glad i came across this post. Last night even i saw similar wonderful R Programming tutorial on youtube so you can check that too for more detailed knowledge on R Programming.https://www.youtube.com/watch?v=rgFVq_Q6VF0

    RispondiElimina
  6. Very Impressive Data Science tutorial. The content seems to be pretty exhaustive and excellent and will definitely help in learning Data Science course. I'm also a learner taken up Data Science training and I think your content has cleared some concepts of mine. While browsing for Data Science tutorials on YouTube i found this fantastic video on Data Science. Do check it out if you are interested to know more.:-https://www.youtube.com/watch?v=gXb9ZKwx29U&t=237s

    RispondiElimina