Est-ce bien Normale ?

mars 2017 statistiques

Est-ce que c'est normal ?

Comment étudier la distribution statistique de données expérimentales.

In [3]:
voir le code

Introduction

On va regarder ici les moyens pratiques d'étudier la distribution statistique d'un ensemble de mesures expérimentalles, et de tester son adéquation avec une loi théorique donnée, typiquement avec la loi Normale (Gaussienne).

Voici trois jeux de données différents, ayant chacun 100 points. Ces données peuvent par exemple être la largeur, la hauteur et le poids d'un échantillons de 100 pièces issus de la même ligne de production. Les variations observées peuvent alors venir du procédé de fabrication ou alors de la méthode de mesure.

In [4]:
voir le code

Il est habituel de caractériser de tels échantillons par la moyenne et l'écart-type, mais ce n'est pas ce qui nous interesse ici. Les valeurs ont été centrées et normées de façon à ce que les 3 ensembles aient la même moyenne (0) et le même écart-type (1). On distingue pourtant une différence entre ces trois jeux de données. On va voir si l'on peut caractériser plus finement cela.

Ce qui nous interesse ici est donc la distribution statistique des valeurs et comment l'étudier experimentalement.

Qu'est qu'une distribution statistique ?

La distribution statistique est pour une variable aléatoire la probabilité d'obtenir cette valeur: $f(x) = P[{\small x=x_0} ]$

Par exemple pour le lancé d'un dé la distribution est uniforme: $P[{\small x=1}]=P[{\small x=2}]=...=P[{\small x=6}]=1/6$

Les distributions peuvent être, entre-autre : - discretes ou (à support) continues - symétriques ou non-symétriques

Il existe un grand nombre de distributions statistiques, plus ou moins complexes: voir Gallery of Distributions, wikipedia ou encore la doc scipy.

La plus célèbre de toutes étant la distribution gaussienne, ou loi Normale: $\exp( - x^2 )$.

In [5]:
voir le code

Utilisation des histogramme ?

Nous souhaitons, à partir de données expérimentales, tracer la distribution statistique sous-jacente. Le moyen le plus connu, et intuitif, est de tracer un histogramme. Ce n'est en réalité pas le plus pratique et rigoureux (avec un ordinateur).

Avec un histogramme on calcul la probabilité, non d'obtenir une valeur unique de $x$ (ce qui n'est pas possible), mais la probabilité d'obtenir la valeur dans un intervalle: $$P_i\left[ b_i < x_m < b_{i+1} \right]$$

Il est alors necessaire de spécifier la taille de cet intervalle. La forme de la distribution obtenue est alors fortement sensible au nombre de valeur tombant dans chaque intervalle. Le graphique suivant illustre cela, avec un jeu de données de 100 points suivant une distribution normale.

In [6]:
voir le code

Pour qu'il ne soit pas necessaire de spécifier un intervalle d'observation, on utilise en pratique la fonction de répartition.

Fonction de répartition

La fonction de répartition $F(x)$ est la probabilité d'obtenir une valeur inférieure à $x$: $$F(x) = P[\, x_m < x \,]$$

$F(x)$ correspond à l'intégrale de la densité de probabilité $f(x)$. Il est très simple d'obtenir cette fonction numériquement. L'astuce est de trier les valeurs de $X$ par ordre croissant.

En effet, la probabilité empirique d'obtenir une mesure inférieure à $x_i$ est le nombre de valeur plus petite que celle-ci divisé par le nombre total de mesures ($N$). Une fois les valeurs triées, le nombre de valeur plus petites est simplement la position (l'indice) de cette valeur. On peut alors normer les indices par le nombre totale de valeurs ($N$).

$$ \begin{align} P[ x< min(X) ] &= 0 \\ P[ x< max(X) ] &= 1 \\ P[ x< x_i ] &= i/N \; avec \; X_i = sort(X)[i] \end{align} $$

ref. wikipedia

In [7]:
fig, ax = plt.subplots(1, 3, sharey=True, figsize=(n_graph*4, 4));
for i, X in enumerate( X_exemples.values() ):
    Xordered = np.sort( X )
    fraction = np.linspace(0, 1, len(X) )
    ax[i].plot( Xordered, fraction ,  '.k' );
    ax[i].set_xticks([-2, 0, 2]);

plt.subplots_adjust(wspace=0.07, hspace=.0);
ax[i].set_yticks([0, 1]);
ax[0].set_ylabel('$P[\,x_m<x\,]$');plt.xlabel('$x$');

Ce sont les mêmes données, mais representées d'une façon différentes. Il n'y a pas eu de calcul (autre qu'une normalisation) ou bien d'informations ajoutées. On constate qu'elles n'ont pas la même distribution.

Rq.: Pour obtenir la densité de probalité $f(x)$, il conviendrait de differencier ces courbes. Du fait du bruit ce n'est pas une bonne idée.

Diagrame Quantile-Quantile

Pour aller plus loin, il existe des moyen de comparer plus finement les distributions entre elles, ou avec des distributions théoriques. Par exemple, la première courbe est t'elle réellement différente de la troisième ?

Voir entre autre sur Wikipedia: Droite de Henry ou Diagramme Quantile-Quantile.

La fonction Scipy probplot permet de réaliser un diagrame Quantile-Quantile par rapport à une distribution de réference, et en même temps qu'une regression linéaire sur les points obtenus (doc).

In [8]:
import scipy.stats as stats
 
X = X_exemples['uniforme']    
stats.probplot(X, dist="norm", plot=plt)
plt.show()

Tests statistiques

Mais est-ce que l'écart observé est-il significatif ? à quoi comparer la différence observée ? Dans l'exemple ci-dessus on voit que le $R^2$ de la regression n'est pas forcement un bon indicateur de la justesse du fit.

Il existe plusieurs test statistique permettant de decider si, oui ou non, la distribution correspond à la distribution théorique spécifiée.

réferences:

Le test de Anderson-Darling est implémenté dans Scipy.stats pour les distributions Normale, exponentielle, logistique, et de Gumbel. (doc)

remarque: Les valeurs critiques sont fonction du nombre de points.

If A2 is larger than these critical values then for the corresponding significance level, the null hypothesis that the data come from the chosen distribution can be rejected. {‘norm’,’expon’,’logistic’,’gumbel’,’gumbel_l’, gumbel_r’, ‘extreme1’}

In [9]:
from scipy.stats import anderson
In [10]:
anderson(X_exemples['uniforme'] , dist='norm')
Out[10]:
AndersonResult(statistic=0.57113873819162109, critical_values=array([ 0.555,  0.632,  0.759,  0.885,  1.053]), significance_level=array([ 15. ,  10. ,   5. ,   2.5,   1. ]))
In [11]:
def testAnderson( X ):
    A2, Vc, Alphas = anderson( X  , dist='norm')

    k = 2
    if A2 > Vc[k]:
        print( "La distribution n'est PAS normale avec un risque de %f pourcent" % Alphas[k])
    else: 
        print( "La distribution EST normale avec un risque de %f pourcent" % Alphas[k])
        
#to do: chercher le risque minimal possible
In [12]:
for key, X in X_exemples.items():
    print( '%s: '%key )
    testAnderson( X )

uniforme: 
La distribution EST normale avec un risque de 5.000000 pourcent
exponential: 
La distribution n'est PAS normale avec un risque de 5.000000 pourcent
normal: 
La distribution EST normale avec un risque de 5.000000 pourcent
In [13]:
Xn = np.random.normal(loc=0.0, scale=1.0, size=200) 
anderson( Xn    , dist='norm')
Out[13]:
AndersonResult(statistic=0.3521176591327162, critical_values=array([ 0.565,  0.644,  0.772,  0.901,  1.071]), significance_level=array([ 15. ,  10. ,   5. ,   2.5,   1. ]))

Test du test sur un mélange de distribution

On peut étudier le comportement du test en l'évaluant sur un jeu de donnée construit comme mélange d'une distribution normale et d'une distribution uniforme. Le graphique suivant montre 5000 résultats de calculs du test d'Anderson sur des données de 400 points. La ligne horizontale rouge est la valeur critique pour un risque de 5%.

In [14]:
voir le code

Et avec des données réels ?

Par curiosité on peut regarder la distribution statistique de données réelles. On trouve l'historique du taux de change Dollar-Euro à l'adresse suivante : https://fred.stlouisfed.org/series/DEXUSEU

C'est une serie temporelle et les valeurs ont une certaine cohérence temporelle. On peut, en revanche, supposer que les variations à hautes fréquences correspondent à un bruit aléatoire. On regarde donc la premier dérivée du taux de change.

In [17]:
voir le code

Regardons la répartition des valeurs:

In [18]:
dYordered = np.sort( dY )
fraction = np.linspace( 0, 1, len(dY) )

plt.figure(figsize=(10, 4) )
plt.plot(  dYordered, fraction,  'k.', alpha=0.6, markersize=1);

plt.xlabel( 'valeur' ); plt.ylabel( ' $P[ x < X0 ]$' ); #$P( x < x^+ )$' );

Est-ce bien une distribution normale ?

In [19]:
testAnderson( dY )

La distribution n'est PAS normale avec un risque de 5.000000 pourcent

Ah ah! ... on peut comparer la distribution empirique avec celles théroriques:

In [23]:
from scipy.special import erf, erfinv

def Laplace(x):
    return .5 + .5*np.sign(x)*( 1 - np.exp( -np.abs(x)*np.sqrt(2) ) )

def Normale(x):
    r = .5*(1 + erf(x/np.sqrt(2) ))
    return r
In [36]:
dYordered_norme = (dYordered -dYordered.mean() )/dYordered.std()

plt.figure(figsize=(10, 4) )
plt.plot(  dYordered_norme, fraction,  'k.', alpha=0.4, markersize=1);

plt.plot(  dYordered_norme, Normale(dYordered_norme),  'r-', alpha=0.6, label='Normale');
plt.plot(  dYordered_norme, Laplace(dYordered_norme),  'b-', alpha=0.7, label='Laplace');

plt.xlabel( 'valeur' ); plt.ylabel( ' $P[ x < X0 ]$' ); #$P( x < x^+ )$' );
plt.xlim([-2, 2]); plt.legend();
In [ ]: