Comment calculer efficacement un écart-type courant

Question

J'ai un tableau de listes de nombres, par exemple :

[0] (0.01, 0.01, 0.02, 0.04, 0.03)
[1] (0.00, 0.02, 0.02, 0.03, 0.02)
[2] (0.01, 0.02, 0.02, 0.03, 0.02)
     ...
[n] (0.01, 0.00, 0.01, 0.05, 0.03)

Je voudrais calculer efficacement la moyenne et l'écart type à chaque indice d'une liste, sur tous les éléments du tableau.

Pour faire la moyenne, j'ai parcouru le tableau en boucle et j'ai additionné la valeur à un index donné d'une liste. A la fin, je divise chaque valeur de ma "liste de moyennes" par n (Je travaille avec une population, pas un échantillon de la population).

Pour calculer l'écart-type, je repasse en boucle, maintenant que j'ai calculé la moyenne.

Je voudrais éviter de parcourir deux fois le tableau, une fois pour la moyenne et une fois pour l'écart-type (après avoir obtenu une moyenne).

Existe-t-il une méthode efficace pour calculer les deux valeurs, en ne parcourant le tableau qu'une seule fois ? Tout code dans un langage interprété (par exemple, Perl ou Python) ou pseudocode est acceptable.

Answer 1

La réponse est d'utiliser l'algorithme de Welford, qui est très clairement défini après les "méthodes naïves" dans :

• Wikipedia : Algorithmes de calcul de la variance

Elle est numériquement plus stable que les collecteurs à deux passages ou la somme simple des carrés en ligne suggérés dans d'autres réponses. La stabilité n'a d'importance que lorsque vous avez beaucoup de valeurs proches les unes des autres, car elles conduisent à ce que l'on appelle la " annulation catastrophique " dans la littérature sur la virgule flottante.

Vous voudrez peut-être aussi revoir la différence entre la division par le nombre d'échantillons (N) et N-1 dans le calcul de la variance (écart au carré). La division par N-1 conduit à une estimation non biaisée de la variance de l'échantillon, alors que la division par N sous-estime en moyenne la variance (parce qu'elle ne tient pas compte de la variance entre la moyenne de l'échantillon et la vraie moyenne).

J'ai écrit deux articles de blog sur le sujet qui donnent plus de détails, notamment sur la manière de supprimer les valeurs précédentes en ligne :

• Calcul en ligne de la moyenne et de la variance d'un échantillon en une seule passe
• Suppression de valeurs dans l'algorithme de Welford pour la moyenne et la variance en ligne

Vous pouvez également jeter un coup d'œil à mon implémentation Java ; la javadoc, les sources et les tests unitaires sont tous en ligne :

• Javadoc : stats.OnlineNormalEstimator
• Fuente: stats.OnlineNormalEstimator.java
• JUnit Source : test.unit.stats.OnlineNormalEstimatorTest.java
• Page d'accueil de LingPipe

Answer 2

La réponse de base est d'accumuler la somme des deux éléments suivants x (appelé 'sum_x1') et x 2 (appelez-la 'sum_x2') au fur et à mesure. La valeur de l'écart-type est alors :

stdev = sqrt((sum_x2 / n) - (mean * mean)) 

où

mean = sum_x / n

Il s'agit de l'écart-type de l'échantillon ; vous obtenez l'écart-type de la population en utilisant 'n' au lieu de 'n - 1' comme diviseur.

Vous pouvez avoir besoin de vous inquiéter de la stabilité numérique de la différence entre deux grands nombres si vous avez affaire à de grands échantillons. Consultez les références externes dans les autres réponses (Wikipedia, etc.) pour plus d'informations.

Answer 3

Ce n'est peut-être pas ce que vous demandiez, mais ... Si vous utilisez un NumPy array, il fera le travail pour vous, efficacement :

from numpy import array

nums = array(((0.01, 0.01, 0.02, 0.04, 0.03),
              (0.00, 0.02, 0.02, 0.03, 0.02),
              (0.01, 0.02, 0.02, 0.03, 0.02),
              (0.01, 0.00, 0.01, 0.05, 0.03)))

print nums.std(axis=1)
# [ 0.0116619   0.00979796  0.00632456  0.01788854]

print nums.mean(axis=1)
# [ 0.022  0.018  0.02   0.02 ]

À propos, cet article de blog et les commentaires sur les méthodes de calcul des moyennes et des variances en une seule passe suscitent une discussion intéressante :

Calcul de la moyenne et de la variance de l'échantillon en ligne en une seule passe

Answer 4

Statistiques::Descriptives est un module Perl très décent pour ces types de calculs :

#!/usr/bin/perl

use strict; use warnings;

use Statistics::Descriptive qw( :all );

my $data = [
    [ 0.01, 0.01, 0.02, 0.04, 0.03 ],
    [ 0.00, 0.02, 0.02, 0.03, 0.02 ],
    [ 0.01, 0.02, 0.02, 0.03, 0.02 ],
    [ 0.01, 0.00, 0.01, 0.05, 0.03 ],
];

my $stat = Statistics::Descriptive::Full->new;
# You also have the option of using sparse data structures

for my $ref ( @$data ) {
    $stat->add_data( @$ref );
    printf "Running mean: %f\n", $stat->mean;
    printf "Running stdev: %f\n", $stat->standard_deviation;
}
__END__

Sortie :

Running mean: 0.022000
Running stdev: 0.013038
Running mean: 0.020000
Running stdev: 0.011547
Running mean: 0.020000
Running stdev: 0.010000
Running mean: 0.020000
Running stdev: 0.012566

Answer 5

Jetez un coup d'œil à PDL (prononcé "piddle !").

Il s'agit du langage de données Perl, conçu pour les mathématiques de haute précision et le calcul scientifique.

Voici un exemple utilisant vos chiffres....

use strict;
use warnings;
use PDL;

my $figs = pdl [
    [0.01, 0.01, 0.02, 0.04, 0.03],
    [0.00, 0.02, 0.02, 0.03, 0.02],
    [0.01, 0.02, 0.02, 0.03, 0.02],
    [0.01, 0.00, 0.01, 0.05, 0.03],
];

my ( $mean, $prms, $median, $min, $max, $adev, $rms ) = statsover( $figs );

say "Mean scores:     ", $mean;
say "Std dev? (adev): ", $adev;
say "Std dev? (prms): ", $prms;
say "Std dev? (rms):  ", $rms;

Ce qui produit :

Mean scores:     [0.022 0.018 0.02 0.02]
Std dev? (adev): [0.0104 0.0072 0.004 0.016]
Std dev? (prms): [0.013038405 0.010954451 0.0070710678 0.02]
Std dev? (rms):  [0.011661904 0.009797959 0.0063245553 0.017888544]

Jetez un coup d'œil à PDL::Primitive pour plus d'informations sur le statsover fonction. Cela semble suggérer que ADEV est l'"écart type".

Cependant, il peut s'agir de PRMS (ce que montre l'exemple de Statistics::Descriptive de Sinan) ou de RMS (que L'exemple NumPy d'ars spectacles). Je suppose qu'un de ces trois doit avoir raison ;-)

Pour plus d'informations sur les PDL, consultez le site :

• pdl.perl.org (page officielle PDL).
• Guide de référence rapide PDL sur PerlMonks
• L'article du Dr Dobb sur le PDL
• PDL Wiki
• Entrée de Wikipédia pour PDL
• Page du projet SourceForge pour PDL