Statistique inférentielle

Screencast 1 : stats descs sur Essence

On importe le fichier contenant les consommations d’essence :

In [1]:
import pandas
essence=pandas.read_csv("essence.txt")

On calcule la moyenne, la variance et l’écart-type de l’échantillon :

In [2]:
import numpy, math
xbar=numpy.mean(essence['conso'])
round(xbar,2)
Out[2]:
31.45
In [3]:
sprime = numpy.std(essence['conso'], ddof=1)
round(sprime,2)
Out[3]:
2.16
In [4]:
n_essence = len(essence)
sprime2 = numpy.var(essence['conso'], ddof=1)
round(sprime2,2)
Out[4]:
4.66

Si on souhaite obtenir la variance “biaisée”, on peut l’obtenir ainsi :

In [5]:
v = numpy.var(essence['conso'])
round(v,2)
Out[5]:
4.63

On peut visualiser l’histogramme :

In [6]:
import matplotlib.pyplot as plt
%matplotlib inline
plt.hist(essence['conso'],bins=range(26,38))
plt.title("histogramme")
Out[6]:
<matplotlib.text.Text at 0x7fc0c6b52d68>

On peut visualiser sur cet histogramme la moyenne empirique :

In [7]:
plt.hist(essence['conso'],bins=range(26,38))
plt.title("histogramme")
plt.axvline(xbar, color='r', linewidth=2)
import matplotlib.lines as mlines
plt.legend(handles=[mlines.Line2D([], [], color='red',markersize=15, label="Moyenne empirique")])
Out[7]:
<matplotlib.legend.Legend at 0x7fc0c6a382e8>

On peut visualiser sur cet histogramme la moyenne empirique et le seuil “métier” :

In [8]:
mu0=31
plt.hist(essence['conso'],bins=range(26,38))
plt.title("histogramme")
plt.axvline(xbar, color='r', linewidth=2)
plt.axvline(mu0, color='b', linewidth=2)
import matplotlib.lines as mlines
plt.legend(handles=[mlines.Line2D([], [], color='red',markersize=15, label="Moyenne empirique"),mlines.Line2D([], [], color='blue',markersize=15, label=u"Seuil testé")])
Out[8]:
<matplotlib.legend.Legend at 0x7fc0c69222b0>

On peut enfin superposer sur cet histogramme la densité gaussienne :

In [9]:
import matplotlib.mlab as mlab
fig, ax = plt.subplots()
n, bins, patches = ax.hist(essence['conso'],normed=True,bins=range(26,38))
bins=numpy.arange(26,37,0.1)
y = mlab.normpdf(bins, xbar,sprime)*100
ax.plot(bins, y, '--')
ax.set_xlabel('Smarts')
ax.set_title(u"Histogramme et densité normale")
Out[9]:
<matplotlib.text.Text at 0x7fc0c690c470>

On constate que l’écart entre l’histogrammme, primo-estimation de la densité de probabilité, et la densité gaussienne sont plutôt proches.

Screencast 2 : IC sur Guerison

On importe le fichier contenant les guérisons ou non-guérisons :

In [10]:
guerison =pandas.read_csv("guerison.txt")

On peut estimer le taux de guerison théorique p :

In [11]:
n_guerison=len(guerison)
n_guerison_gueris= sum(guerison[guerison["guerison"]==1]["guerison"])
p_estim= n_guerison_gueris/float(n_guerison)

Si on souhaite encadrer le taux de guérison p avec une probabilité de 1α=95%, on obtient alors comme intervalle de confiance (Φ1α2=Φ0.9751.96) :

[1672161.96167216(1167216)216;167216+1.96167216(1167216)216].

Si on lance “manuellement” les calculs au niveau de test 5% :

In [12]:
alpha = 0.05
from scipy.stats import norm
icinf =p_estim-norm.ppf(1-alpha/2)*math.sqrt(p_estim*(1-p_estim)/n_guerison)
round(icinf,2)
Out[12]:
0.71999999999999997
In [13]:
icsup = p_estim+norm.ppf(1-alpha/2)*math.sqrt(p_estim*(1-p_estim)/n_guerison)
round(icsup,2)
Out[13]:
0.82999999999999996

On obtient alors : [0.72 ; 0.82]=[72% ; 83%] .

On constate que la largeur de l’intervalle n’est pas négligeable, mais n’oublions pas qu’il n’y a que 216 individus dans l’échantillon.

En pratique, le data analyst peut obtenir simplement cet intervalle à l’aide de la commande proportion_confint :

In [14]:
from statsmodels.stats.proportion import proportion_confint
proportion_confint(n_guerison_gueris,n_guerison, alpha=alpha, method='normal')
Out[14]:
(0.71729807581993277, 0.82899822047636351)

On aurait également pu obtenir un intervalle de confiance “exact”, basé sur la loi binomiale :

In [15]:
from statsmodels.stats.proportion import proportion_confint
proportion_confint(n_guerison_gueris,n_guerison, alpha=alpha, method='binom_test')
Out[15]:
(0.7110082209040385, 0.8245118012017411)

On constate que les résultats diffèrent très peu de l’intervalle de confiance asymptotique.

Enfin, si on avait choisi un niveau de confiance plus faible, 1α=90% par exemple, on aurait obtenu un intervalle de confiance plus étroit :

In [16]:
from statsmodels.stats.proportion import proportion_confint
proportion_confint(n_guerison_gueris,n_guerison, alpha=0.1, method='normal')
Out[16]:
(0.72627728992871754, 0.82001900636757874)

On obtient ici [0.726;0.820]=[72.6%;82.0%], la largeur de l’intervalle a bel et bien diminué !

Screencast 3 : IC sur moyenne Essence

Si on souhaite encadrer la consommation d’essence moyenne (théorique) μ avec une probabilité de 95%, on obtient alors comme intervalle de confiance tn1,1α2=t127,0.9751.97 et s2.16:

[31.451.972.16128;31.45+1.972.16128]

Remarquons ici que l’hypothèse gaussienne n’était pas obligatoire ici, en effet l’échantillon est de taille suffisamment importante (supérieure à 30).

Si on lance “manuellement” les calculs au niveau de test 5% :

In [17]:
import scipy
icinf=xbar-scipy.stats.t.ppf(1-alpha/2,n_essence-1)*sprime/math.sqrt(n_essence)
round(icinf,2)
Out[17]:
31.07
In [18]:
icsup=xbar+scipy.stats.t.ppf(1-alpha/2,n_essence-1)*sprime/math.sqrt(n_essence)
round(icsup,2)
Out[18]:
31.829999999999998

On obtient alors : [31.07 ; 31.83] .

Screencast 4 : IC sur variance Essence

Si on souhaite encadrer la variance (théorique) de la consommation d’essence σ2 avec une probabilité de 95%, on obtient alors comme intervalle de confiance (χ2n1,α2=χ2127,0.02597.7,χ2n1,α2=χ2127,0.975160.1ets24.66) : [127×4.66160.1 ; 127×4.6697.7] .

Si on lance “manuellement” les calculs au niveau de test 5% :

In [19]:
icinf = (n_essence-1)*sprime2/scipy.stats.chi2.ppf(1-alpha/2,n_essence-1)
round(icinf,2)
Out[19]:
3.7000000000000002
In [20]:
icsup = (n_essence-1)*sprime2/scipy.stats.chi2.ppf(alpha/2,n_essence-1)
round(icsup,2)
Out[20]:
6.0599999999999996

On obtient alors : [3.70 ; 6.06] .

Remarquons encore ici que l’hypothèse gaussienne n’était pas obligatoire ici, en effet l’échantillon est de taille suffisamment importante (supérieure à 30).

Screencast 5 : test sur Guerison

On teste : {H0:p=p0H1:p>p0

avec p0=0.75.

On considère que l’hypothèse gaussienne est acceptable ici, en effet : np0(1p0)=216(0.75)(10.75)=40.5>5 .

Pour α=5%, on a :

Φ1α1.64 p0+Φ1αp0(1p0)n0.798.

On ne rejette pas H0 au niveau de test 5% car : x¯=0.7730.798p0+Φ1αp0(1p0)n.

La p-valeur vaut : pvaleur=1Φ(nx¯p0p0(1p0))1Φ(0.786)0.216.

On constate que le niveau de test retenu devrait au moins être égal à 21,6% pour rejeter H0 ! Accepter de se tromper dans plus de 20% des cas lorsque rejette l’hypothèse nulle n’est pas rencontré en pratique.

Au final, le laboratoire ne peut pas conclure avec un niveau de test raisonnable que le nouveau médicament est meilleur que celui déjà sur le marché.

In [21]:
import statsmodels
statsmodels.stats.proportion.binom_test(n_guerison_gueris,n_guerison, prop=0.75, alternative='larger')
Out[21]:
0.24183479554062609

Screencast 7 : test sur moyenne Essence

On teste : {H0:μ=μ0H1:μμ0 avec μ0=0.31.

Pour α=5, on a :

tn1,1α2=t127,0.9751.97|nx¯μ0s||12831.45¯312.16|2.35

On rejette H0 au niveau de test 5% car : |nx¯μ0s|=2.35>1.97tn1,1α2.

La p-valeur vaut :

pvaleur=P(|τ(n1)|nx¯μ0s)=2(1P(τ(n1)nx¯μ0s))2(1P(τ(127)2.35))0.02.

On constate qu’on rejette bien l’hypothèse nulle au niveau de test 5% (0.02<0.05) mais pas au niveau de test 1% (0.02≥1%).

En pratique, le data analyst peut utiliser la commande ttest_1samp pour résoudre ce test :

In [22]:
scipy.stats.ttest_1samp(essence["conso"],31)
Out[22]:
Ttest_1sampResult(statistic=2.3543575167619162, pvalue=0.020088330843337596)

Screencast 8 : test sur variance Essence

Considérons par exemple le test suivant : {H0:σ2=σ02H1:σ2σ02

avec σ02=4.5.

Screencast 9 : test comparaison sur Iris

Considérons dans le (fameux) jeu de données historique des iris de Fisher la longueur des pétales d’iris des variétés et , sachant que nous disposons d’échantillons de 50 individus de chaque variété.

In [23]:
from sklearn import datasets
Iris = datasets.load_iris()
In [24]:
def function(x):
    return Iris['target_names'][x]
iris=pandas.DataFrame(Iris["data"])
iris.columns = Iris['feature_names']
iris["Species"]=Iris['target']
iris["Species"]=iris["Species"].apply(function)
In [25]:
iris
Out[25]:
sepal length (cm) sepal width (cm) petal length (cm) petal width (cm) Species
0 5.1 3.5 1.4 0.2 setosa
1 4.9 3.0 1.4 0.2 setosa
2 4.7 3.2 1.3 0.2 setosa
3 4.6 3.1 1.5 0.2 setosa
4 5.0 3.6 1.4 0.2 setosa
5 5.4 3.9 1.7 0.4 setosa
6 4.6 3.4 1.4 0.3 setosa
7 5.0 3.4 1.5 0.2 setosa
8 4.4 2.9 1.4 0.2 setosa
9 4.9 3.1 1.5 0.1 setosa
10 5.4 3.7 1.5 0.2 setosa
11 4.8 3.4 1.6 0.2 setosa
12 4.8 3.0 1.4 0.1 setosa
13 4.3 3.0 1.1 0.1 setosa
14 5.8 4.0 1.2 0.2 setosa
15 5.7 4.4 1.5 0.4 setosa
16 5.4 3.9 1.3 0.4 setosa
17 5.1 3.5 1.4 0.3 setosa
18 5.7 3.8 1.7 0.3 setosa
19 5.1 3.8 1.5 0.3 setosa
20 5.4 3.4 1.7 0.2 setosa
21 5.1 3.7 1.5 0.4 setosa
22 4.6 3.6 1.0 0.2 setosa
23 5.1 3.3 1.7 0.5 setosa
24 4.8 3.4 1.9 0.2 setosa
25 5.0 3.0 1.6 0.2 setosa
26 5.0 3.4 1.6 0.4 setosa
27 5.2 3.5 1.5 0.2 setosa
28 5.2 3.4 1.4 0.2 setosa
29 4.7 3.2 1.6 0.2 setosa
... ... ... ... ... ...
120 6.9 3.2 5.7 2.3 virginica
121 5.6 2.8 4.9 2.0 virginica
122 7.7 2.8 6.7 2.0 virginica
123 6.3 2.7 4.9 1.8 virginica
124 6.7 3.3 5.7 2.1 virginica
125 7.2 3.2 6.0 1.8 virginica
126 6.2 2.8 4.8 1.8 virginica
127 6.1 3.0 4.9 1.8 virginica
128 6.4 2.8 5.6 2.1 virginica
129 7.2 3.0 5.8 1.6 virginica
130 7.4 2.8 6.1 1.9 virginica
131 7.9 3.8 6.4 2.0 virginica
132 6.4 2.8 5.6 2.2 virginica
133 6.3 2.8 5.1 1.5 virginica
134 6.1 2.6 5.6 1.4 virginica
135 7.7 3.0 6.1 2.3 virginica
136 6.3 3.4 5.6 2.4 virginica
137 6.4 3.1 5.5 1.8 virginica
138 6.0 3.0 4.8 1.8 virginica
139 6.9 3.1 5.4 2.1 virginica
140 6.7 3.1 5.6 2.4 virginica
141 6.9 3.1 5.1 2.3 virginica
142 5.8 2.7 5.1 1.9 virginica
143 6.8 3.2 5.9 2.3 virginica
144 6.7 3.3 5.7 2.5 virginica
145 6.7 3.0 5.2 2.3 virginica
146 6.3 2.5 5.0 1.9 virginica
147 6.5 3.0 5.2 2.0 virginica
148 6.2 3.4 5.4 2.3 virginica
149 5.9 3.0 5.1 1.8 virginica

150 rows × 5 columns

In [26]:
versi = iris[iris["Species"]=="versicolor"]["petal length (cm)"]
virgi = iris[iris["Species"]=="virginica"]["petal length (cm)"]

On teste tout d’abord l’égalité des variances à l’aide de la commande :

In [27]:
scipy.stats.bartlett(versi,virgi)
Out[27]:
BartlettResult(statistic=1.2490262252046771, pvalue=0.2637385445079472)

La p-valeur valant 0.26, on ne rejette donc pas l’égalité des variances au niveau de test 5%.

On teste ensuite l’égalité des moyennes à l’aide de la commande :

In [28]:
scipy.stats.ttest_ind(versi,virgi, equal_var=True)
Out[28]:
Ttest_indResult(statistic=-12.603779441384985, pvalue=3.1788195478061495e-22)

On obtient une p-valeur égale à environ 3.2×10−22. On constate donc que l’hypothèse d’égalité des longueurs moyennes de pétales est (très facilement) rejetée à un niveau de test de 5%.

Au-delà de l’aspect gaussien que nous n’avions pas vérifié rigoureusement (on pourrait l’éprouver via des tests d’adéquation à une loi), on rejette que les longueurs des pétales d’iris des variétés et suivent la même distribution, on a en effet rejeté l’hypothèse d’égalité des moyennes.

Screencast 10 : test Khi-deux sur Weldon

Weldon a réalisé n=26306 lancers de 12 dés à 6 faces et a obtenu les résultats suivants :

In [29]:
df=pandas.DataFrame([[0,1,2,3,4,5,6,7,8,9,10,11,12],[185,1149,3265,5475,6114,5194,3067,1331,403,105,14,4,0]])
df.index=["Nombre de 5 ou 6","Effectif"]
df
Out[29]:
0 1 2 3 4 5 6 7 8 9 10 11 12
Nombre de 5 ou 6 0 1 2 3 4 5 6 7 8 9 10 11 12
Effectif 185 1149 3265 5475 6114 5194 3067 1331 403 105 14 4 0

Soit Xi la v.a.r. à valeurs dans {0,…,12} désignant le nombre de faces comportant un cinq ou un six lors du i-ème lancer.

Dans le cas où les dés sont équilibrés : XiB(13).

On teste si les dés sont équilibrés :

In [30]:
n =[185,1149,3265,5475,6114,5194,3067,1331,403,105,14,4,0]
p0 = scipy.stats.binom.pmf(range(13), 12, 1/3.0)
In [31]:
from scipy.stats import chi2_contingency
chi2_contingency([p0*float(sum(n))*100,n])
Out[31]:
(40.742379820746294,
 5.410516619909016e-05,
 12,
 array([[  2.02573723e+04,   1.21582848e+05,   3.34457039e+05,
           5.57461402e+05,   6.27099152e+05,   5.01979123e+05,
           2.92857954e+05,   1.25526959e+05,   3.92143653e+04,
           8.72959493e+03,   1.30770657e+03,   1.21582685e+02,
           4.90092871e+00],
        [  2.02573723e+02,   1.21582848e+03,   3.34457039e+03,
           5.57461402e+03,   6.27099152e+03,   5.01979123e+03,
           2.92857954e+03,   1.25526959e+03,   3.92143653e+02,
           8.72959493e+01,   1.30770657e+01,   1.21582685e+00,
           4.90092871e-02]]))

Le test n’est pas valide car les effectifs ne sont pas suffisants pour les valeurs 10, 11 et 12. On agrège ces modalités et on relance les calculs :

In [32]:
n= [185,1149,3265,5475,6114,5194,3067,1331,403,105,18]
p0 =numpy.array(list(scipy.stats.binom.pmf(range(10), 12, 1/3.0))+[sum(scipy.stats.binom.pmf(range(10,13), 12, 1/3.0))])
In [33]:
n
Out[33]:
[185, 1149, 3265, 5475, 6114, 5194, 3067, 1331, 403, 105, 18]
In [34]:
p0
Out[34]:
array([ 0.00770735,  0.04624408,  0.12717122,  0.21195203,  0.23844604,
        0.19075683,  0.11127482,  0.04768921,  0.01490288,  0.00331175,
        0.0005438 ])
In [35]:
chi2_contingency([p0*float(sum(n))*100,n])
Out[35]:
(35.130118599082756,
 0.00011863880683778717,
 10,
 array([[  2.02573723e+04,   1.21582848e+05,   3.34457039e+05,
           5.57461402e+05,   6.27099152e+05,   5.01979123e+05,
           2.92857954e+05,   1.25526959e+05,   3.92143653e+04,
           8.72959493e+03,   1.43419018e+03],
        [  2.02573723e+02,   1.21582848e+03,   3.34457039e+03,
           5.57461402e+03,   6.27099152e+03,   5.01979123e+03,
           2.92857954e+03,   1.25526959e+03,   3.92143653e+02,
           8.72959493e+01,   1.43419018e+01]]))

On rejette très largement l’hypothèse de dés équilibrés au niveau de test 5%.

Screencast 11 : test Kolmogorov-Smirnov sur Essence

On peut tester l’adéquation de la consommation d’essence à une loi normale à l’aide de Kolmogorov-Smirnov :

In [36]:
from scipy.stats import ks_2samp
In [37]:
ks_2samp(essence["conso"],list(numpy.random.normal(numpy.mean(essence["conso"]), numpy.std(essence["conso"]), 1000)))
Out[37]:
Ks_2sampResult(statistic=0.049562499999999954, pvalue=0.9375919473764821)

On ne peut donc pas rejetter l’hypothèse de normalité au niveau de test 5%.

On aurait également pu lancer un test de Shapiro-Wilk :

In [38]:
scipy.stats.shapiro(essence["conso"])
Out[38]:
(0.9871523380279541, 0.2743317186832428)

La conclusion est identique au niveau de test 5%.

Notons que le test de Shapiro-Wilk est recommandé pour tester la normalité dans le cas de petits échantillons.

In [ ]: