Découvrir Scilab/Statistiques

4. Statistiques

Pour l'analyse statistique de données, on peut se reporter au document

(en) Anna Bassi et Giovanni Borzi, « Data Analysis and Statistics: a Scilab data mining tutorial » [PDF], sur Openeering.

Tendance centrale

Considérons une matrice X. On peut calculer les moyennes suivantes :

mean(X) : calcule la moyenne de tous les éléments de X ;
mean(X, "r") : donne un vecteur ligne (row) contenant les moyennes colonne par colonne ;
mean(X, "c") : donne un vecteur colonne contenant les moyennes ligne par ligne.

On peut calculer la moyenne géométrique avec geomean(), et la moyenne harmonique avec harmean(), avec la même syntaxe.

De la même manière, on peut calculer les médianes avec median(X), median(X, "r") et median(X, "c").

La commande trimmean(X, p) permet de calculer la moyenne en éliminant les p pourcents extrêmes (les p/2 % valeurs les plus élevées, et les p/2 % valeurs les plus faibles). On peut également faire le calcul par ligne ou par colonne — trimmean(X, p, "r") ou trimmean(X, p, "c"). notons que l'on a :

mean(X) = trimmean(X, 0) ;
median(X) = trimmean(X, 100).

Dispersion

La commande

strange(X)

calcule l'écart (range) entre la valeur maximale et la valeur minimale de la matrice. On peut calculer les valeurs pour chaque colonne ou chaque ligne avec strange(X, "r") et strange(X, "c"). La commande

mad(X)

calcule la déviation absolue moyenne (mean absolute deviation)

{\dfrac {1}{n}}\sum |x_{i}-{\bar {x}}|

on peut de même la calculer pour chaque colonne ou ligne.

La commande

iqr(X)

calcule l'écart interquartile (interquartile range).

La commande stdev calcule l'écart type (standard deviation) et la commande variance la variance des éléments de la matrice (entière, par ligne ou par colonne). Par défaut, la variance est sans biais : normalisée par m × n – 1 pour une matrice de m × n éléments, par m – 1 par colonne et par n – 1 par ligne. Si l'on calcule la variance par ligne ou par colonne, on peut ajouter l'option 1 pour calculer le moment d'ordre 2 (normalisation par m ou n), par exemple

mom2 = variance(X, "r", 1) // pour une matrice non-vecteur
mom2 = variance(X, "*", 1) // pour un vecteur

L'écart type est toujours sans biais. Si l'on veut l'estimateur normalisé par m ou n, il faut utiliser

estim_ecart_type = sqrt(variance(X, "*", 1)) // pour une matrice non-vecteur

Normalisation

La commande center(X) permet de centrer la matrice (équivalent de X - mean(X)). La commande wcenter(X) donne la matrice centrée réduite (équivalent de 1/stdev(X)*(X - mean(X))).

Fréquences et quartiles

La commande tabul(X) calcule les fréquences d'apparition des valeurs dans la matrice X. Le résultat est classé par défaut par ordre décroissant des valeurs ; pour avoir un classement croissant, on ajoute le paramètre "i" (increasing order) :

tabul(X, "i")

La commande quart(X) calcule les quartiles et les renvoie dans un vecteur de dimension 3. On peut calculer les quartiles pour chaque colonne ou chaque ligne, avec quart(X, "r") et quart(X, "c").

On peut aussi travailler avec des valeurs et leur fréquence. Si la matrice X contient les valeurs et la matrice F leur fréquence — F(i, j) est la fréquence de X(i, j) —, alors

meanf(X, F) donne la moyenne ;
stdevf(X, F) donne l'écart type ;
variancef(X, F) donne la variance.

Notons que F n'est pas nécessairement normalisé à 1 ou à 100 ; il peut par exemple s'agir d'un nombre d'occurrences.

Exemple

X = int(10*rand(100,1, "normal")); // génère 100 valeurs entières avec une loi normale centrée d'écart type 10
F = tabul(X); // calcule les fréquences d'apparition
x = F(:, 1); // extraction de la liste des valeurs
f = F(:, 2); // extraction de la liste des fréquences

quartiles = quart(X); // calcule les quartiles de la population

moyenne = meanf(x, f); // calcul de la moyenne pondérée par les fréquences

Considérons maintenant deux vecteurs X et Y de même dimension, et la matrice F des fréquences de X et Y : F(i, j) est la fréquence de X(i) & Y(j). Alors :

covar(X, Y, F) donne la covariance ;
correl(X, Y, F) donne le coefficient de corrélation.

Exemple

n = 10; // nombre de valeurs
X = 1:n; // abscisses
Y = 2*X + 5 + 3*rand(X); // ordonnées : droite "bruitée"
plot(X, Y); // tracé des données
F = eye(n, n); // couplage (x, y) : chaque couple est présent une fois, (X(i) & Y(i)) = 1, (X(i) & Y(j)) = 0

covariance1 = covar(X, Y, F);
covariance2 = 1/n*sum(X.*Y) - 1/n^2*sum(X)*sum(Y); // calcul "à la main" pour comparer

correlation1 = correl(X, Y, F);
correlation2 = covariance2/sqrt(variance(X, "*", 1))/sqrt(variance(Y, "*", 1)); // calcul "à la main"

Fonctions de distribution et de répartition de lois statistiques

Par défaut, Scilab propose une loi de distribution, pour la loi binomiale, et onze lois de répartition (probabilités cumulatives).

Pour la loi binomiale, la fonction de répartition pour n épreuves (nombre entier) avec une probabilité p (nombre réel entre 0 et 1) est donnée par

f = binomial(p, n)

La variable f est un vecteur de n + 1 composantes, avec

{\mathtt {f}}(k+1)={\begin{pmatrix}n\\k\end{pmatrix}}\cdot p^{k}\cdot (1-p)^{n-k}=\mathrm {C} _{n}^{k}\cdot p^{k}\cdot (1-p)^{n-k}

Notez que les coefficients binomiaux C^k_n s'obtiennent eux avec la commande b = nchoosek(n, k).

La fonction de répartition s'obtient avec cdfbin() :

P =  = cdfbin("PQ", k, n, p, 1-p)

donne la probabilité P d'avoir entre 0 et k succès pour n tirages avec une probabilité de réussite p. On peut d'ailleurs le comparer avec

f = binomial(p, n)
c = sum(f(1:k+1))

De manière générale, les fonctions de répartition commencent par cdf… (cumulative distribution function), par exemple :

P = cdfchi("PQ", x, k) : loi du χ² à k degrés de liberté ;
P = cdfnor("PQ", x, mu, sigma) : loi normale d'espérance mu et d'écart type sigma ;
P = cdfpoi("PQ", k, lambda) : loi de Poisson de paramètre lambda ;
P = cdft("PQ", t, k) : loi t de Student à k degrés de liberté ;
… voir Scilab help >> Statistics > Cumulated Distribution Functions

Notons que pour chaque cas, on peut aussi avoir la valeur complémentaire Q = 1 – P avec la syntaxe

[P, Q] = cdf…("PQ", …)

Exemple

Ci-dessous, nous traçons la fonction de répartition de la loi normale centrée réduite (μ = 0, σ = 1). Notons que pour des calculs vectoriels, l'espérance et l'écart type sont des vecteurs de même dimension que X.

X = -3:0.1:3;
mu = zeros(X);
sigma = ones(X);
[P, Q] = cdfnor("PQ", X, mu, sigma);
plot(X', [P', Q'])

Les mêmes fonctions permettent de rechercher les quantiles ou fractiles :

k = cdfbin("S", n, P, 1-P, p, 1-p) : nombre de succès k sur n épreuves ayant pour paramètre p, donnant une probabilité cumulée P ;
x = cdfchi("X", k, P, 1-P) : valeur de x pour laquelle on a une probabilité cumulée P avec une loi du χ² à k degrés de liberté ;
x = cdfnor("X", mu, sigma, P, 1-P) : valeur de x pour laquelle on a une probabilité cumulée P avec une loi normale d'espérance mu et d'écart type sigma ;
k = cdfpoi("S", lambda, P, 1-P) : valeur de k pour laquelle on a une probabilité cumulée P avec une loi de Poisson de paramètre lambda ;
t = cdft("T", k, P, 1-P) : valeur de t pour laquelle on a une probabilité cumulée P avec une loi t de Student à k degrés de liberté ;
…

Exemple

Ci dessous, nous déterminons le quantile pour un risque α bilatéral de 5 %, c'est-à-dire que l'on rejette 0,25 % des cas les plus hauts et les plus bas, pour une loi normale centrée réduite puis pour une loi t de Student à 30 degrés de liberté

-->x = cdfnor("X", 0, 1, 0.975, 0.025)
 x  =
 
    1.959964
 
-->t = cdft("T", 30, 0.975, 0.025)
 t  =
 
    2.0422725

Enfin, ces fonctions permettent de déterminer les paramètres des loi pour atteindre un quantile visé :

loi binomiale :
- n = cdfbin("Xn", p, 1-p, P, 1-P, k) : nombre d'épreuves ayant pour paramètre p afin d'avoir un quantile P valant k,
- [p, q] = cdfbin("PrOmpr", P, 1-P, k, n) : paramètre p de chaque épreuve afin d'avoir un quantile P pour k succès parmi n ;
loi du χ² :
k = cdfchi("Df", P, 1-P, x) : valeur du nombre de degrés de liberté (degrees of freedom) k pour avoir un quantile P valant x ;
loi normale :
- mu = cdfnor("Mean", sigma, P, 1-P, x) : valeur de l'espérance (moyenne) pour avoir un quantile P valant x si l'on a un écart type sigma,
- sigma = cdfnor("Std", P, 1-P, X, mu) : valeur de l'écart type (standard deviation) pour avoir un quantile P valant x si l'on a une espérance de mu ;
loi de Poisson :
lambda = cdfpoi("Xlam", P, 1-P, k) : valeur du paramètre λ pour avoir un quantile P valant k ;
loi t de Student :
k = cdft("Df", P, 1-P, t) : nombre de degrés de liberté k pour avoir un quantile P valant t.

Générateurs pseudo-aléatoires

Les générateurs de nombres aléatoires permettent de tester des programmes ou bien de faire des simulations par la méthode de Monte Carlo.

Scilab possède une fonction rand() qui par défaut génère un nombre entre 0 et 1 avec une loi de distribution uniforme. Par exemple, pour simuler un jet de dé à six faces, on peut utiliser :

de = 1 + int(6*rand())

La fonction peut générer une matrice m × n avec les syntaxes

rand(m, n)
rand(X)

où X est une matrice m × n. Notons que si l'on veut générer n valeurs, il faut bien taper rand(n, 1) ou rand(1, n).

On peut lui faire générer une valeur selon une loi normale centrée réduite avec la clef "normal" :

rand(m, n, "normal")

Exemple

Supposons par exemple que l'on a une poutre faite dans un acier dont la résistance intrinsèque, la limite d'élasticité R_e, est garantie au minimum à 235 MPa. En réalité, R_e suit une loi normale d'espérance μ_Re = 269 MPa et d'écart type σ_Re = 17 MPa. Dans les conditions d'utilisation, l'acier de cette poutre doit résister à une contrainte S (stress) de 214 MPa ; en fait, les conditions de chargement sont variables et suivent une loi normale de paramètres μ_S = 214 MPa et σ_S = 15 MPa.

On veut déterminer statistiquement dans combien de cas la contrainte S est supérieure à la résistance R_e.

n = 1e6; // nombre d'événements
R = 269 + 17*rand(n, 1, "normal"); // résistance
S = 214 + 14*rand(n, 1, "normal"); // contrainte
Mbool = S>R; // comparaison des valeurs (matrice de booléens)
moyenne = sum(Mbool)/n; // les "vrai" comptent pour 1 et les "faux" pour 0
disp(string(100*moyenne)+" %") // affichage du résultat

La commande grand() (get random number) permet d'avoir des nombres aléatoires selon de nombreuses lois, par exemple :

Y = grand(m, n, "bin", N, p) : matrice m × n de nombres aléatoires tirés selon une loi binomiale de N essais avec une probabilité p chacun ;
Y = grand(m, n, "nor", mu, sigma) : pour une loi normale d'espérance mu et d'écart type sigma ;
Y = grand(m, n, "unf", a, b) : pour une loi uniforme sur [a ; b[ (b exclus) ;
Y = grand(m, n, "uin", a, b) : pour une loi uniforme sur [a ; b] (b inclus) ;
…

Dans l'exemple précédent, on aurait donc :

R = grand(n, 1, "nor", 369, 17); // résistance
S = grand(n, 1, "nor", 314, 14); // contrainte

Autres lois statistiques

Des modules complémentaires Atoms permettent d'avoir accès à d'autres loi.

Par exemple le module CASCI donne accès à de nombreuses fonctions de distribution et de répartition, notamment des lois log-normale (pdflognormal(), cdflognormal()), exponentielle (pdfexponential(), cdfexponential()) ou de Weibull (pdfweibull(), cdfweibull()). Ainsi :

p = pdfnormal(x, mu, sigma); // densité d'une loi normale en x
y = cdflognormal(x, lambda); // répartition d'une loi log-normale d'espérance mu et d'écart type sigma
y = cdfexponential(x, mu, sigma); // répartition d'une loi exponentielle de paramètre lambda
y = cdfweibull(x, k, lambda, theta); // répartition d'une loi de Weibull
   // de paramètre de forme k, de paramètre d'échelle lambda et de paramètre de position thêta

Le module fournit également le calcul des quantiles avec les fonctions idf… (inverse cumulative distribution functions), et des générateurs de nombres aléatoires avec les fonctions rnd…. Par exemple :

x = idfweibull(y, k, lambda, theta); // quantile y d'une loi de Weibull
   // de paramètre de forme k, de paramètre d'échelle lambda et de paramètre de position thêta
X = rndweibull(n, k, lambda, theta); // génère n nombres aléatoires selon la loi de Weibull

Module Stixbox

Le module Stixbox définit également un certain nombre de fonctions statistiques utiles. Citons par exemple la fonction quantile() qui permet de calculer les quantiles expérimentaux et la fonction qqplot() qui permet de tracer un diagramme quantile-quantile (comme une droite de Henry). Par exemple :

n = 100;
X = rand(1, n, "norm"); // tirage aléatoire de n valeurs
// avec une loi normale centrée réduite

p = 0.1:0.1:0.9;
q = 1 - p;

qexp = quantile(X, p); // déciles expérimentaux

for i = 1:size(p, "*") // déciles théorique pour une loi normale centrée réduite
    qth(i) = cdfnor("X", 0, 1, p(i), q(i));
end

scf(0);
clf();
qqplot(qth', qexp); // diagramme quantile-quantile
xtitle("diagramme quantiles-quantiles", "quantiles théoriques",...
    "quantiles expérimentaux");

Notes

Calcul numérique < ↑ > Interpolation, extrapolation et lissage