Jak skutecznie obliczyć biegowe odchylenie standardowe?

Podstawową odpowiedzią jest zebranie sumy obu x (nazwij to 'sum_x1') i x² (nazwij to "sum_x2"). Wartość odchylenia standardowego wynosi wtedy:

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

Gdzie

mean = sum_x / n

Jest to odchylenie standardowe próbki; otrzymujesz odchylenie standardowe populacji za pomocą " n "zamiast" n - 1 " jako dzielnika.

Być może będziesz musiał martwić się o stabilność liczbową biorąc różnicę między dwoma dużymi liczbami, jeśli masz do czynienia z duże próbki. Przejdź do zewnętrznych odniesień w innych odpowiedziach (Wikipedia itp.), Aby uzyskać więcej informacji.

Może nie o to prosiłeś, ale ..... Jeśli używasz tablicy numpy, będzie ona wykonywać pracę za Ciebie, efektywnie:

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 wpisie na blogu jest ciekawa dyskusja i komentarze na temat jednokrotnych metod obliczania środków i wariancji:

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

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

Użycie:

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

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

Moduł runstats Pythona służy właśnie do tego typu rzeczy. Zainstaluj runstaty z PyPI:

pip install runstats

Podsumowania Runstats mogą wytworzyć średnią, wariancję, odchylenie standardowe, skośność i kurtozę w jednym przejściu danych. Możemy użyć tego do stworzenia Twojej "uruchomionej" 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 oparte są na metodzie Knutha i Welforda obliczania odchylenia standardowego w jednym przebiegu, opisanej w Art of Computer Programming, Vol 2, p. 232, 3rd edycja. Zaletą tego są liczbowo stabilne i dokładne wyniki.

Disclaimer: jestem autorem modułu runstats Pythona.

Spójrz na PDL (wymawiane "piddle!").

Jest to język danych Perla, który jest przeznaczony do matematyki o wysokiej precyzji i obliczeń naukowych.

Oto przykład z wykorzystaniem Twoich liczb....

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 wytwarza:

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]

Zobacz PDL::Primitive, Aby uzyskać więcej informacji na temat funkcji statsover. Wydaje się to sugerować, że ADEV jest "odchyleniem standardowym".

Jednak może PRMS (które pokazują statystyki Sinana::opisowy przykład) lub RMS (który pokazuje przykład NumPy ARS). Chyba jedna z tych trzech musi mieć rację; -)

Więcej informacji na temat PDL można znaleźć na stronie:

• pdl.perl.org (Oficjalna strona PDL).
• PDL quick reference guide on PerlMonks
• artykuł Dr. Dobba o PDL
• PDL Wiki
• wpis w Wikipedii dla PDL
• strona projektu Sourceforge dla PDL

Statystyki:: opisowy jest bardzo przyzwoitym modułem 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__

Wyjście:

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

Jak duża jest Twoja tablica? Jeśli nie jest to zillions elementów długi, nie martw się o zapętlenie przez niego dwa razy. Kod jest prosty i łatwy do przetestowania.

Moim preferowaniem byłoby użycie rozszerzenia numpy array maths do konwersji tablicy tablic do tablicy numpy 2D i uzyskania odchylenia standardowego bezpośrednio:

>>> 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 Pythona, Czytaj dalej...

Jeśli tablica jest

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

Wtedy standard odchylenie 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 zdecydujesz się na pętlę przez tablicę tylko raz, sumy bieżące mogą być łączone.

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 rozwiązanie do rozumienia listy powyżej.

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

Możesz zajrzeć do artykułu Wikipedii na temat odchylenia standardowego , w szczególności sekcji o szybkich metodach obliczeniowych.

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

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 poniższa odpowiedź: Czy pandy / scipy / numpy zapewniają skumulowaną funkcję odchylenia standardowego? Moduł Python Pandas zawiera metodę obliczania biegowego lub skumulowanego odchylenia standardowego . Do tego będziesz musiał przekonwertować dane na ramkę danych pandy (lub serię, jeśli jest to 1D), ale są do tego funkcje.

Oto "jednolinijkowy", rozłożony na wiele linii, w stylu programowania funkcyjnego:

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

Lubię wyrażać aktualizację 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 tak:

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)

Należy zauważyć, że jest to obliczanie wariancji próbki (1 / N), A nie bezstronne oszacowanie wariancji populacji (która wykorzystuje współczynnik normalizacji 1 / (N-1)). W przeciwieństwie do innych odpowiedzi, zmienna var, czyli śledzenie wariancji biegowej, nie rośnie proporcjonalnie do liczby próbek. Przez cały czas jest to tylko wariancja zestawu próbek widzianych do tej pory (nie ma ostatecznego "dzielenia przez n" w uzyskiwaniu wariancji).

W klasie 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ż dla 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)