Jak skutecznie obliczyć bieżące odchylenie standardowe?

87

Mam tablicę list liczb, np:

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

Chciałbym efektywnie obliczyć średnią i odchylenie standardowe dla każdego indeksu listy, we wszystkich elementach tablicy.

Aby określić średnią, przeglądałem tablicę w pętli i sumowałem wartość dla danego indeksu listy. Na koniec każdą wartość z mojej „listy średnich” dzielę przez n(pracuję na populacji, a nie na próbie z populacji).

Aby obliczyć odchylenie standardowe, ponownie przechodzę przez pętlę, gdy mam już obliczoną średnią.

Chciałbym uniknąć dwukrotnego przechodzenia przez tablicę, raz dla średniej i raz dla SD (po tym, jak mam średnią).

Czy istnieje skuteczna metoda obliczania obu wartości, przechodząc przez tablicę tylko raz? Każdy kod w języku interpretowanym (np. Perl lub Python) lub pseudokodzie jest w porządku.

python  perl  statistics 
— Alex Reynolds
źródło
7
Inny język, ale ten sam algorytm: stackoverflow.com/questions/895929/…
— dmckee --- ex-moderator kitten
Dzięki, sprawdzę ten algorytm. Brzmi jak to, czego potrzebuję.
— Alex Reynolds
Dzięki za wskazanie mi właściwej odpowiedzi, dmckee. Chciałbym dać Ci znacznik wyboru „najlepsza odpowiedź”, jeśli chcesz poświęcić chwilę na dodanie swojej odpowiedzi poniżej (jeśli chcesz uzyskać punkty).
— Alex Reynolds
1
Istnieje również kilka przykładów na rosettacode.org/wiki/Standard_Deviation
— glenn jackman
1
Wikipedia ma implementację Pythona en.wikipedia.org/wiki/…
— Hamish Grubijan

Odpowiedzi:

116

Odpowiedzią jest użycie algorytmu Welforda, który jest bardzo jasno zdefiniowany po „naiwnych metodach” w:

Jest bardziej stabilny numerycznie niż prosta suma kolektorów kwadratów z dwoma przebiegami lub online sugerowana w innych odpowiedziach. Stabilność naprawdę ma znaczenie tylko wtedy, gdy masz wiele wartości, które są blisko siebie, ponieważ prowadzą one w literaturze do tak zwanego „ katastroficznego anulowania ”.

Możesz również odświeżyć różnicę między dzieleniem przez liczbę próbek (N) i N-1 w obliczeniu wariancji (odchylenie kwadratowe). Dzielenie przez N-1 prowadzi do nieobciążonego oszacowania wariancji z próby, podczas gdy dzielenie przez N średnio niedoszacowuje wariancję (ponieważ nie uwzględnia wariancji między średnią z próby a średnią prawdziwą).

Napisałem dwa wpisy na blogu na ten temat, które zawierają więcej szczegółów, w tym jak usunąć poprzednie wartości online:

Możesz także rzucić okiem na moje narzędzie Java; testy javadoc, źródła i testy jednostkowe są dostępne online:

— Bob Carpenter
źródło
1
+1, za zajęcie się usuwaniem wartości z algorytmu
— Welforda
3
Dobra odpowiedź, +1 za przypomnienie czytelnikowi różnicy między odchyleniem standardowym populacji a odchyleniem próbnym.
— Assad Ebrahim
Po powrocie do tego pytania po tylu latach chciałem tylko podziękować za poświęcenie czasu na udzielenie świetnej odpowiedzi.
— Alex Reynolds
76

Podstawową odpowiedzią jest skumulowanie sumy obu x (nazwij to „sum_x1”) i x 2 (nazwij to „suma_x2”) na bieżąco. Wartość odchylenia standardowego wynosi zatem:

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

gdzie

mean = sum_x / n

To jest odchylenie standardowe próbki; odchylenie standardowe populacji otrzymujesz, używając „n” zamiast „n - 1” jako dzielnika.

Jeśli masz do czynienia z dużymi próbkami, być może będziesz musiał martwić się o stabilność numeryczną podczas obliczania różnicy między dwiema dużymi liczbami. Aby uzyskać więcej informacji, przejdź do odnośników zewnętrznych w innych odpowiedziach (Wikipedia itp.).

— Jonathan Leffler
źródło
To właśnie zamierzałem zasugerować. To najlepszy i najszybszy sposób, przy założeniu, że błędy precyzji nie stanowią problemu.
— Ray Hidayat,
2
Zdecydowałem się na Algorytm Welforda, ponieważ działa on bardziej niezawodnie przy tym samym obciążeniu obliczeniowym.
— Alex Reynolds
2
To jest uproszczona wersja odpowiedzi i może dawać nierzeczywiste wyniki w zależności od danych wejściowych (tj. Gdy sum_x2 <sum_x1 * sum_x1). Aby zapewnić poprawny rzeczywisty wynik, przejdź do `sd = sqrt (((n * sum_x2) - (sum_x1 * sum_x1)) / (n * (n - 1)))
— Dan Tao
2
@Dan wskazuje na ważny problem - powyższa formuła rozkłada się na x> 1, ponieważ w końcu bierzesz sqrt liczby ujemnej. Podejście Knutha to: sqrt ((sum_x2 / n) - (średnia * średnia)) gdzie średnia = (suma_x / n).
— G__
1
@UriLoya - nie powiedziałeś nic o tym, jak obliczasz wartości. Jeśli jednak używasz intw C do przechowywania sumy kwadratów, napotkasz problemy z przepełnieniem podanych wartości.
— Jonathan Leffler
38

Oto dosłowne tłumaczenie implementacji algorytmu Welforda w czystym Pythonie ze strony 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())

Stosowanie:

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}')
— Marc Liyanage
źródło
9
To powinna być akceptowana odpowiedź, ponieważ jest jedyną poprawną i pokazującą algorytm w odniesieniu do Knutha.
— Johan Lundberg
26

Być może nie to, o co prosiłeś, ale ... Jeśli używasz tablicy numpy, wykona pracę za Ciebie, wydajnie:

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 ]

Nawiasem mówiąc, w tym poście na blogu znajduje się kilka interesujących dyskusji i komentarzy na temat jednoprzebiegowych metod obliczania średnich i wariancji:

— ars
źródło
14

Przez moduł Pythona Runstats jest tylko dla tego rodzaju rzeczy. Zainstaluj runstats z PyPI:

pip install runstats

Podsumowania Runstats mogą generować średnią, wariancję, odchylenie standardowe, skośność i kurtoozę w jednym przebiegu danych. Możemy to wykorzystać do stworzenia Twojej "działającej" wersji.

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

Podsumowania statystyczne są oparte na metodzie Knutha i Welforda do obliczania odchylenia standardowego w jednym przebiegu, zgodnie z opisem w Art of Computer Programming, tom 2, str. 232, wydanie trzecie. Zaletą tego są stabilne numerycznie i dokładne wyniki.

Zastrzeżenie: jestem autorem modułu runstats w języku Python.

— GrantJ
źródło
Niezły moduł. Byłoby interesujące, gdyby nie było Statisticsma .popmetodę więc statystyka walcowanie można również obliczyć.
— Gustavo Bezerra
@GustavoBezerra runstatsnie prowadzi wewnętrznej listy wartości, więc nie jestem pewien, czy to możliwe. Ale prośby o ściągnięcie są mile widziane.
— GrantJ,
8

Statistics :: Descriptive to bardzo przyzwoity moduł Perla do tego typu obliczeń:

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

Wynik:

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
— Sinan Ünür
źródło
8

Spójrz na PDL (wymawiane „piddle!”).

To jest język danych Perl, który jest przeznaczony do obliczeń matematycznych i naukowych o wysokiej precyzji.

Oto przykład wykorzystujący twoje dane ...

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;


Który produkuje:

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]


Zajrzyj do PDL :: Primitive, aby uzyskać więcej informacji na temat funkcji statsover . Wydaje się to sugerować, że ADEV jest „odchyleniem standardowym”.

Jednak może to być PRMS (co pokazuje Sinan's Statistics :: opisowy przykład) lub RMS (który pokazuje przykład z NumPy Arsa). Chyba jeden z tych trzech musi mieć rację ;-)

Aby uzyskać więcej informacji na temat języka PDL, zobacz:

— draegtun
źródło
1
To nie są bieżące obliczenia.
— Jake
3

Jak duża jest twoja tablica? O ile nie ma zilionów elementów, nie martw się o dwukrotne zapętlenie. Kod jest prosty i łatwy do przetestowania.

Wolałbym użyć rozszerzenia numpy array maths do konwersji tablicy tablic na tablicę numpy 2D i bezpośredniego uzyskania odchylenia standardowego:

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

Jeśli to nie jest opcja i potrzebujesz czystego rozwiązania w języku Python, czytaj dalej ...

Jeśli twoja tablica to 

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

Wtedy odchylenie standardowe wynosi:

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

Jeśli jesteś zdecydowany wykonać pętlę przez tablicę tylko raz, sumy bieżące można łączyć.

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

To nie jest tak eleganckie, jak powyższe rozwiązanie do rozumienia listy.

— Stephen Simmons
źródło
W rzeczywistości mam do czynienia z milionami liczb, co motywuje mnie do znalezienia skutecznego rozwiązania. Dzięki!
— Alex Reynolds
nie chodzi o to, jak duży jest zestaw danych, chodzi o to, jak CZĘSTO, muszę wykonać 3500 różnych obliczeń odchylenia standardowego ponad 500 elementów w każdym obliczeniu na sekundę
— PirateApp
1

Oto „jedna linijka”, podzielona na wiele linii, w funkcjonalnym stylu programowania:

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)))
— user541686
źródło
1 
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
— Anuraag
źródło
1

Aktualizację chcę wyrazić w ten sposób:

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)

aby funkcja jednoprzebiegowa wyglądała następująco:

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)

zwróć uwagę, że jest to obliczanie wariancji próby (1 / N), a nie nieobciążonej oceny wariancji populacji (która wykorzystuje współczynnik normalizacji 1 / (N-1)). W odróżnieniu od pozostałych odpowiedzi, zmienna, varczyli śledzenie bieżącej wariancji, nie rośnie proporcjonalnie do liczby próbek. Przez cały czas jest to po prostu wariancja dotychczasowego zbioru próbek (nie ma końcowego „dzielenia przez n” przy uzyskiwaniu wariancji).

Na zajęciach wyglądałoby to tak:

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

Działa to również w przypadku próbek ważonych:

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)
— Dave
źródło
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.