M4202C - Statistique inférentielle
TP 4 - Tests de conformité

2021/2022 - A. Ridard

Avant de commencer les tests

Importation des modules

In [1]:
import os # pour gérer les répertoires
import numpy as np # pour quelques fonctions mathématiques et pour stocker les données
import scipy.stats as sps # pour toutes les fonctions statistiques (moyenne, variance, tests, ...)
import matplotlib.pyplot as plt

Lecture des données dans un fichier .txt

In [2]:
os.chdir('D:\\Users\\ridard\\Desktop') # on change le répertoire de travail
f=open('data_1.txt','r') # on ouvre le fichier en mode lecture
lst=f.readlines()
x=''
for ligne in lst:
    x=x+ligne
f.close()
x
Out[2]:
'1695\t2202\t2722\t2534\t2494\t2648\t2298\t1997\t2118\t2767\n2867\t2391\t2029\t2121\t2105\t2565\t2652\t2497\t2822\t2713\n2014\t2350\t2343\t2398\t2505\t2630\t2169\t2661\t2325\t2031\n2683\t2328\t2710\t2417\t2264\t2299\t2531\t2423\t2592\t2577\n2568\t1992\t2872\t2603\t2415\t2072\t2475\t2089\t2140\t2720\n'

Stockage des données dans un tableau numpy

In [3]:
y=x.replace('\n','\t')
lst=y[:-1].split('\t')
lstc=lst[:]
for i in range(len(lst)):
    lstc[i]=eval(lst[i])
data=np.array(lstc)
data
Out[3]:
array([1695, 2202, 2722, 2534, 2494, 2648, 2298, 1997, 2118, 2767, 2867,
       2391, 2029, 2121, 2105, 2565, 2652, 2497, 2822, 2713, 2014, 2350,
       2343, 2398, 2505, 2630, 2169, 2661, 2325, 2031, 2683, 2328, 2710,
       2417, 2264, 2299, 2531, 2423, 2592, 2577, 2568, 1992, 2872, 2603,
       2415, 2072, 2475, 2089, 2140, 2720])

Description statistique des données

In [4]:
sps.describe(data)
Out[4]:
DescribeResult(nobs=50, minmax=(1695, 2872), mean=2408.6599999999999, variance=73987.208571428564, skewness=-0.3432569718756301, kurtosis=-0.5784175026497564)

NB : la variance fournie est la variance empirique corrigée

Test de normalité pour savoir si un échantillon peut être supposé gaussien

Première approche (graphique) : l'histogramme ressemble-t-il à une courbe de Gauss ?

In [14]:
stat = sps.describe(data)
n, moyEmp, varEmpC = stat.nobs, stat.mean, stat.variance
valmin, valmax = stat.minmax

x = np.linspace(valmin, valmax, 1000)
yTh = sps.norm(moyEmp, np.sqrt(varEmpC)).pdf(x)
plt.plot(x, yTh, 'r', label='Courbe de Gauss')

# Plusieurs possibilités pour déterminer le nombre de classes d'un histogramme

# Règle de Sturges
# nbCl = int(1 + np.log2(n))

# Règle de Rule (sans mauvais jeu de mots)
# nbCl = int(2*n**(1/3))

# Choix simple de la racine carrée
nbCl = int(np.sqrt(n))

plt.hist(data, bins=nbCl, normed=1, edgecolor='black', label='Histogramme')
plt.legend(loc='best')
plt.show()

Test de normalité à partir des coefficients d'asymétrie (skewness) et d'applatissement (kurtosis)

In [6]:
sps.normaltest(data)
Out[6]:
NormaltestResult(statistic=1.815472334971608, pvalue=0.40343650369431905)

La première composante de la réponse (statistic) est la réalisation $w$ de la variable de décision du test $W$
La deuxième (pvalue) est la probabilité sous $H_0$ (l'échantillon est gaussien) d'avoir $|W|>|w|$, c'est la p-value bilatérale

On peut alors rejeter $H_0$ dès que la p-value est inférieure à l'erreur de première espèce $\alpha$, mais attention il s'agit là d'un test bilatéral

NB : pour récupérer la p-value unilatérale, il suffit de diviser par deux la p-value fournie par Python !

Test de normalité de Shapiro-Wilk

In [7]:
sps.shapiro(data)
Out[7]:
(0.9690803289413452, 0.2122964859008789)

Test de conformité de la moyenne d'une gaussienne

Directement

In [8]:
sps.ttest_1samp(data,2400) # cf. exercice du cours
Out[8]:
Ttest_1sampResult(statistic=0.22512543915704525, pvalue=0.822817950327982)

A la main

D'abord, on calcule la réalisation $w$ de la variable de décision $W$

In [9]:
w=(moyEmp-2400)/np.sqrt(varEmpC/n)
w
Out[9]:
0.22512543915704525

Ensuite, on calcule la valeur critique c'est à dire l'extrémité de la zone de rejet (test unilatéral à droite avec une erreur de première espèce égale à 5%)

In [10]:
alpha=0.05
dl=n-1
valCrit=sps.t.ppf(1-alpha,dl)
valCrit
Out[10]:
1.6765508919142629

On décide

In [11]:
if w>valCrit:
    print('On rejette H0')
else :
    print('Aucune raison de rejeter H0')

Aucune raison de rejeter H0

On peut aussi retrouver la p-value bilatérale fournie par Python

In [12]:
pValBil=(1-sps.t.cdf(np.abs(w),dl))*2
pValBil
Out[12]:
0.822817950327982

Test de conformité de la variance d'une gaussienne

Ecrire un script "à la main" permettant de répondre à l'exercice du cours (toujours le même) concernant la variance.

In [ ]:
# Réponse