Algoritma online untuk simpangan absolut rata-rata dan kumpulan data besar

Saya punya sedikit masalah yang membuat saya panik. Saya harus menulis prosedur untuk proses akuisisi online dari rangkaian waktu multivarian. Pada setiap interval waktu (misalnya 1 detik), saya mendapatkan sampel baru, yang pada dasarnya adalah vektor titik mengambang ukuran N. Operasi yang perlu saya lakukan agak rumit:

1. Untuk setiap sampel baru, saya menghitung persentase untuk sampel itu (dengan menormalkan vektor sehingga elemen akan berjumlah 1).

2. Saya menghitung vektor persentase rata-rata dengan cara yang sama, tetapi menggunakan nilai masa lalu.

3. Untuk setiap nilai masa lalu, saya menghitung deviasi absolut dari vektor persentase yang terkait dengan sampel tersebut dengan vektor persentase rata-rata global yang dihitung pada langkah 2. Dengan cara ini, deviasi absolut selalu berupa angka antara 0 (ketika vektor sama dengan rata-rata). vektor) dan 2 (ketika itu benar-benar berbeda).

4. Dengan menggunakan rata-rata penyimpangan untuk semua sampel sebelumnya, saya menghitung deviasi absolut rata-rata, yang lagi-lagi angka antara 0 dan 2.

5. Saya menggunakan deviasi absolut rata-rata untuk mendeteksi apakah sampel baru kompatibel dengan sampel lain (dengan membandingkan deviasi absolutnya dengan deviasi absolut rata-rata dari seluruh rangkaian yang dihitung pada langkah 4).

Karena setiap kali sampel baru dikumpulkan perubahan rata-rata global (dan juga berarti deviasi absolut juga berubah), adakah cara untuk menghitung nilai ini tanpa memindai seluruh data yang ditetapkan beberapa kali? (satu kali untuk menghitung persentase rata-rata global, dan satu kali untuk mengumpulkan penyimpangan absolut). Ok, saya tahu sangat mudah untuk menghitung rata-rata global tanpa memindai seluruh rangkaian, karena saya hanya perlu menggunakan vektor sementara untuk menyimpan jumlah setiap dimensi, tetapi bagaimana dengan deviasi absolut rata-rata? Perhitungannya termasuk abs()operator, jadi saya perlu mengakses semua data masa lalu!

Terima kasih atas bantuan Anda.

Jika Anda dapat menerima beberapa ketidakakuratan, masalah ini dapat diselesaikan dengan mudah dengan menghitung binning . Yaitu, pilih beberapa angka (katakanlah M = 1000 ), kemudian inisialisasi beberapa bilangan bulat B i , j untuk i = 1 ... M dan j = 1 ... N , di mana N adalah ukuran vektor, sebagai nol. Kemudian ketika Anda melihat k th pengamatan vektor percentual, kenaikan B i , j jika j th unsur vektor ini adalah antara $M$$M = 1000$$B_{i,j}$$i = 1\ldots M$$j = 1\ldots N$$N$$k$$B_{i,j}$$j$ dan i / M , mengulangelemen N dari vektor. (Saya mengasumsikan vektor input Anda non-negatif, sehingga ketika Anda menghitung 'persentase' Anda, vektor-vektor itu berada dalam kisaran [ 0 , 1 ] .)$(i-1)/M$$i/M$$N$$[0,1]$

Kapan saja, Anda dapat memperkirakan vektor rata-rata dari tempat sampah, dan deviasi absolut rata-rata. Setelah mengamati vektor-vektor , elemen j dari rata-rata diperkirakan oleh ˉ X j = $K$$j$danelemen ke-jdari deviasi absolut rata-rata diperkirakan sebesar

sunting : ini adalah kasus spesifik dari pendekatan yang lebih umum di mana Anda sedang membangun estimasi kepadatan empiris. Ini bisa dilakukan dengan polinomial, splines, dll, tetapi pendekatan binning adalah yang paling mudah untuk dijelaskan dan diimplementasikan.

Saya telah menggunakan pendekatan berikut di masa lalu untuk menghitung penyimpangan absolusi dengan cukup efisien (perhatikan, ini adalah pendekatan programmer, bukan ahli statistik, jadi pasti ada trik pintar seperti shabbychef yang mungkin lebih efisien).

PERINGATAN: Ini bukan algoritma online. Itu membutuhkan O(n)memori. Selain itu, ia memiliki kinerja kasus terburuk O(n), untuk kumpulan data seperti [1, -2, 4, -8, 16, -32, ...](yaitu sama dengan perhitungan penuh). [1]

Namun, karena masih berkinerja baik dalam banyak kasus penggunaan, mungkin ada baiknya memposting di sini. Misalnya, untuk menghitung penyimpangan absolut dari 10.000 angka acak antara -100 dan 100 ketika setiap item tiba, algoritma saya membutuhkan waktu kurang dari satu detik, sedangkan perhitungan ulang penuh membutuhkan waktu lebih dari 17 detik (pada mesin saya, akan bervariasi per mesin dan sesuai dengan input data). Anda perlu mempertahankan seluruh vektor dalam memori, yang mungkin menjadi kendala untuk beberapa penggunaan. Garis besar algoritma adalah sebagai berikut:

1. Alih-alih memiliki vektor tunggal untuk menyimpan pengukuran masa lalu, gunakan tiga antrian prioritas yang diurutkan (sesuatu seperti tumpukan min / maks). Ketiga daftar mempartisi input menjadi tiga: item lebih besar dari rata-rata, item kurang dari rata-rata dan item sama dengan rata-rata.
2. (Hampir) setiap kali Anda menambahkan item perubahan berarti, jadi kami perlu melakukan partisi ulang. Yang penting adalah sifat partisi yang diurutkan yang berarti bahwa alih-alih memindai setiap item dalam daftar untuk partisi ulang, kita hanya perlu membaca item-item yang kita pindahkan. Sementara dalam kasus terburuk ini masih akan memerlukan O(n)operasi pindah, untuk banyak kasus penggunaan ini tidak demikian.
3. Dengan menggunakan beberapa pembukuan yang cerdas, kita dapat memastikan bahwa penyimpangan dihitung dengan benar setiap saat, saat partisi ulang dan ketika menambahkan item baru.

Beberapa kode contoh, dengan python, ada di bawah ini. Perhatikan bahwa itu hanya memungkinkan item ditambahkan ke daftar, tidak dihapus. Ini dapat dengan mudah ditambahkan, tetapi pada saat saya menulis ini saya tidak perlu melakukannya. Alih-alih mengimplementasikan antrian prioritas sendiri, saya telah menggunakan daftar sortir dari paket blist Daniel Stutzbach yang sangat baik , yang menggunakan B + Tree secara internal.

Pertimbangkan kode ini dilisensikan di bawah lisensi MIT . Ini belum dioptimalkan atau dipoles secara signifikan, tetapi telah bekerja untuk saya di masa lalu. Versi baru akan tersedia di sini . Beri tahu saya jika Anda memiliki pertanyaan, atau temukan bug.

from blist import sortedlist
import operator

class deviance_list:
    def __init__(self):
        self.mean =  0.0
        self._old_mean = 0.0
        self._sum =  0L
        self._n =  0  #n items
        # items greater than the mean
        self._toplist =  sortedlist()
        # items less than the mean
        self._bottomlist = sortedlist(key = operator.neg)
        # Since all items in the "eq list" have the same value (self.mean) we don't need
        # to maintain an eq list, only a count
        self._eqlistlen = 0

        self._top_deviance =  0
        self._bottom_deviance =  0

    @property
    def absolute_deviance(self):
        return self._top_deviance + self._bottom_deviance

    def append(self,  n):
        # Update summary stats
        self._sum += n
        self._n +=  1
        self._old_mean =  self.mean
        self.mean =  self._sum /  float(self._n)

        # Move existing things around
        going_up = self.mean > self._old_mean
        self._rebalance(going_up)

        # Add new item to appropriate list
        if n >  self.mean:
            self._toplist.add(n)
            self._top_deviance +=  n -  self.mean
        elif n == self.mean: 
            self._eqlistlen += 1
        else:
            self._bottomlist.add(n)
            self._bottom_deviance += self.mean -  n


    def _move_eqs(self,  going_up):
        if going_up:
            self._bottomlist.update([self._old_mean] *  self._eqlistlen)
            self._bottom_deviance += (self.mean - self._old_mean) * self._eqlistlen
            self._eqlistlen = 0
        else:
            self._toplist.update([self._old_mean] *  self._eqlistlen)
            self._top_deviance += (self._old_mean - self.mean) * self._eqlistlen
            self._eqlistlen = 0


    def _rebalance(self, going_up):
        move_count,  eq_move_count = 0, 0
        if going_up:
            # increase the bottom deviance of the items already in the bottomlist
            if self.mean !=  self._old_mean:
                self._bottom_deviance += len(self._bottomlist) *  (self.mean -  self._old_mean)
                self._move_eqs(going_up)


            # transfer items from top to bottom (or eq) list, and change the deviances
            for n in iter(self._toplist):
                if n < self.mean:
                    self._top_deviance -= n -  self._old_mean
                    self._bottom_deviance += (self.mean -  n)
                    # we increment movecount and move them after the list
                    # has finished iterating so we don't modify the list during iteration
                    move_count +=  1
                elif n == self.mean:
                    self._top_deviance -= n -  self._old_mean
                    self._eqlistlen += 1
                    eq_move_count +=  1
                else:
                    break
            for _ in xrange(0,  move_count):
                self._bottomlist.add(self._toplist.pop(0))
            for _ in xrange(0,  eq_move_count):
                self._toplist.pop(0)

            # decrease the top deviance of the items remain in the toplist
            self._top_deviance -= len(self._toplist) *  (self.mean -  self._old_mean)
        else:
            if self.mean !=  self._old_mean:
                self._top_deviance += len(self._toplist) *  (self._old_mean -  self.mean)
                self._move_eqs(going_up)
            for n in iter(self._bottomlist): 
                if n > self.mean:
                    self._bottom_deviance -= self._old_mean -  n
                    self._top_deviance += n -  self.mean
                    move_count += 1
                elif n == self.mean:
                    self._bottom_deviance -= self._old_mean -  n
                    self._eqlistlen += 1
                    eq_move_count +=  1
                else:
                    break
            for _ in xrange(0,  move_count):
                    self._toplist.add(self._bottomlist.pop(0))
            for _ in xrange(0,  eq_move_count):
                self._bottomlist.pop(0)

            # decrease the bottom deviance of the items remain in the bottomlist
            self._bottom_deviance -= len(self._bottomlist) *  (self._old_mean -  self.mean)


if __name__ ==  "__main__":
    import random
    dv =  deviance_list()
    # Test against some random data,  and calculate result manually (nb. slowly) to ensure correctness
    rands = [random.randint(-100,  100) for _ in range(0,  1000)]
    ns = []
    for n in rands: 
        dv.append(n)
        ns.append(n)
        print("added:%4d,  mean:%3.2f,  oldmean:%3.2f,  mean ad:%3.2f" %
              (n, dv.mean,  dv._old_mean,  dv.absolute_deviance / dv.mean))
        assert sum(ns) == dv._sum,  "Sums not equal!"
        assert len(ns) == dv._n,  "Counts not equal!"
        m = sum(ns) / float(len(ns))
        assert m == dv.mean,  "Means not equal!"
        real_abs_dev = sum([abs(m - x) for x in ns])
        # Due to floating point imprecision, we check if the difference between the
        # two ways of calculating the asb. dev. is small rather than checking equality
        assert abs(real_abs_dev - dv.absolute_deviance) < 0.01, (
            "Absolute deviances not equal. Real:%.2f,  calc:%.2f" %  (real_abs_dev,  dv.absolute_deviance))

[1] Jika gejalanya menetap, kunjungi dokter Anda.

$X$$X$$X$$s$$s \sqrt{2 / \pi}$

MAD (x) hanyalah dua perhitungan median bersamaan, yang masing-masing dapat dilakukan secara online melalui algoritma binmedian .

Anda dapat menemukan kertas terkait serta kode C dan FORTRAN online di sini .

(Ini hanya penggunaan trik pintar di atas trik pintar Shabbychef, untuk menghemat memori).

Tambahan:

Ada sejumlah metode multi-pass lama untuk menghitung kuantil. Pendekatan yang populer adalah memelihara / memperbarui reservoir pengamatan berukuran deterministik yang dipilih secara acak dari sungai dan menghitung kuantil secara rekursif (lihat ulasan ini ) pada reservoir ini. Pendekatan ini (dan terkait) digantikan oleh yang diusulkan di atas.

Berikut ini memberikan perkiraan yang tidak akurat, meskipun ketidaktepatan akan tergantung pada distribusi data input. Ini adalah algoritma online, tetapi hanya mendekati penyimpangan absolut. Ini didasarkan pada algoritma yang terkenal untuk menghitung varians online, dijelaskan oleh Welford pada 1960-an. Algoritma-nya, diterjemahkan ke dalam R, terlihat seperti:

M2 <- 0
mean <- 0
n <- 0

var.online <- function(x){
    n <<- n + 1
    diff <- x - mean
    mean <<- mean + diff / n
    M2 <<- M2 + diff * (x - mean)
    variance <- M2 / (n - 1)
    return(variance)
}

Kerjanya sangat mirip dengan fungsi varians bawaan R:

set.seed(2099)
n.testitems <- 1000
n.tests <- 100
differences <- rep(NA, n.tests)
for (i in 1:n.tests){
        # Reset counters
        M2 <- 0
        mean <- 0
        n <- 0

        xs <- rnorm(n.testitems)
        for (j in 1:n.testitems){
                v <- var.online(xs[j])
        }

        differences[i] <- abs(v - var(xs))

}
summary(differences)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.000e+00 2.220e-16 4.996e-16 6.595e-16 9.992e-16 1.887e-15 

Memodifikasi algoritma untuk menghitung penyimpangan absolut hanya melibatkan sqrtpanggilan tambahan . Namun, sqrtmemperkenalkan ketidakakuratan yang tercermin dalam hasil:

absolute.deviance.online <- function(x){
    n <<- n + 1
    diff <- x - mean
    mean <<- mean + diff / n
    a.dev <<- a.dev + sqrt(diff * (x - mean))
    return(a.dev)
}

Kesalahan, dihitung seperti di atas, jauh lebih besar daripada perhitungan varians:

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.005126 0.364600 0.808000 0.958800 1.360000 3.312000 

Namun, tergantung pada kasus penggunaan Anda, besarnya kesalahan ini mungkin dapat diterima.