Introduction aux tests de normalité statistiques avec Python

Le test de Shapiro-Wilk

Le test de Shapiro–Wilk, publié en 1965 par Samuel Shapiro et Martin Wilk, teste l'hypothèse nulle selon laquelle un échantillon $x_1, ..., x_n$ est issu d'une population normalement distribuée (voir article wikipedia). Ce test est directement disponible dans le module stats de scipy (voir Test de Shapiro-Wilk avec scipy) :

from scipy.stats import shapiro
shapiro(x, a=None, reta=False)

Paramètres d'entrée :

• x : matrice de données de dimension n.
• a (optionnel): matrice de dimension $n/2$
  Tableau de paramètres internes utilisés dans le calcul. Si ceux-ci ne sont pas donnés, ils seront calculés en interne.
• reta (optionnel): booléen, valeur par défaut: False. Retourne les paramètres internes utilisés dans le calcul.

Paramètres de sortie :

• W : float (test statistique)
• p-value : float (p-value pour le test d'hypothèse)
• a (optionnel): matrice, si reta = True.

Le test d'Anderson

Le test d'Anderson-Darling (voir article wikipedia) est une variante du test de Kolmogorov-Smirnov pour l'hypothèse nulle où un échantillon est tiré d'une population qui suit une distribution particulière. Pour le test d'Anderson-Darling, les valeurs critiques dépendent de la distribution testée (normale, exponentielle, logistique, ou Gumbel (Extreme Value Type I) distributions). Le test d'Anderson est également implementé dans le module stats de scipy:
(voir Test d'Anderson avec scipy)

from scipy.stats import anderson
anderson(x, dist='norm')

Paramètres d'entrée :

• x : matrice de données
• dist : {‘norm’,’expon’,’logistic’,’gumbel’,’extreme1’}(optionnel)
  type de distribution testée. La valeur par défaut est ‘norm’ (Note: ‘extreme1’ est synonyme de ‘gumbel’.

Paramètres de sortie :

• A2 : float, Le test d'Anderson-Darling
• critical : list, Les valeurs critiques par rapport à la distribution testée.
• sig : list
  Les niveaux de signification pour les valeurs critiques correspondantes en pourcentages. La fonction retourne des valeurs-seuil pour un ensemble différent de niveaux de signification en fonction de la distribution testée.

Notes :

Les valeurs critiques fournies sont pour les niveaux de signification suivants (Si A2 est supérieure à ces valeurs critiques alors pour le niveau de signification correspondante, l'hypothèse nulle que les données proviennent de la distribution choisie peut être rejetée}:

Valeurs critiques | | | | | |
------------- | ------------- | ------------- | ------------- | ------------- | ------------- | -------------
normal/exponenential | 15% | 10% | 5% | 2.5% | 1% |
logistic | 25% | 10% | 5% | 2.5% | 1% | 0.5%
Gumbel | 25% | 10% | 5% | 2.5% | 1% |

Les test D’Agostino (normaltest)

Cette fonction teste l'hypothèse nulle que l'échantillon provient d'une distribution normale. Il est basé sur D'Agostino et [R231], [R232] Le test de Pearson qui combine l'inclinaison et l'aplatissement de produire un test d'omnibus de la normalité.
Page Scipy sur le Test D’Agostino (normaltest)

scipy.stats.mstats.normaltest(a, axis=0)

Paramètres d'entrée :

• a : matrice de données
  axis : int or None
  Si Non, la matrice est traitée comme un ensemble unique de données, quel que soit sa forme. Sinon, dimension du tableau est testée.

Paramètres de sortie :

• k2 : float or array
  s^2 + k^2, où s est le z-score retourné par skewtest et k est le z-score retourné par kurtosistest.
• p-value : float or array
  Un chi2 faces carré probabilité pour le test d'hypothèses.

Exemple

Analyses graphiques

Transformation de Box-Cox

source: Cours pdf sur les Tests de normalité

Téléchargez le fichier de données: [attachment:198]

Téléchargez le script Python [attachment:199]

Execution du script: python NormalityTests.py

Source Code

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
from random import shuffle
from random import choice
from scipy.stats import norm
from scipy.stats import pearsonr
from scipy.stats import shapiro
from scipy.stats import kstest
from scipy.stats import anderson
import scipy.stats as stats

x = np.loadtxt("DataFileNormalityTest.txt", unpack=True)
x = np.sort(x)

print '#++++++++++++++++++++ Henry Plot (Plot 1) ++++++++++++++++++++#'

NbData = x.shape
QuantilesTheoriques = np.zeros(NbData[0])
for i in np.arange(NbData[0]):
    #print x[i],i+1, (i+1-0.375)/(50+0.25), norm.ppf( (i+1-0.375)/(50+0.25) )
    QuantilesTheoriques[i] = norm.ppf( (i+1-0.375)/(50+0.25) )

print 'mean', np.mean(x)
print 'std', np.std(x)
print 'std  ddof=1', np.std(x, ddof=1)
print 'shapiro', shapiro(x)
print 'kstest', kstest(x, 'norm')
print 'anderson', anderson(x)

fig = plt.figure()
plt.plot(x,QuantilesTheoriques,'o', markersize=2.5)
plt.xlabel(u'Observed quantile (X)')
plt.ylabel(u'Theoretical quantiles (U)')
plt.title("Henry Plot")
plt.text(1000, 2, "r = "+str( round(r_row,4) ))
plt.grid()
plt.savefig('HenryPlot.png', bbox_inches='tight')

print '#++++++++++++++++++++ Box Cox (Plot 2) ++++++++++++++++++++#'

lambda_tab = np.arange(-2,2.1,0.1)
lambda_tab_shape = lambda_tab.shape
pearsonr_tab = np.zeros( lambda_tab_shape[0] )
j = 0
for l in lambda_tab:
    x2 = np.zeros(NbData[0])
    if abs(l) < 0.001:
        for i in np.arange(NbData[0]):
            x2[i] = np.log(x[i])   
    else:
        for i in np.arange(NbData[0]):
            x2[i] = ( x[i]**l - 1 ) / l    
    r_row, p_value = pearsonr(x2,QuantilesTheoriques)
    pearsonr_tab[j] = r_row
    j = j + 1

print 'r max:',  pearsonr_tab.argmax(), pearsonr_tab.max(), pearsonr_tab[pearsonr_tab.argmax()]   
print 'l ', lambda_tab[pearsonr_tab.argmax()]

fig = plt.figure()
plt.plot(lambda_tab,pearsonr_tab,'o', markersize=2.5)
plt.plot(lambda_tab,pearsonr_tab,'b-', markersize=2.5)
plt.xlabel(r'Lambda $\lambda$')
plt.ylabel(u'r(theoretical quantiles : observed quantile)')
plt.title("Box-Cox Normality Plot")
plt.grid()
plt.savefig('BoxCox.png', bbox_inches='tight')

print '#++++++++++++++++++++ Henry Plot (Plot 3) ++++++++++++++++++++#'

l = lambda_tab[pearsonr_tab.argmax()]
if abs(l) < 0.001:
    for i in np.arange(NbData[0]):
        x2[i] = np.log(x[i])   
else:
    for i in np.arange(NbData[0]):
        x2[i] = ( x[i]**l - 1 ) / l

print 'mean', np.mean(x2)
print 'var', np.var(x2)
print 'std', np.std(x2)
print 'std  ddof=1', np.std(x2, ddof=1)
print 'Median', np.median(x2)
print 'shapiro', shapiro(x2)
print 'kstest', kstest(x2, 'norm',N=50)
print 'anderson', anderson(x2, dist='norm')
print 'normaltest', stats.normaltest(x2)

fig = plt.figure()
plt.plot(x2,QuantilesTheoriques,'o', markersize=2.5)
plt.xlabel(u'Observed quantile (X)')
plt.ylabel(u'Theoretical quantiles (U)')
plt.title("Henry Plot" + r' ( $\lambda$ = ' + str( lambda_tab[pearsonr_tab.argmax()] ) + ' )' )
plt.text(9, 2, "r = "+str( round(pearsonr_tab.max(),4) ))
plt.grid()
plt.savefig('HenryPlot2.png', bbox_inches='tight')

Référence

Liste non exhaustive des pages web consultées lors de la rédaction de cet article: