Comment calculer efficacement un écart type courant?

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)

Ce que je voudrais faire, c'est calculer efficacement la moyenne et l'écart type à chaque index d'une liste, sur tous les éléments du tableau.

Pour faire le moyen, j'ai parcouru le tableau et additionné la valeur à un index donné d'une liste. À 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 faire l'écart type, je boucle à nouveau, maintenant que j'ai la moyenne calculée.

Je voudrais éviter de parcourir le tableau deux fois, une fois pour la moyenne et une fois pour la SD (après avoir 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 pseudo-code convient.

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:

• Wikipédia: algorithmes de calcul de la variance

Il est plus stable numériquement que les collecteurs de somme de carrés simples à deux passes ou en ligne suggérés dans d'autres réponses. La stabilité n'a vraiment d'importance que lorsque vous avez beaucoup de valeurs qui sont proches les unes des autres car elles conduisent à ce que l'on appelle «l' annulation catastrophique » dans la littérature en virgule flottante.

Vous voudrez peut-être également revenir sur 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 à partir de l'échantillon, tandis que la division par N en moyenne sous-estime la variance (car elle ne prend pas en compte la variance entre la moyenne de l'échantillon et la moyenne réelle).

J'ai écrit deux entrées de blog sur le sujet qui vont dans plus de détails, y compris comment supprimer les valeurs précédentes en ligne:

• Calcul de la moyenne et de la variance de l'échantillon en ligne en un seul passage
• Suppression de valeurs dans l'algorithme de Welford pour la moyenne et la variance en ligne

Vous pouvez également jeter un oeil à mon implémentation Java; les tests javadoc, source et unitaires sont tous en ligne:

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

La réponse de base est d'accumuler la somme de x (appelez-la «somme_x1») et de x ² (appelez-la «somme_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 devrez peut-être vous soucier de la stabilité numérique de la différence entre deux grands nombres si vous avez affaire à de grands échantillons. Accédez aux références externes dans d'autres réponses (Wikipédia, etc.) pour plus d'informations.

Voici une traduction littérale pure Python de l'implémentation de l'algorithme de Welford à partir de http://www.johndcook.com/standard_deviation.html :

https://github.com/liyanage/python-modules/blob/master/running_stats.py

import math

class RunningStats:

    def __init__(self):
        self.n = 0
        self.old_m = 0
        self.new_m = 0
        self.old_s = 0
        self.new_s = 0

    def clear(self):
        self.n = 0

    def push(self, x):
        self.n += 1

        if self.n == 1:
            self.old_m = self.new_m = x
            self.old_s = 0
        else:
            self.new_m = self.old_m + (x - self.old_m) / self.n
            self.new_s = self.old_s + (x - self.old_m) * (x - self.new_m)

            self.old_m = self.new_m
            self.old_s = self.new_s

    def mean(self):
        return self.new_m if self.n else 0.0

    def variance(self):
        return self.new_s / (self.n - 1) if self.n > 1 else 0.0

    def standard_deviation(self):
        return math.sqrt(self.variance())

Usage:

rs = RunningStats()
rs.push(17.0)
rs.push(19.0)
rs.push(24.0)

mean = rs.mean()
variance = rs.variance()
stdev = rs.standard_deviation()

print(f'Mean: {mean}, Variance: {variance}, Std. Dev.: {stdev}')

Peut-être pas ce que vous demandiez, mais ... Si vous utilisez un tableau numpy, 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, il y a une discussion intéressante dans cet article de blog et des commentaires sur les méthodes en un seul passage pour calculer les moyennes et les variances:

• http://lingpipe-blog.com/2009/03/19/computing-sample-mean-variance-online-one-pass/

Le module runstats Python est juste pour ce genre de chose. Installez runstats depuis PyPI:

pip install runstats

Les résumés Runstats peuvent produire la moyenne, la variance, l'écart type, l'asymétrie et l'aplatissement en un seul passage de données. Nous pouvons l'utiliser pour créer votre version "en cours d'exécution".

from runstats import Statistics

stats = [Statistics() for num in range(len(data[0]))]

for row in data:

    for index, val in enumerate(row):
        stats[index].push(val)

    for index, stat in enumerate(stats):
        print 'Index', index, 'mean:', stat.mean()
        print 'Index', index, 'standard deviation:', stat.stddev()

Les résumés statistiques sont basés sur la méthode de Knuth et Welford pour calculer l'écart-type en une seule passe, comme décrit dans Art of Computer Programming, Vol 2, p. 232, 3e édition. L'avantage de ceci est des résultats numériquement stables et précis.

Avertissement: je suis l'auteur du module runstats Python.

Statistics :: Descriptive 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__

Production:

C:\Temp> g
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

Jetez un œ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 œil à PDL :: Primitive pour plus d'informations sur la fonction statsover . Cela semble suggérer que l'ADEV est «l'écart type».

Cependant, c'est peut-être PRMS (que montre l'exemple Statistics :: Descriptive de Sinan) ou RMS (ce que montre l'exemple NumPy d'ars). Je suppose que l'un de ces trois doit avoir raison ;-)

Pour plus d'informations sur PDL, consultez:

• pdl.perl.org (page officielle PDL).
• Guide de référence rapide PDL sur PerlMonks
• Article du Dr Dobb sur PDL
• Wiki PDL
• Entrée Wikipedia pour PDL
• Page de projet Sourceforge pour PDL

Quelle est la taille de votre tableau? À moins qu'il ne soit long de plusieurs milliards d'éléments, ne vous inquiétez pas de le parcourir deux fois. Le code est simple et facilement testé.

Ma préférence serait d'utiliser l' extension maths du tableau numpy pour convertir votre tableau de tableaux en un tableau 2D numpy et obtenir directement l'écart type:

>>> x = [ [ 1, 2, 4, 3, 4, 5 ], [ 3, 4, 5, 6, 7, 8 ] ] * 10
>>> import numpy
>>> a = numpy.array(x)
>>> a.std(axis=0) 
array([ 1. ,  1. ,  0.5,  1.5,  1.5,  1.5])
>>> a.mean(axis=0)
array([ 2. ,  3. ,  4.5,  4.5,  5.5,  6.5])

Si ce n'est pas une option et que vous avez besoin d'une solution pure Python, continuez à lire ...

Si votre baie est

x = [ 
      [ 1, 2, 4, 3, 4, 5 ],
      [ 3, 4, 5, 6, 7, 8 ],
      ....
]

Alors l'écart type est:

d = len(x[0])
n = len(x)
sum_x = [ sum(v[i] for v in x) for i in range(d) ]
sum_x2 = [ sum(v[i]**2 for v in x) for i in range(d) ]
std_dev = [ sqrt((sx2 - sx**2)/N)  for sx, sx2 in zip(sum_x, sum_x2) ]

Si vous êtes déterminé à parcourir votre tableau une seule fois, les sommes en cours peuvent être combinées.

sum_x  = [ 0 ] * d
sum_x2 = [ 0 ] * d
for v in x:
   for i, t in enumerate(v):
   sum_x[i] += t
   sum_x2[i] += t**2

Ce n'est pas aussi élégant que la solution de compréhension de liste ci-dessus.

Vous pouvez consulter l'article de Wikipédia sur l'écart type , en particulier la section sur les méthodes de calcul rapide.

Il y a aussi un article que j'ai trouvé qui utilise Python, vous devriez pouvoir utiliser le code sans trop de changement: Messages subliminaux - Exécution des écarts standard .

Je pense que ce numéro vous aidera. Écart-type

Voici un "one-liner", réparti sur plusieurs lignes, dans un style de programmation fonctionnelle:

def variance(data, opt=0):
    return (lambda (m2, i, _): m2 / (opt + i - 1))(
        reduce(
            lambda (m2, i, avg), x:
            (
                m2 + (x - avg) ** 2 * i / (i + 1),
                i + 1,
                avg + (x - avg) / (i + 1)
            ),
            data,
            (0, 0, 0)))

n=int(raw_input("Enter no. of terms:"))

L=[]

for i in range (1,n+1):

    x=float(raw_input("Enter term:"))

    L.append(x)

sum=0

for i in range(n):

    sum=sum+L[i]

avg=sum/n

sumdev=0

for j in range(n):

    sumdev=sumdev+(L[j]-avg)**2

dev=(sumdev/n)**0.5

print "Standard deviation is", dev

Comme le décrit la réponse suivante: pandas / scipy / numpy fournit-il une fonction d'écart type cumulatif? Le module Python Pandas contient une méthode pour calculer l' écart type courant ou cumulé . Pour cela, vous devrez convertir vos données en un dataframe pandas (ou une série si c'est 1D), mais il existe des fonctions pour cela.

J'aime exprimer la mise à jour de cette façon:

def running_update(x, N, mu, var):
    '''
        @arg x: the current data sample
        @arg N : the number of previous samples
        @arg mu: the mean of the previous samples
        @arg var : the variance over the previous samples
        @retval (N+1, mu', var') -- updated mean, variance and count
    '''
    N = N + 1
    rho = 1.0/N
    d = x - mu
    mu += rho*d
    var += rho*((1-rho)*d**2 - var)
    return (N, mu, var)

pour qu'une fonction en un seul passage ressemble à ceci:

def one_pass(data):
    N = 0
    mu = 0.0
    var = 0.0
    for x in data:
        N = N + 1
        rho = 1.0/N
        d = x - mu
        mu += rho*d
        var += rho*((1-rho)*d**2 - var)
        # could yield here if you want partial results
   return (N, mu, var)

notez qu'il s'agit de calculer la variance de l'échantillon (1 / N), et non l'estimation sans biais de la variance de la population (qui utilise un facteur de normalisation 1 / (N-1)). Contrairement aux autres réponses, la variable, varqui suit la variance courante ne croît pas proportionnellement au nombre d'échantillons. En tout temps, c'est juste la variance de l'ensemble d'échantillons vu jusqu'ici (il n'y a pas de «division finale par n» pour obtenir la variance).

Dans une classe, cela ressemblerait à ceci:

class RunningMeanVar(object):
    def __init__(self):
        self.N = 0
        self.mu = 0.0
        self.var = 0.0
    def push(self, x):
        self.N = self.N + 1
        rho = 1.0/N
        d = x-self.mu
        self.mu += rho*d
        self.var += + rho*((1-rho)*d**2-self.var)
    # reset, accessors etc. can be setup as you see fit

Cela fonctionne également pour les échantillons pondérés:

def running_update(w, x, N, mu, var):
    '''
        @arg w: the weight of the current sample
        @arg x: the current data sample
        @arg mu: the mean of the previous N sample
        @arg var : the variance over the previous N samples
        @arg N : the number of previous samples
        @retval (N+w, mu', var') -- updated mean, variance and count
    '''
    N = N + w
    rho = w/N
    d = x - mu
    mu += rho*d
    var += rho*((1-rho)*d**2 - var)
    return (N, mu, var)

Voici un exemple pratique de la façon dont vous pourriez implémenter un écart-type en cours d'exécution avec python et numpy:

a = np.arange(1, 10)
s = 0
s2 = 0
for i in range(0, len(a)):
    s += a[i]
    s2 += a[i] ** 2 
    n = (i + 1)
    m = s / n
    std = np.sqrt((s2 / n) - (m * m))
    print(std, np.std(a[:i + 1]))

Cela imprimera l'écart-type calculé et un écart-type de contrôle calculé avec numpy:

0.0 0.0
0.5 0.5
0.8164965809277263 0.816496580927726
1.118033988749895 1.118033988749895
1.4142135623730951 1.4142135623730951
1.707825127659933 1.707825127659933
2.0 2.0
2.29128784747792 2.29128784747792
2.5819888974716116 2.581988897471611

J'utilise simplement la formule décrite dans ce fil:

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