Statistique et probabilité

Ce tutorial montre comment génerer puis analyser des réalisations aléatoire en Python.

Home Statistics .ipynb
In [1]:
%matplotlib inline
from numpy import *
from scipy.stats import *
from matplotlib.pyplot import *

Dans ce tutorial, nous allons apprendre à utiliser le couteau suisse Python pour la génération et l'analyse des variables aléatoires. De base, Python ne connait qu'un nombre limité de distributions aléatoires. Toutefois, couplé avec le module stats de Scipy, Python rivalise fièrement avec n'importe quel outil numérique spécialisé en statistique comme R, Matlab ou Scilab. Il suffit en effet de jeter un coup d'oeil rapide à la liste des distributions implémentées dans Scipy pour se rendre compte que Python n'a pas à rougir face à ses concurrents !

Cas des Variables Aléatoires Discrêtes

Une variable aléatoire discrête est une variable qui prend un nombre dénombrable de valeurs. La loi de probabilité d'une variable aléatoire discrête $X$ est donnée par :

$$P(X=x)$$

où $x \in \chi$ correspond à une valeur possible pour $X$ et $\sum_{x \in \chi}P(X=x)=1$. Dans ce tutorial, nous allons prendre comme exemple le cas de la loi binomiale $\mathcal{B}(n,p)$ de paramêtres $n=10$ et $p=0.4$. La loi de probabilité d'une variable aléatoire $X\sim \mathcal{B}(n,p)$ est donnée pour tout $0\le k \le n$ par:

$$P(X=k)=C_{n}^{k}p^k (1-p)^{n-k}$$

Pour générer une réalisation de $N=1000$ échantillons suivant une loi binomiale $\mathcal{B}(10,0.4)$, nous allons utiliser l'objet binom() du module stats de Scipy (voir documentation). Le constructeur de cet objet prend deux paramètres: $n$ et $p$. La génération d'une réalisation de $N$ échantillons s'obtient ensuite en appelant la méthode rvs() avec l'argument size=1000.

In [2]:
loi = binom(10,0.4)      # creation de l'objet
X = loi.rvs(size=1000) # réalisation de 1000 échantillons
print(X[:10])         # affichage des 10 premiers élements
[5 2 5 4 3 3 6 1 6 6]

Affichage de l'histogramme

Pour vérifier que notre réalisation suit bien une loi binomiale, nous allons afficher l'histogramme des $1000$ échantillons. Pour obtenir l'histogramme, nous allons utiliser la fonction hist() de Matplotlib (voir Documentation). Cette fonction prend en entrée les échantillons $X$. Cette fonction peut également prendre un argument optionnel bins permettant de spécifier les intervalles dans lesquels nous souhaitons évaluer l'histogramme. Pour spécifier ces intervalles, nous allons utiliser la fonction arange(min,max) de Python. Cette fonction permet de générer une suite de nombres allant de min à max-1 par pas de $1$.

In [3]:
bins = arange(-0.5,10.5)
print(bins)
hist(X,bins);
[-0.5  0.5  1.5  2.5  3.5  4.5  5.5  6.5  7.5  8.5  9.5]

Remarquons que l'affichage des courbes avec Matplotlib est facilement "customisable". D'une part, il est possible de modifier l'appel de la fonction hist pour préciser que nous souhaitons explicitement normaliser la somme des élements à $1$ (cas d'une distribution aléatoire). D'autre part, il est possible de changer l'apparence graphique de l'histogramme en modifiant la couleur, la taille des traits et les annotations.

In [4]:
hist(X,bins=bins,density=True, color='w', edgecolor="b")
xlabel("$x$")
ylabel("$P[X=x]$");

Nous allons maintenant superposer à cette courbe la loi théorique de notre variable aléatoire via un diagramme en bâton. La loi théorique s'obtient en appelant la méthode pmf() de nottre objet loi. Cette méthode prend comme argument les différentes valeurs de $x$. L'affichage du diagramme en bâton s'obtient ensuite en utilisant la commande stem de Matplotlib (voir Documentation).

In [5]:
x = np.arange(11)
hist(X,bins=bins,density=True, color='w', edgecolor="b",label="theorique")
stem(x,loi.pmf(x),"r",markerfmt='ro',label="theo")
xlabel("$x$")
ylabel("$P[X=x]$")
legend();

C'est plutôt pas mal ! A noter qu'il est possible d'obtenir un meilleur fit en augmentant le nombre d'échantillons $N$.

Affichage de la Fonction de Répartition

La fonction de répartition théorique $F$ d'une loi aléatoire est définie par la fonction:

$$F(t)=P(X\le t)$$

où $t\in \mathbb{R}$ et $F(t)\in [0,1]$. En Python, Il est possible de calculer explicitement cette fonction de répartition en appelant la fonction cumsum(X). Pour aller plus vite, il est également possible d'afficher directement la fonction de répartition en ajoutant l'argument cumulative=True à la fonction hist de Matplotlib. A noter que dans le cas des variables discrêtes, nous allons évaluer $F(t)$ pour $t=x-\epsilon$ où $\epsilon$ est une valeur très petite.

In [6]:
eps = 10**(-10)
bins = arange(0-eps,10-eps)
hist(X,bins=bins,density=True,cumulative=True,color='w', edgecolor="b",)
xlabel("$t$")
ylabel("$F(t)$");

Pour superposer la fonction de répartition théorique, nous allons utiliser la méthode cdf() de notre objet loi. La diagramme en escalier s'obtient ensuite en utilisant la fonction step de Matplotlib.

In [7]:
hist(X,bins=bins,density=True,cumulative=True,color='w', edgecolor="b",label="exp")
step(bins,loi.cdf(bins),"r",linestyle='--', label="theo") 
xlabel("$t$")
ylabel("$F(t)$")
legend();