Bagaimana cara menghitung deviasi standar yang sedang berjalan secara efisien?

Saya memiliki serangkaian daftar angka, misalnya:

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

Yang ingin saya lakukan adalah menghitung mean dan deviasi standar secara efisien pada setiap indeks daftar, di semua elemen array.

Untuk melakukan maksudnya, saya telah melakukan perulangan melalui array dan menjumlahkan nilai pada indeks tertentu dari sebuah daftar. Pada akhirnya, saya membagi setiap nilai dalam "daftar rata-rata" saya dengan n(Saya bekerja dengan populasi, bukan sampel dari populasi).

Untuk melakukan deviasi standar, saya mengulang lagi, sekarang setelah rata-rata dihitung.

Saya ingin menghindari melalui array dua kali, sekali untuk mean dan kemudian sekali untuk SD (setelah saya punya mean).

Apakah ada metode yang efisien untuk menghitung kedua nilai, hanya melalui larik sekali? Kode apa pun dalam bahasa yang ditafsirkan (misalnya Perl atau Python) atau pseudocode baik-baik saja.

Jawabannya adalah dengan menggunakan algoritme Welford, yang didefinisikan dengan sangat jelas setelah "metode naif" di:

• Wikipedia: Algoritma untuk menghitung varians

Ini lebih stabil secara numerik daripada jumlah kolektor kotak dua lintasan atau online sederhana yang disarankan dalam tanggapan lain. Stabilitas hanya benar-benar penting jika Anda memiliki banyak nilai yang dekat satu sama lain karena mengarah pada apa yang dikenal sebagai " pembatalan bencana " dalam literatur floating point.

Anda mungkin juga ingin mempelajari perbedaan antara membagi dengan jumlah sampel (N) dan N-1 dalam perhitungan varians (deviasi kuadrat). Membagi dengan N-1 mengarah pada estimasi varians yang tidak bias dari sampel, sedangkan membagi dengan N rata-rata meremehkan varians (karena tidak memperhitungkan varians antara mean sampel dan mean sebenarnya).

Saya menulis dua entri blog tentang topik yang menjelaskan lebih detail, termasuk cara menghapus nilai sebelumnya secara online:

• Menghitung Rata-Rata Sampel dan Varians Online dalam Satu Langkah
• Menghapus Nilai dalam Algoritma Welford untuk Rata-rata dan Varians Online

Anda juga dapat melihat implement Java saya; tes javadoc, sumber, dan unit semuanya online:

• Javadoc: stats.OnlineNormalEstimator
• Sumber: stats.OnlineNormalEstimator.java
• Sumber JUnit: test.unit.stats.OnlineNormalEstimatorTest.java
• Halaman Beranda LingPipe

Jawaban dasarnya adalah mengakumulasikan kedua x (sebut saja 'sum_x1') dan x ² (sebut saja 'sum_x2') saat Anda melanjutkan. Nilai deviasi standar kemudian:

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

dimana

mean = sum_x / n

Ini adalah deviasi standar sampel; Anda mendapatkan simpangan baku populasi menggunakan 'n' alih-alih 'n - 1' sebagai pembagi.

Anda mungkin perlu mengkhawatirkan stabilitas numerik dalam mengambil selisih antara dua bilangan besar jika Anda berurusan dengan sampel besar. Buka referensi eksternal di jawaban lain (Wikipedia, dll) untuk informasi lebih lanjut.

Berikut adalah terjemahan Python murni literal dari implementasi algoritma Welford dari 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())

Pemakaian:

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

Mungkin bukan yang Anda tanyakan, tapi ... Jika Anda menggunakan array numpy, ini akan bekerja untuk Anda, secara efisien:

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 ]

Ngomong-ngomong, ada beberapa diskusi menarik di posting blog ini dan komentar tentang metode sekali jalan untuk menghitung sarana dan varians:

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

The Python RUNSTATS Modul ini hanya hal semacam ini. Instal runstats dari PyPI:

pip install runstats

Ringkasan runstats dapat menghasilkan mean, variance, standar deviasi, skewness, dan kurtosis dalam satu kali lintasan data. Kami dapat menggunakan ini untuk membuat versi "berjalan" Anda.

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

Ringkasan statistik didasarkan pada metode Knuth dan Welford untuk menghitung deviasi standar dalam satu kali jalan seperti yang dijelaskan dalam Art of Computer Programming, Vol 2, hal. 232, edisi ke-3. Manfaatnya adalah hasil yang stabil dan akurat secara numerik.

Penafian: Saya adalah pembuat modul runstats Python.

Statistics :: Descriptive adalah modul Perl yang sangat layak untuk jenis kalkulasi berikut:

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

Keluaran:

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

Coba lihat PDL (dibaca "piddle!").

Ini adalah Bahasa Data Perl yang dirancang untuk matematika presisi tinggi dan komputasi ilmiah.

Berikut ini contoh penggunaan gambar Anda ....

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;

Yang menghasilkan:

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]

Lihat PDL :: Primitive untuk informasi lebih lanjut tentang fungsi statsover . Ini sepertinya memberi kesan bahwa ADEV adalah "deviasi standar".

Namun mungkin PRMS (yang ditunjukkan oleh Sinan's Statistics :: Contoh Deskriptif) atau RMS (yang ditunjukkan oleh contoh NumPy ars). Saya kira salah satu dari ketiganya pasti benar ;-)

Untuk informasi PDL selengkapnya, lihat:

• pdl.perl.org (halaman resmi PDL).
• Panduan referensi cepat PDL di PerlMonks
• Artikel Dr. Dobb tentang PDL
• PDL Wiki
• Entri Wikipedia untuk PDL
• Halaman proyek Sourceforge untuk PDL

Seberapa besar array Anda? Kecuali jika panjangnya ziliunan elemen, jangan khawatir tentang mengulanginya dua kali. Kode ini sederhana dan mudah diuji.

Preferensi saya adalah menggunakan ekstensi matematika numpy array untuk mengubah array array Anda menjadi array 2D numpy dan mendapatkan deviasi standar secara langsung:

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

Jika itu bukan pilihan dan Anda membutuhkan solusi Python murni, teruslah membaca ...

Jika array Anda adalah

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

Maka standar deviasinya adalah:

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

Jika Anda ditentukan untuk mengulang melalui array Anda hanya sekali, jumlah yang berjalan dapat digabungkan.

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

Ini tidak seanggun solusi pemahaman daftar di atas.

Anda dapat melihat artikel Wikipedia tentang Standar Deviasi , khususnya bagian tentang metode penghitungan Cepat.

Ada juga artikel yang saya temukan yang menggunakan Python, Anda harus dapat menggunakan kode di dalamnya tanpa banyak perubahan: Subliminal Messages - Running Standard Deviations .

Saya rasa masalah ini akan membantu Anda. Simpangan baku

Berikut adalah "satu baris", tersebar di beberapa baris, dalam gaya pemrograman fungsional:

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

Seperti yang dijelaskan oleh jawaban berikut: Apakah pandas / scipy / numpy menyediakan fungsi deviasi standar kumulatif? Modul Python Pandas berisi metode untuk menghitung deviasi standar kumulatif atau berjalan . Untuk itu Anda harus mengubah data Anda menjadi pandas dataframe (atau seri jika 1D), tetapi ada fungsi untuk itu.

Saya ingin mengungkapkan pembaruan dengan cara ini:

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)

sehingga fungsi one-pass akan terlihat seperti ini:

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)

perhatikan bahwa ini menghitung varians sampel (1 / N), bukan estimasi tidak bias dari varians populasi (yang menggunakan faktor normalisasi 1 / (N-1)). Berbeda dengan jawaban lainnya, variabel,, varyaitu pelacakan running varians tidak tumbuh sebanding dengan jumlah sampel. Sepanjang waktu itu hanya varians dari himpunan sampel yang dilihat sejauh ini (tidak ada "pembagian dengan n" akhir untuk mendapatkan varians).

Di kelas akan terlihat seperti ini:

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

Ini juga berfungsi untuk sampel berbobot:

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)

Berikut adalah contoh praktis bagaimana Anda dapat mengimplementasikan deviasi standar yang sedang berjalan dengan python dan numpy:

a = np.arange(1, 10)
s = 0
s2 = 0
for i in range(0, len(a)):
    s += a[i]
    s2 += a[i] ** 2 
    n = (i + 1)
    m = s / n
    std = np.sqrt((s2 / n) - (m * m))
    print(std, np.std(a[:i + 1]))

Ini akan mencetak deviasi standar yang dihitung dan standar deviasi cek dihitung dengan numpy:

0.0 0.0
0.5 0.5
0.8164965809277263 0.816496580927726
1.118033988749895 1.118033988749895
1.4142135623730951 1.4142135623730951
1.707825127659933 1.707825127659933
2.0 2.0
2.29128784747792 2.29128784747792
2.5819888974716116 2.581988897471611

Saya hanya menggunakan rumus yang dijelaskan di utas ini:

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