Jak skutecznie obliczyć bieżące odchylenie standardowe?

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.

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

• Wikipedia: Algorytmy obliczania wariancji

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:

• Obliczanie średniej i wariancji próbki online w jednym przebiegu
• Usuwanie wartości w algorytmie Welforda dla średniej i wariancji online

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

• Javadoc: stats.OnlineNormalEstimator
• Źródło: stats.OnlineNormalEstimator.java
• Źródło JUnit: test.unit.stats.OnlineNormalEstimatorTest.java
• Strona główna LingPipe

Podstawową odpowiedzią jest skumulowanie sumy obu x (nazwij to „sum_x1”) i x ² (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.).

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

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:

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

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.

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

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:

• pdl.perl.org (oficjalna strona PDL).
• Skrócona instrukcja obsługi PDL w PerlMonks
• Artykuł dr Dobba na temat PDL
• PDL Wiki
• Wpis w Wikipedii dotyczący PDL
• Strona projektu Sourceforge dla PDL

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.

Możesz zajrzeć do artykułu Wikipedii o odchyleniu standardowym , w szczególności do sekcji o szybkich metodach obliczeniowych.

Jest też artykuł, który znalazłem, który używa Pythona, powinieneś być w stanie użyć w nim kodu bez większych zmian: Wiadomości podprogowe - uruchamianie odchyleń standardowych .

Myślę, że ten problem ci pomoże. Odchylenie standardowe

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

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

Jak opisuje następująca odpowiedź: Czy pandy / scipy / numpy zapewniają skumulowaną funkcję odchylenia standardowego? Moduł Python Pandas zawiera metodę obliczania bieżącego lub skumulowanego odchylenia standardowego . W tym celu będziesz musiał przekonwertować dane na ramkę danych pandy (lub serię, jeśli jest 1D), ale są do tego funkcje.

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)