Comment calculer efficacement un écart type de fonctionnement?

la réponse de base est d'accumuler la somme de x (appelez-la 'sum_x1') et x ² (appelez-le' 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

c'est 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 vous devez vous inquiéter de la stabilité numérique de prendre la différence entre deux grands nombres si vous avez affaire à de grands échantillons. Allez aux références externes dans les autres réponses (Wikipedia, etc) pour plus d'informations.

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 ]

soit dit en passant, il y a une discussion intéressante dans ce billet de blog et des commentaires sur les méthodes mono-pass pour calculer les moyens et les variances:

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

Voici une traduction littérale en Python pur de L'implémentation de L'algorithme de Welford de http://www.johndcook.com/standard_deviation.html :

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

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())

Utilisation:

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

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

le Python runstats Module est pour ce genre de chose. Installer "runstats" de PyPI:

pip install runstats

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

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 Knuth et Welford pour le calcul standard déviation en une seule passe comme décrit dans l'Art de la Programmation Informatique, Vol 2, p. 232, 3e édition. L'avantage de ceci est numériquement stable et des résultats exacts.

clause de non-responsabilité: je suis l'auteur du module Python runstats.

regarder PDL (prononcé "piddle!").

C'est le langage de données Perl qui est conçu pour les mathématiques de haute précision et l'informatique scientifique.

voici un exemple en 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;

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]

Regardez PDL:: primitif pour plus d'informations sur la fonction statsover . Cela semble indiquer que L'ADEV est l ' "écart-type".

cependant, il peut PRMS (ce que les statistiques de Sinan::exemple descriptif montrent) ou RMS (qui ars num Py exemple montre). Je suppose que l'une de ces trois doit être droit ;-)

pour plus D'informations sur PDL, consultez:

• pdl.perl.org (page officielle de la PDL).
• PDL guide de référence rapide sur PerlMonks
• article du Dr Dobb sur PDL
• PDL Wiki
• entrée de Wikipedia pour les PDL
• page du projet sur Sourceforge pour les PDL

statistiques:: 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__

sortie:

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

Quelle est la taille de votre réseau? À moins que ce ne soit des millions d'éléments de long, ne vous inquiétez pas de le traverser en boucle deux fois. Le code est simple et facile à tester.

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

>>> 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 Python pure, continuez à lire...

si votre tableau 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é à boucler votre tableau une seule fois, les sommes en cours d'exécution 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.

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

vous pouvez consulter L'article de Wikipedia sur é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 y utiliser le code sans grand changement: Subliminal Messages - Running Standard Deviations .

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 la réponse suivante décrit: Ne pandas/scipy/numpy fournir une norme cumulative fonction d'écart? Le module Python Pandas contient une méthode pour calculer l'exécution ou écart-type cumulatif . Pour cela, vous devez convertir vos données en une base de données pandas (ou une série si elle est 1D), mais il y a des fonctions pour cela.

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)))

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)

de sorte qu'une fonction à un seul passage ressemblerait à 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 du calcul de la variance de l'échantillon (1/N), et non de l'estimation non biaisée de la variance de la population (qui utilise un facteur de normalisation 1/(N-1)). Contrairement aux autres réponses, la variable var , qui suit la variance courante , ne croît pas en proportion du nombre d'échantillons. À tous temps il s'agit simplement de la variance de l'ensemble des échantillons vus jusqu'à présent (il n'y a pas de "division par n" finale dans l'obtention de la variance).

dans une classe il 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 aussi 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)