Algorytmy „on-line” (iteracyjne) do estymacji mediany statystycznej, modu, skośności, kurtozy?

Czy istnieje algorytm do szacowania mediany, trybu, skośności i / lub kurtozy zbioru wartości, ale NIE wymaga to jednoczesnego przechowywania wszystkich wartości w pamięci?

Chciałbym obliczyć podstawowe statystyki:

• mean: średnia arytmetyczna
• wariancja: średnia kwadratów odchyleń od średniej
• odchylenie standardowe: pierwiastek kwadratowy z wariancji
• mediana: wartość, która oddziela większą połowę liczb od mniejszej połowy
• tryb: najczęstsza wartość znaleziona w zestawie
• skośność: tl; dr
• kurtosis: tl; dr

Podstawowymi formułami do obliczania któregokolwiek z nich są arytmetyka podstawowa i znam je. Istnieje również wiele bibliotek statystyk, które je implementują.

Moim problemem jest duża liczba (miliardy) wartości w zestawach, które obsługuję: pracując w Pythonie, nie mogę po prostu sporządzić listy lub mieszania z miliardami elementów. Nawet jeśli napisałem to w C, tablice zawierające miliardy elementów nie są zbyt praktyczne.

Dane nie są posortowane. Jest wytwarzany losowo, w locie, przez inne procesy. Rozmiar każdego zestawu jest bardzo zmienny, a rozmiary nie będą znane z góry.

Dowiedziałem się już, jak całkiem dobrze radzić sobie ze średnią i wariancją, iterując po każdej wartości w zestawie w dowolnej kolejności. (Właściwie w moim przypadku biorę je w kolejności, w jakiej są generowane). Oto algorytm, którego używam, dzięki uprzejmości http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#On-line_algorithm :

• Zainicjuj trzy zmienne: count, sum i sum_of_squares
• Dla każdej wartości:
  • Liczba przyrostów.
  • Dodaj wartość do zsumowania.
  • Dodaj kwadrat wartości do sum_of_squares.
• Podzielić sumę przez liczbę, przechowując jako zmienną średnią.
• Podzielić sumę_kwadratów przez liczbę, przechowując jako zmienną średnią_kwadratów.
• Średnia kwadratowa, przechowywana jako kwadrat_średniej.
• Odejmij square_of_mean od mean_of_squares, przechowując jako wariancję.
• Średnia wyjściowa i wariancja.

Ten algorytm „on-line” ma słabe punkty (np. Problemy z dokładnością, ponieważ sum_of_squares szybko rośnie niż zakres liczb całkowitych lub precyzja typu float), ale zasadniczo daje mi to, czego potrzebuję, bez konieczności przechowywania każdej wartości w każdym zestawie.

Ale nie wiem, czy istnieją podobne techniki szacowania dodatkowych statystyk (mediana, mod, skośność, kurtoza). Mógłbym żyć z obciążonym estymatorem lub nawet metodą, która do pewnego stopnia ogranicza dokładność, o ile pamięć wymagana do przetwarzania wartości N jest znacznie mniejsza niż O (N).

Wskazanie mi istniejącej biblioteki statystyk również pomoże, jeśli biblioteka ma funkcje obliczania jednej lub więcej z tych operacji „on-line”.

Skośność i kurtozy

Dla algorytmów on-line skośności i kurtozy (wzdłuż linii wariancji), zobacz na tej samej stronie wiki tutaj równoległych algorytmów statystycznych wyższej chwili.

Mediana

Mediana jest trudna bez posortowanych danych. Jeśli wiesz, ile masz punktów danych, w teorii wystarczy posortować je tylko częściowo, np. Za pomocą algorytmu selekcji . Jednak to nie pomaga zbytnio przy miliardach wartości. Sugerowałbym użycie liczników częstotliwości, zobacz następną sekcję.

Mediana i tryb z licznikami częstotliwości

Jeśli są to liczby całkowite, policzyłbym częstotliwości , prawdopodobnie odcinając najwyższe i najniższe wartości poza jakąś wartość, jeśli jestem pewien, że nie ma już znaczenia. W przypadku liczb zmiennoprzecinkowych (lub zbyt wielu liczb całkowitych) prawdopodobnie utworzyłbym segmenty / przedziały, a następnie zastosowałbym to samo podejście, co w przypadku liczb całkowitych. (Przybliżony) tryb i obliczanie mediany niż staje się łatwe, na podstawie tabeli częstotliwości.

Zmienne losowe z rozkładem normalnym

Jeśli ma rozkład normalny , użyłbym średniej próby populacji , wariancji , skośności i kurtozy jako estymatorów maksymalnego prawdopodobieństwa dla małego podzbioru. Algorytmy (on-line) do ich obliczania, już teraz. Np. Czytaj kilkaset tysięcy lub milionów punktów danych, aż błąd oszacowania stanie się wystarczająco mały. Po prostu upewnij się, że wybierasz losowo ze swojego zestawu (np. Nie wprowadzasz odchylenia, wybierając pierwsze 100 000 wartości). To samo podejście można również zastosować do estymacji trybu i mediany dla przypadku normalnego (w obu przypadkach średnia z próby jest estymatorem).

Dalsze komentarze

Wszystkie powyższe algorytmy mogą być uruchamiane równolegle (w tym wiele algorytmów sortowania i selekcji, np. QuickSort i QuickSelect), jeśli to pomaga.

Zawsze zakładałem (z wyjątkiem sekcji dotyczącej rozkładu normalnego), że mówimy o momentach próbkowania, medianie i trybie, a nie o estymatorach momentów teoretycznych przy znanym rozkładzie.

Ogólnie rzecz biorąc, próbkowanie danych (tj. Patrzenie tylko na podzbiór) powinno być całkiem skuteczne, biorąc pod uwagę ilość danych, o ile wszystkie obserwacje są realizacjami tej samej zmiennej losowej (mają te same rozkłady) i momentów, mody i mediana faktycznie istnieje dla tego rozkładu. Ostatnie zastrzeżenie nie jest nieszkodliwe. Na przykład średnia (i wszystkie wyższe momenty) dla rozkładu Cauchy'ego nie istnieją. W takim przypadku średnia z próbki „małego” podzbioru może znacznie odbiegać od średniej z całej próbki.

Używam tych przyrostowych / rekurencyjnych średnich i median estymatorów, które używają stałej pamięci:

mean += eta * (sample - mean)
median += eta * sgn(sample - median)

gdzie eta to mały parametr szybkości uczenia się (np. 0,001), a sgn () to funkcja signum, która zwraca jedną z wartości {-1, 0, 1}. (Użyj stałej eta, jeśli dane są niestacjonarne i chcesz śledzić zmiany w czasie; w przeciwnym razie dla źródeł stacjonarnych możesz użyć czegoś takiego jak eta = 1 / n dla estymatora średniej, gdzie n to liczba próbek widzianych w ten sposób daleko ... niestety wydaje się, że nie działa to w przypadku estymatora mediany.)

Ten typ estymatora średniej przyrostowej wydaje się być używany wszędzie, np. W regułach uczenia się sieci neuronowych bez nadzoru, ale wersja mediany wydaje się znacznie mniej powszechna, pomimo jej zalet (odporność na wartości odstające). Wydaje się, że wersja mediany mogłaby służyć jako zamiennik dla estymatora średniej w wielu zastosowaniach.

Bardzo chciałbym zobaczyć estymator trybu przyrostowego o podobnej formie ...

AKTUALIZACJA

Właśnie zmodyfikowałem estymator przyrostowej mediany, aby oszacować dowolne kwantyle. Ogólnie rzecz biorąc, funkcja kwantyla ( http://en.wikipedia.org/wiki/Quantile_function ) podaje wartość, która dzieli dane na dwa ułamki: p i 1-p. Następujący szacuje tę wartość w sposób przyrostowy:

quantile += eta * (sgn(sample - quantile) + 2.0 * p - 1.0)

Wartość p powinna mieścić się w granicach [0,1]. To zasadniczo przesuwa symetryczne wyjście funkcji sgn () {-1,0,1} tak, aby przechylało się w jedną stronę, dzieląc próbki danych na dwa pojemniki o nierównej wielkości (ułamki p i 1-p danych są mniejsze / większe niż szacunek kwantylowy). Zauważ, że dla p = 0,5, sprowadza się to do estymatora mediany.

Zaimplementowałem algorytm P-Square do dynamicznego obliczania kwantyli i histogramów bez przechowywania obserwacji w zgrabnym module Pythona, który napisałem o nazwie LiveStats . Powinien dość skutecznie rozwiązać twój problem. Biblioteka obsługuje wszystkie wspomniane statystyki z wyjątkiem trybu. Nie znalazłem jeszcze satysfakcjonującego rozwiązania do estymacji trybu.

Ryan, obawiam się, że nie robią średniej i wariancji prawo ... To wymyślił kilka tygodni temu tutaj . Jedną z mocnych stron wersji online (która faktycznie nosi nazwę metody Welforda) jest fakt, że jest ona wyjątkowo dokładna i stabilna, zobacz dyskusję tutaj . Jedną z mocnych stron jest fakt, że nie musisz przechowywać całkowitej sumy ani całkowitej sumy kwadratów ...

Nie przychodzi mi do głowy żadne podejście do trybu i mediany w trybie on-line, które wydaje się wymagać od razu przeanalizowania całej listy. Ale może się zdarzyć, że podobne podejście niż to dla wariancji i średniej będzie działać również w przypadku skośności i kurtozy ...

Cytowany w pytaniu artykuł Wikipedii zawiera formuły obliczania skośności i kurtozy on-line.

W przypadku trybu - jak sądzę - nie ma możliwości zrobienia tego w trybie online. Czemu? Załóżmy, że wszystkie wartości wprowadzonych danych są inne niż ostatnia, która powiela poprzednią. W takim przypadku musisz zapamiętać wszystkie wartości już widziane na wejściu, aby wykryć, że ostatnia wartość powiela wartość widzianą wcześniej i sprawia, że ​​jest ona najczęstsza.

Mediana jest prawie taka sama - aż do ostatniego wejścia nie wiadomo, jaka wartość stanie się medianą, jeśli wszystkie wartości wejściowe będą różne, ponieważ może to być przed lub po bieżącej medianie. Jeśli znasz długość danych wejściowych, możesz znaleźć medianę bez przechowywania wszystkich wartości w pamięci, ale nadal będziesz musiał przechowywać wiele z nich (chyba około połowy), ponieważ zła sekwencja wejściowa może znacznie przesunąć medianę w w drugiej połowie może mieć dowolną wartość z pierwszej połowy mediany.

(Zwróć uwagę, że odnoszę się tylko do dokładnych obliczeń).

Jeśli masz miliardy punktów danych, nie jest prawdopodobne, że potrzebujesz dokładnych odpowiedzi, w przeciwieństwie do zamkniętych odpowiedzi. Ogólnie rzecz biorąc, jeśli masz miliardy punktów danych, podstawowy proces, który je generuje, będzie prawdopodobnie podlegał pewnej statystycznej stacjonarności / ergodyczności / właściwościach mieszania. Może również mieć znaczenie, czy spodziewasz się, że dystrybucje będą w miarę ciągłe, czy nie.

W takich okolicznościach istnieją algorytmy dla on-line, małej pamięci, estymacji kwantyli (mediana jest szczególnym przypadkiem kwantyla 0,5), a także modów, jeśli nie potrzebujesz dokładnych odpowiedzi. To jest aktywne pole statystyk.

przykład oszacowania kwantyli: http://www.computer.org/portal/web/csdl/doi/10.1109/WSC.2006.323014

przykład oszacowania trybu: Bickel DR. Solidne estymatory modu i skośności danych ciągłych. Statystyka obliczeniowa i analiza danych. 2002; 39: 153–163. doi: 10.1016 / S0167-9473 (01) 00057-3.

Są to aktywne pola statystyki obliczeniowej. Wchodzisz w dziedziny, w których nie ma jednego najlepszego dokładnego algorytmu, ale różnorodność z nich (tak naprawdę estymatory statystyczne), które mają różne właściwości, założenia i wydajność. To matematyka eksperymentalna. Na ten temat są prawdopodobnie setki, a nawet tysiące artykułów.

Ostatnie pytanie brzmi, czy naprawdę potrzebujesz samych skośności i kurtozy, czy też bardziej prawdopodobne jest, że niektóre inne parametry mogą być bardziej wiarygodne przy charakteryzowaniu rozkładu prawdopodobieństwa (zakładając, że masz rozkład prawdopodobieństwa!) Czy spodziewasz się Gaussa?

Czy masz sposoby na czyszczenie / wstępne przetwarzanie danych, aby były w większości gaussowskie? (na przykład kwoty transakcji finansowych są często nieco gaussowskie po wzięciu logarytmów). Czy spodziewasz się skończonych odchyleń standardowych? Czy spodziewasz się grubych ogonów? Czy zależy Ci na ilościach w ogonach czy w masie?

Wszyscy powtarzają, że nie da się tego zrobić w trybie online, ale to po prostu nieprawda. Oto artykuł opisujący algorytm do rozwiązania tego właśnie problemu, wynaleziony w 1982 roku przez Michaela E. Fischera i Stevena L. Salzberga z Uniwersytetu Yale. Z artykułu:

Algorytm znajdowania większości wykorzystuje jeden ze swoich rejestrów do tymczasowego przechowywania pojedynczego elementu ze strumienia; ten element jest aktualnym kandydatem na element większościowy. Drugi rejestr to licznik zainicjalizowany na 0. Dla każdego elementu strumienia prosimy algorytm o wykonanie następującej procedury. Jeśli licznik odczytuje 0, zainstaluj bieżący element strumienia jako nowego kandydata większościowego (zastępując każdy inny element, który może już znajdować się w rejestrze). Następnie, jeśli bieżący element pasuje do kandydata większościowego, zwiększ licznik; w przeciwnym razie zmniejsz licznik. W tym momencie cyklu, jeśli widziana do tej pory część strumienia ma element większościowy, ten element znajduje się w rejestrze kandydującym, a licznik ma wartość większą niż 0. A jeśli nie ma elementu większości? Bez drugiego przejścia przez dane - co nie jest możliwe w środowisku strumieniowym - algorytm nie zawsze może udzielić jednoznacznej odpowiedzi w tej sytuacji. Obiecuje jedynie prawidłowe zidentyfikowanie elementu większościowego, jeśli taki istnieje.

Można go również rozszerzyć, aby znaleźć górne N z większą ilością pamięci, ale powinno to rozwiązać problem w trybie.

Ostatecznie, jeśli nie masz parametrycznej wiedzy a priori o rozkładzie, myślę, że musisz przechowywać wszystkie wartości.

To powiedziawszy, jeśli nie masz do czynienia z jakąś patologiczną sytuacją, środek zaradczy (Rousseuw i Bassett 1990) może być wystarczająco dobry do twoich celów.

Po prostu polega na obliczeniu mediany partii median.

mediany i trybu nie można obliczyć online przy użyciu tylko stałej dostępnej przestrzeni. Jednakże, ponieważ mediana i mod są i tak bardziej „opisowe” niż „ilościowe”, można je oszacować, np. Poprzez próbkowanie zbioru danych.

Jeśli w dłuższej perspektywie dane mają rozkład normalny, możesz po prostu użyć średniej do oszacowania mediany.

Możesz również oszacować medianę za pomocą następującej techniki: ustal medianę oszacowania M [i] dla każdego, powiedzmy, 1 000 000 wpisów w strumieniu danych, tak że M [0] jest medianą pierwszego miliona wpisów, M [1] mediana drugiego miliona pozycji itd. Następnie użyj mediany M [0] ... M [k] jako estymatora mediany. To oczywiście oszczędza miejsce i możesz kontrolować, ile chcesz użyć miejsca, „dostrajając” parametr 1 000 000. Można to również uogólnić rekurencyjnie.

OK stary, spróbuj tych:

dla c ++:

double skew(double* v, unsigned long n){
    double sigma = pow(svar(v, n), 0.5);
    double mu = avg(v, n);

    double* t;
    t = new double[n];

    for(unsigned long i = 0; i < n; ++i){
        t[i] = pow((v[i] - mu)/sigma, 3);
    }

    double ret = avg(t, n);

    delete [] t;
    return ret;
}

double kurt(double* v, double n){
    double sigma = pow(svar(v, n), 0.5);
    double mu = avg(v, n);

    double* t;
    t = new double[n];

    for(unsigned long i = 0; i < n; ++i){
        t[i] = pow( ((v[i] - mu[i]) / sigma) , 4) - 3;
    }

    double ret = avg(t, n);

    delete [] t;
    return ret;
}

jeśli mówisz, że możesz już obliczyć wariancję próbki (svar) i średnią (avg), kierujesz je do swoich funkcji, aby to zrobić.

Spójrz też na przybliżenie Pearsona. na tak dużym zbiorze danych byłoby całkiem podobnie. 3 (średnia - mediana) / odchylenie standardowe masz medianę jako max - min / 2

dla trybu zmiennoprzecinkowego nie ma znaczenia. zazwyczaj umieszcza się je w pojemnikach o znacznej wielkości (np. 1/100 * (max - min)).

Ten problem został rozwiązany przez Pebay i wsp .:

https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf

Zwykle używałbym wiader, które mogą być adaptacyjne. Rozmiar łyżki powinien odpowiadać wymaganej dokładności. Następnie, gdy przychodzi każdy punkt danych, dodajesz jeden do liczby odpowiedniego segmentu. Powinny one dać proste przybliżenie mediany i kurtozy, licząc każdy segment jako jego wartość ważoną liczbą.

Jedynym problemem może być utrata rozdzielczości w liczbach zmiennoprzecinkowych po miliardach operacji, tj. Dodanie jednej nie zmienia już wartości! Aby to obejść, jeśli maksymalny rozmiar wiadra przekracza pewien limit, można usunąć dużą liczbę ze wszystkich obliczeń.

for j in range (1,M):
    y=np.zeros(M) # build the vector y
    y[0]=y0

    #generate the white noise
    eps=npr.randn(M-1)*np.sqrt(var)

    #increment the y vector
    for k in range(1,T):
        y[k]=corr*y[k-1]+eps[k-1]

    yy[j]=y

list.append(y)