Wie kann eine laufende Standardabweichung effizient berechnet werden?

Ich habe eine Reihe von Listen mit Zahlen, z.

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

Ich möchte den Mittelwert und die Standardabweichung an jedem Index einer Liste über alle Array-Elemente hinweg effizient berechnen.

Um den Mittelwert zu ermitteln, habe ich das Array durchlaufen und den Wert an einem bestimmten Index einer Liste summiert. Am Ende teile ich jeden Wert in meiner "Durchschnittsliste" durch n(ich arbeite mit einer Population, nicht mit einer Stichprobe aus der Population).

Um die Standardabweichung zu machen, durchlaufe ich erneut, nachdem ich den Mittelwert berechnet habe.

Ich möchte vermeiden, das Array zweimal durchzugehen, einmal für den Mittelwert und dann einmal für die SD (nachdem ich einen Mittelwert habe).

Gibt es eine effiziente Methode zur Berechnung beider Werte, bei der das Array nur einmal durchlaufen wird? Jeder Code in einer interpretierten Sprache (z. B. Perl oder Python) oder Pseudocode ist in Ordnung.

Die Antwort ist die Verwendung des Welford-Algorithmus, der nach den "naiven Methoden" in:

• Wikipedia: Algorithmen zur Varianzberechnung

Es ist numerisch stabiler als die in anderen Antworten vorgeschlagene einfache Summe von Quadratsammlern mit zwei Durchgängen oder online. Die Stabilität ist nur dann wirklich wichtig, wenn Sie viele Werte haben, die nahe beieinander liegen, da sie in der Gleitkomma-Literatur zu einer sogenannten " katastrophalen Löschung " führen.

Möglicherweise möchten Sie auch den Unterschied zwischen der Division durch die Anzahl der Stichproben (N) und N-1 in der Varianzberechnung (quadratische Abweichung) auffrischen. Die Division durch N-1 führt zu einer unvoreingenommenen Schätzung der Varianz aus der Stichprobe, während die Division durch N im Durchschnitt die Varianz unterschätzt (da die Varianz zwischen dem Stichprobenmittelwert und dem wahren Mittelwert nicht berücksichtigt wird).

Ich habe zwei Blogeinträge zu diesem Thema geschrieben, in denen weitere Details behandelt werden, darunter das Löschen früherer Werte online:

• Online-Berechnung von Stichprobenmittelwert und Varianz in einem Durchgang
• Löschen von Werten in Welfords Algorithmus für Online-Mittelwert und Varianz

Sie können sich auch mein Java-Gerät ansehen. Die Javadoc-, Quell- und Unit-Tests sind alle online:

• Javadoc: stats.OnlineNormalEstimator
• Quelle: stats.OnlineNormalEstimator.java
• JUnit Quelle: test.unit.stats.OnlineNormalEstimatorTest.java
• LingPipe-Startseite

Die grundlegende Antwort besteht darin, die Summe von x (nennen Sie es 'sum_x1') und x ² (nennen Sie es 'sum_x2') zu akkumulieren, während Sie gehen. Der Wert der Standardabweichung ist dann:

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

wo

mean = sum_x / n

Dies ist die Standardabweichung der Stichprobe. Sie erhalten die Populationsstandardabweichung mit 'n' anstelle von 'n - 1' als Divisor.

Möglicherweise müssen Sie sich über die numerische Stabilität der Differenz zwischen zwei großen Zahlen Gedanken machen, wenn Sie mit großen Stichproben arbeiten. Weitere Informationen finden Sie in den externen Referenzen in anderen Antworten (Wikipedia usw.).

Hier ist eine wörtliche reine Python-Übersetzung der Implementierung des Welford-Algorithmus von 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())

Verwendung:

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

Vielleicht nicht das, was Sie gefragt haben, aber ... Wenn Sie ein numpy-Array verwenden, erledigt es die Arbeit effizient für Sie:

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 ]

Übrigens gibt es in diesem Blog-Beitrag einige interessante Diskussionen und Kommentare zu One-Pass-Methoden zur Berechnung von Mitteln und Abweichungen:

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

Das Python-Runstats-Modul ist genau für diese Art von Dingen gedacht . Installieren Sie Runstats von PyPI:

pip install runstats

Runstats-Zusammenfassungen können den Mittelwert, die Varianz, die Standardabweichung, die Schiefe und die Kurtosis in einem einzigen Datenlauf erzeugen. Wir können dies verwenden, um Ihre "laufende" Version zu erstellen.

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

Statistikzusammenfassungen basieren auf der Knuth- und Welford-Methode zur Berechnung der Standardabweichung in einem Durchgang, wie in Art of Computer Programming, Vol. 2, p. 232, 3. Auflage. Der Vorteil davon sind numerisch stabile und genaue Ergebnisse.

Haftungsausschluss: Ich bin der Autor des Python-Runstats-Moduls.

Statistics :: Descriptive ist ein sehr anständiges Perl-Modul für diese Art von Berechnungen:

#!/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__

Ausgabe:

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

Schauen Sie sich PDL an (ausgesprochen "Piddle!").

Dies ist die Perl-Datensprache, die für hochpräzise Mathematik und wissenschaftliches Rechnen entwickelt wurde.

Hier ist ein Beispiel mit Ihren Zahlen ....

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;

Welches produziert:

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]

Werfen Sie einen Blick auf PDL :: Primitive für weitere Informationen über die statsover Funktion. Dies scheint darauf hinzudeuten, dass ADEV die "Standardabweichung" ist.

Es kann sich jedoch um PRMS (das Sinans Statistics :: Descriptive-Beispiel zeigt) oder RMS (das NsPy-Beispiel von ars zeigt) handeln. Ich denke einer dieser drei muss stimmen ;-)

Weitere PDL-Informationen finden Sie unter:

• pdl.perl.org (offizielle PDL-Seite).
• PDL-Kurzanleitung zu PerlMonks
• Dr. Dobbs Artikel über PDL
• PDL Wiki
• Wikipedia-Eintrag für PDL
• Sourceforge-Projektseite für PDL

Wie groß ist dein Array? Machen Sie sich keine Sorgen, wenn Sie nicht zig Millionen Elemente lang sind. Der Code ist einfach und leicht zu testen.

Ich würde es vorziehen, die mathematische Erweiterung des numpy- Arrays zu verwenden, um Ihr Array von Arrays in ein numpy-2D-Array zu konvertieren und die Standardabweichung direkt zu erhalten:

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

Wenn dies keine Option ist und Sie eine reine Python-Lösung benötigen, lesen Sie weiter ...

Wenn Ihr Array ist

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

Dann ist die Standardabweichung:

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

Wenn Sie entschlossen sind, Ihr Array nur einmal zu durchlaufen, können die laufenden Summen kombiniert werden.

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

Dies ist bei weitem nicht so elegant wie die oben beschriebene Lösung zum Listenverständnis.

Sie können sich den Wikipedia-Artikel über Standardabweichung ansehen , insbesondere den Abschnitt über schnelle Berechnungsmethoden.

Ich habe auch einen Artikel gefunden, der Python verwendet. Sie sollten in der Lage sein, den darin enthaltenen Code ohne große Änderungen zu verwenden: Unterschwellige Nachrichten - Ausführen von Standardabweichungen .

Ich denke, dieses Problem wird Ihnen helfen. Standardabweichung

Hier ist ein "Einzeiler", der über mehrere Zeilen verteilt ist und einen funktionalen Programmierstil aufweist:

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

Wie die folgende Antwort beschreibt: Bietet Pandas / Scipy / Numpy eine kumulative Standardabweichungsfunktion? Das Python Pandas-Modul enthält eine Methode zur Berechnung der laufenden oder kumulativen Standardabweichung . Dafür müssen Sie Ihre Daten in einen Pandas-Datenrahmen (oder eine Serie, wenn es sich um 1D handelt) konvertieren, aber dafür gibt es Funktionen.

Ich möchte das Update folgendermaßen ausdrücken:

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)

so dass eine One-Pass-Funktion folgendermaßen aussehen würde:

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)

Beachten Sie, dass dies die Stichprobenvarianz (1 / N) berechnet, nicht die unvoreingenommene Schätzung der Populationsvarianz (die einen 1 / (N-1) -Normalisierungsfaktor verwendet). Im Gegensatz zu den anderen Antworten varwächst die Variable, die die laufende Varianz verfolgt, nicht proportional zur Anzahl der Stichproben. Zu allen Zeiten ist es nur die Varianz des bisher gesehenen Satzes von Stichproben (es gibt keine endgültige "Division durch n", um die Varianz zu erhalten).

In einer Klasse würde es so aussehen:

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

Dies funktioniert auch für gewichtete Proben:

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)