Python pour le calcul scientifique/Statistiques

Le module NumPy fournit des fonctions statistiques.

Rappelons que dorénavant les programmes commencent tous par :

#!/usr/bin/python3

import numpy as np
import matplotlib.pyplot as plt

Méthodes de matrices

Rappel : la classe des matrices (ndarray) dispose de méthodes permettant de calculer des statistiques sur les éléments des matrices :

• .min() : minimum des valeurs ;
• .max() : maximum des valeurs ;
• .ptp() : écart amplitude « max – min » (peak to peak) ;
• .mean() : moyenne ;
• .std() : écart type (standard deviation).

Statistiques descriptives

NumPy fournit la fonction np.quantile() qui détermine les quantiles avec la syntaxe :

np.quantile(M, q)

où M est une matrice (ou une liste, un n-uplet, bref un itérable de nombres) et q est un quantile ou un vecteur de quantiles sous la forme d'un nombre entre 0 et 1. Par exemple, pour avoir les quartiles :

np.quantile(M, [0, 0.25, 0.5, 0.75, 1])

Si un des éléments de la matrice est un NaN, le résultat est un NaN. Pour éviter cela, on peut utiliser la fonction np.nanquantile() qui ignore les NaN.

Les fonctions np.percentile() et np.nanpercentile() donnent les centiles ; on indique alors le centile que l'on veut sous la forme d'un nombre entier entre 0 et 100.

Nous disposons également des fonctions suivantes :

• np.amin(), np.nanmin() : minimum ;
• np.amax(), np.nanmax() : maximum ;
• np.ptp() : amplitude.

Notons que pour toutes les fonctions, il est possible d'indiquer l'axe selon lequel on effectue le calcul. Par exemple,

M = np.arange(9).reshape(3, 3)
# [[0, 1, 2]
#  [3, 4, 5]
#  [6, 7, 8]]
print(np.quantile(M, [0.25, 0.5], 0))
# [[1.5 2.5 3.5]  : 1er quartile des colonnes
#  [3.  4.  5. ]] : médiane des colonnes
print(np.quantile(M, [0.25, 0.5], 1))
# [[0.5 3.5 6.5]  : 1er quartile des lignes
#  [1.  4.  7. ]] : médiane des lignes

Indicateurs de position et de dispersion

NumPy fournit les indicateurs de tendance centrale suivants :

• np.mean() et np.nanmean() : moyenne ;
• np.average() et np.nanaverage() : moyenne pondérée ; la syntaxe est np.average(M, axe, poids) ou bien np.average(M, weights = poids) (voir ci-après) ;
• np.median() et np.nanmedian() : médiane.

En terme de performances, la fonction np.mean() est équivalente à la méthode M.mean() et à la fonction np.average() sans poids (ce qui équivaut donc à la moyenne). En revanche, la fonction np.nanmean() est plus lente, de même que la fonction np.average() lorsque l'on utilise des poids (même s'ils sont tous égaux à 1).

Comme précédemment, on peut indiquer l'axe (si l'on veut évaluer les valeurs par colonne ou par ligne). Pour la moyenne pondérée, on utilise une matrice poids P de même dimension que la matrice de valeurs M, P[i, j] étant le poids associé à la valeur M[i, j]. Par exemple, si l'on veut évaluer la moyenne pondérée pour toutes les valeurs de M (pas d'axe), on peut écrire une des deux solutions suivantes :

np.average(M, None, P) # l'axe est le 2e paramètre
np.average(M, weights = P)

NumPy fournit également les indicateurs de dispersion suivants :

• np.std(), np.nanstd() : écart type (standard deviation) ;
• np.var(), np.np.nanvar() : variance.

Fréquence, histogramme

On peut générer une matrice aléatoire avec les fonctions np.random.rand(), qui utilise une loi uniforme sur [0 ; 1], et np.random.randn() qui utilise une loi normale centrée réduite.

Lorsque l'on dispose d'une série de données aléatoires, qu'elles aient été mesurées ou bien générées par une fonction aléatoire, on peut ensuite les mettre dans des classes (bins). Les classes sont définies par un vecteur (ou une liste, un n-uplet) [c1, c2,c3, …, cn]. La classe 0 désigne les valeurs inférieures à c₁ ; la classe 1 désigne les valeurs c₁ ≤ x < c₂ ; la classe n (n == len(classes)) désigne les valeurs supérieures à c_n.

La fonction np.digitize() indique dans quelle classe se trouve un nombre. Si on lui donne une matrice (ou un vecteur, une liste, un n-uplet), il renvoie une matrice de même dimension, np.digitize(M, classes)[i, j] étant la classe de l'élément M[i, j]. Par exemple :

M = np.random.rand(10)

classes = (0, 0.2, 0.4, 0.6, 0.8, 1)

print(M)
print(np.digitize(M, classes))

On peut changer la « largeur » des inégalités avec le paramètre right = True : right, « droite », est vrai (true) lorsque l'inégalité large est à droite, la classe i désigne les valeurs c_i < x ≤ c_{i + 1}. La syntaxe est alors np.digitize(M, classes, True) ou bien np.digitize(M, classes, right = True). Le sens de l'inégalité ne dépend pas du sens du vecteur de classes ; ainsi, right désigne la valeur supérieure de l'intervalle même si le vecteur de classes est classé par ordre décroissant.

Histogramme tracé avec Python/Numpy/Matplotlib.

La fonction np.histogram() détermine n classes de même taille et renvoie deux vecteurs : le premier contient le nombre d'éléments dans chaque classe et le second décrit les classes (c'est-à-dire les bornes des classes).

Pour tracer l'histogramme, nous disposons de la fonction plt.hist() :

M = np.random.randn(50)
plt.plot(M, np.ones_like(M), "|")
plt.hist(M, bins=10, density=1)

La fonction np.bincount() travaille sur les listes de nombres entiers. Elle renvoie un n-uplet, np.bincount(M)[i] est le nombre de fois que le nombre i revient dans la matrice M — rappel, le premier élément du n-uplet np.bincount(M) a l'indice 0 donc correspond au nombre de fois que le nombre 0 apparaît dans la matrice M.

Lois de probabilités

Avec NumPy

NumPy fournit des générateurs de nombres pseudo-aléatoires, avec le module random.

Pour définir un générateur de nombres (pseudo-)aléatoires (GNA), on l'initialise avec la fonction

gna = np.random.default_rng()

on peut lui indiquer une graine, qui est un nombre entier, idéalement grand. Par exemple :

gna = np.random.default_rng()

On peut ensuite utiliser ce générateur pour effectuer des tirages selon différentes lois de probabilité :

• .random() : distribution uniforme dans [0 ; 1[ ;
• .standard_normal() : distribution normale centrée réduite ;
• .integers(a, b) : nombres entiers répartis uniformément sur [a ; b [.

Par exemple :

# loi uniforme
A = gna.random() # un nombre aléatoire
B = gna.random(5) # vecteur de 5 nombres aléatoires ; équivalent à gna.random(size = 5)

# loi normale
A = gna.standard_normal() # un nombre aléatoire
B = gna.standard_normal(5) # vecteur de 5 nombres aléatoires ; équivalent à gna.standard_normal(size = 5)

# loi uniforme sur des entiers
A = gna.integers(1, 10) # un nombre aléatoire entre 1 et 9 inclus ; équivalent à gna.integers(low = 1, high = 10)
A = gna.integers(1, 10, 5) # 5 nombre aléatoires entre 1 et 9 inclus ; équivalent à gna.integers(low = 1, high = 10, size = 10)

NumPy fournit des générateurs pseudo-aléatoires pour les lois suivantes :

.beta(a, b), .binomial(n, p), .chisquare(df), dirichlet(alpha), .exponential(), .f(dnum, dfden) (F), .gamma(shape), .geometric(p), .gumbel(), .hypergeometric(ngood, nbad, nsample), .laplace(), .logisticl(), .lognormal(), .logseries(p), .multinomial(n, pvals), .multivariate_hypergeometric(colors, nsample), .multivariate_normal(mean, cov), .negative_binomial(n, p), .noncentral_chisquare(df, nonc), .noncentral_f(dfnum, dfden, nonc), .normal(), .pareto(a), .poisson(), .power(a), .rayleigh(), .standatd_cauchy(), .standard_exponential(), .standard_normal(), .standard_t(df) (Student), .triangular(left, mode, right), .uniform(low, high, size), .vonmises(mu, kappa), .wald(mean, scale) (Gauss inverse), .weibull(a) et .zipf(a).

Nous n'avons indiqué que les paramètres obligatoires. Elles possèdent toues des paramètres optionnels, dont size (1 par défaut) qui donne la taille du vecteur généré.

Voir

• « Random Generator > Distributions », sur NumPy.org

Le module np.random permet également :

• de prendre un ou plusieurs éléments au hasard dans une matrice M : .choice(M)
gna.choice(M) pour un élément, gna.choice(M, n) ou gna.choice(M, size=n) pour n éléments ;
• de mélanger un matrice M :
  • .shuffle(M) mélange la matrice M,
  • .permutation(M) crée une copie mélangée de la matrice M ;
• on peut indiquer l'axe de mélange (axe 0 par défaut), par exemple pour l'axe 1 : .shuffle(M, 1), .shuffle(M, axis = 1), .permutation(M, 1), .permutation(M, axis=1).

Pour assurer une transition simple avec Matlab, NumPy fournit les générateurs pseudo-aléatoires suivants :

• np.random.rand() pour la loi uniforme sur [0 ; 1[ ; on peut indiquer les dimension d'une matrice (ou hypermatrice) en paramètres, par exemple np.random.rand(a, b) pour une matrice a × b ;
• np.random.randn() pour la loi normale centrée réduite ; on peut de même indiquer les dimension d'une matrice (ou hypermatrice) en paramètres, par exemple np.random.rand(a, b) pour une matrice a × b.

Avec SciPy

Le module scipy.stats fournit un grand nombre de lois de probabilités. Ces lois sont des objets ayant toutes les mêmes méthodes. Par exemple, on dispose des lois suivantes :

• lois dicrètes :
  • bernoulli : loi de Bernoulli,
  • scipy.stats.binom : loi binomiale,
  • scipy.stats.poisson : loi de Poisson,
  • scipy.stats.randint : loi discrète uniforme ;
• loi continues :
  • scipy.stats.chi2 : loi du χ² (khi carré),
  • scipy.stats.norm : loi normale,
  • scipy.stats.t : loi de Student,
  • scipy.stats.uniform : loi continue uniforme ;

et des méthodes suivantes :

• .mean() : calcule la moyenne, l'espérance de la loi ;
• .median() : calcule la médiane de la loi ;
• .var() : calcule la variance de la loi ;
• .std() : calcule l'écart type (standard deviation) de la loi ;
• .moment() : calcule les moments de la loi ;
• .rvs() : effectue des tirages aléatoires (random variates) ;
• .pdf() : fonction densité de probabilité (probability density function) pour les lois continues ;
• .mdf() : fonction de masse (mass density function) pour les lois discrètes ;
• .cdf() : fonction de répartition (cumulative density function) ;
• .ppf() : fonction quantile (percent-point function), inverse de .cdf() ;
• .fit() : calcule les paramètres de la loi correspondant au mieux à un échantillon, au sens du maximum de vraisemblance ; uniquement pour les lois continues.

Par exemple, pour la loi normale, on a :

${\mathtt {scipy.stats.norm.mean(m,sigma)}}=m$

${\mathtt {scipy.stats.norm.pdf(x,m,sigma)}}={\frac {1}{\sigma {\sqrt {2\pi }}}}\exp \left(-{\frac {1}{2}}\left({\frac {x-m}{\sigma }}\right)^{2}\right)$

${\mathtt {scipy.stats.norm.cdf(x,m,sigma)}}=\int _{-\infty }^{x}{\frac {1}{\sigma {\sqrt {2\pi }}}}\exp \left(-{\frac {1}{2}}\left({\frac {t-m}{\sigma }}\right)^{2}\right)\mathrm {d} t$

Exemple de mise en œuvre :

from scipy import stats

# moments mvsk : mean (moyenne), variance, skewness (asymétrie), kurtosis
moy, var, asy, kurt = stats.norm.stats(moments="mvsk")
print(moy, var, asy, kurt)
# 0.0 1.0 0.0 0.0

print(stats.norm.pdf(1), ";", stats.norm.cdf(0.5), ";", stats.norm.ppf(0.99))
# Pour la fonction normale centrée réduite : densité de probabilité en 1 ; fonction de répartition en 0,5 ; dernier centile
# 0.24197072451914337 ; 0.6914624612740131 ; 2.3263478740408408

x = np.linspace(-6, 6, 100)
y1 = stats.norm.pdf(x) # Densité de probabilité de la loi normale centrée réduite
y2 = stats.norm.cdf(x) # Fonction de répartition de cette loi

plt.plot(x, y1, label="densité")
plt.plot(x, y2, label="répartition")
plt.legend()

Par défaut, les lois continues sont centrées réduites (de moyenne nulle et d'écart type unité). Les objets disposent des paramètres loc (location, position) et scale (échelle) ; par exemple pour un tirage aléatoire de 100 échantillons avec une loi normale de moyenne 10 et d'écart type 5 :

stats.norm.rvs(100, loc=5, scale=10)

Il est possible de fixer ces paramètres pour une loi, de la « geler », par exemple :

from scipy import stats

varAl = stats.norm(loc=5, scale=10) # loi « gelée » (loi normale d'espérance 5 et d'écart type 10)
n = 100
x = varAl.rvs(n)
plt.plot(x, np.ones_like(x), "|")
plt.hist(x)
print(varAl.mean(), varAl.std())

ou encore

from scipy import stats

n = 100
varAl = stats.poisson(mu=50) # loi de Poison d'espérance mu
x = varAl.rvs(size=n)
plt.plot(x, np.ones_like(x), "|")
plt.hist(x)
print(varAl.mean(), varAl.std())

res = stats.norm.fit(x)
print(res)

Diagramme quantile-quantile avec Python Scipy et Matplotlib.

L'exemple suivant consiste à vérifier qu'un échantillon suit bien une loi normale centrée réduite en traçant sa droite de Henry :

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

X = stats.norm.rvs(size=100)

p = np.arange(0.1, 1, 0.1)

qexp = np.quantile(X, p)
qth = stats.norm.ppf(p)
extr = stats.norm.ppf([0.1, 0.9])

plt.plot(extr, extr, "--r")
plt.plot(qth, qexp, "*")
plt.title("diagramme quantiles-quantiles")
plt.xlabel("quantiles théoriques")
plt.ylabel("quantiles expérimentaux")

plt.savefig("diagramme_qq_python_matplotlib.svg")

Pour les lois discrète, il est nécessaire d'indiquer leurs paramètres. Par exemple :

• pour la loi de Bernoulli, la probabilité p : stats.bernoulli(p) ;
• pour la loi binomiale, les facteurs de forme n et p : stats.binom(n, p) ;
• pour la loi de Poisson : le paramètre μ (mu) : stats.poisson(mu) ; l est également possible de décaler la loi, stats.poisson(mu, loc) ;
• pour la loi discrète uniforme, les valeurs minimales et maximales : stats.randint(low, high).

Ressources

• « Statistical functions (scipy.stats) », sur Numpy and Scipy Documentation (consulté le 22 mai 2019)

Notes et références

---

Polynômes < ↑ > Interpolation, extrapolation et lissage