Algorytmy" On-line " (iterator) do szacowania mediany statystycznej, trybu, skośności, kurtozy?

Używam tych przyrostowych / rekurencyjnych estymatorów średnich i median, które używają stałego przechowywania:

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

Gdzie eta jest małym parametrem szybkości uczenia się (np. 0,001), a sgn () jest funkcją signum, która zwraca jeden z {-1, 0, 1}. (Użyj stałej eta Jeśli dane są niestacjonarne I chcesz śledzić zmiany w czasie; w przeciwnym razie dla stacjonarnych źródeł możesz użyć czegoś w rodzaju eta=1/N dla estymatora średniego, gdzie n jest liczbą próbki widziane do tej pory... niestety, wydaje się, że nie działa to dla estymatora mediany.)

Ten typ przyrostowego estymatora średniej wydaje się być stosowany wszędzie, np. w nienadzorowanych regułach uczenia się sieci neuronowych, ale wersja mediana wydaje się znacznie mniej powszechna, pomimo swoich korzyści (odporność na wartości odstające). Wydaje się, że wersja mediana może być stosowany jako zamiennik dla estymatora średniej w wielu zastosowaniach.

Chciałbym zobaczyć Estymator trybu przyrostowego o podobnej formie...

UPDATE

Właśnie zmodyfikowałem Przyrostowy Estymator mediany, aby oszacować dowolne kwantyle. Ogólnie rzecz biorąc, funkcja kwantyl ( http://en.wikipedia.org/wiki/Quantile_function ) podaje wartość, która dzieli dane na dwa ułamki: p I 1-p. poniższy oszacowuje tę wartość stopniowo:

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} aby pochylić się w kierunku jednej strony, dzielenie próbek danych na dwa pojemniki o nierównej wielkości (ułamki p i 1 - P danych są mniejsze niż / większe niż oszacowanie kwantylu, odpowiednio). Należy zauważyć, że dla p=0,5 zmniejsza się to do estymatora mediany.

Zaimplementowałem P-Kwadratowy algorytm do dynamicznego obliczania Kwantyli i histogramów bez przechowywania obserwacji w zgrabnym module Pythona, który napisałem o nazwie LiveStats. Powinno to dość skutecznie rozwiązać twój problem. Biblioteka obsługuje wszystkie wymienione statystyki z wyjątkiem trybu. Nie znalazłem jeszcze zadowalającego rozwiązania dla oszacowania trybu.

Ryan, obawiam się, że nie robisz tego źle i źle... To pojawiło się kilka tygodni temu tutaj . A jednym z mocnych punktów wersji online (która w rzeczywistości nosi nazwę metody Welforda) jest fakt, że jest ona szczególnie dokładna i stabilna, patrz dyskusja Tutaj . Jednym z mocnych punktów jest fakt, że nie trzeba przechowywać całkowitą sumę lub całkowitą sumę kwadratów...

Nie przychodzi mi do głowy żadne podejście on-line do trybu i mediany, które wydaje się, że wymaga rozważenia całej listy na raz. Ale może być bardzo dobrze, że podobne podejście niż to dla wariancji i średniej będzie działać również dla skewness i kurtosis...

Cytowany w pytaniu artykuł Wikipedii zawiera wzory kalkuracji skośności i kurtozy on-line.

Dla mode-chyba-nie da się tego zrobić on-line. Dlaczego? Załóżmy, że wszystkie wartości Twojego wejścia są inne niż ostatnia, która duplikuje poprzednią. W tym przypadku musisz zapamiętać wszystkie wartości już widoczne na wejściu, aby wykryć, że ostatnia wartość powiela wartość widzianą i czyni ją najczęstszą.

Dla mediany jest prawie tak samo - do ostatniego wejścia nie wiesz, jaka wartość stanie się medianą, jeśli wszystkie wartości wejściowe są różne, ponieważ może być przed lub po bieżącej medianie. Jeśli znasz długość wejścia, możesz znaleźć medianę bez przechowywania wszystkich wartości w pamięci, ale nadal będziesz musiał przechowywać wiele z nich (myślę, że około połowy), ponieważ zła Sekwencja wejściowa może mocno przesunąć medianę w drugiej połowie, prawdopodobnie dokonując dowolnej wartości z pierwszej połowy. mediana.

(zauważ, że odwołuję się tylko do dokładnych obliczeń.)

Jeśli masz miliardy punktów danych, to nie jest prawdopodobne, że potrzebujesz dokładnych odpowiedzi, w przeciwieństwie do bliskich odpowiedzi. Ogólnie rzecz biorąc, jeśli masz miliardy punktów danych, podstawowy proces, który je generuje, prawdopodobnie będzie posłuszny jakiejś stateczności statystycznej / ergodyczności / właściwości mieszania. Również może mieć znaczenie, czy spodziewasz się, że dystrybucje będą w miarę ciągłe, czy nie.

W tych okolicznościach istnieją algorytmy do estymacji on-line, niskiej pamięci, kwantyli (mediana jest szczególnym przypadkiem 0,5 kwantyla), a także trybów, jeśli nie potrzebujesz dokładnych odpowiedzi. Jest to aktywna dziedzina statystyki.

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

Przykład estymacji trybu: Bickel DR. solidne estymatory trybu i skośności danych ciągłych. Statystyka obliczeniowa i analiza danych. 2002;39:153–163. podoba mi się! do obserwowanych nr: 64737

Te są aktywnymi dziedzinami statystyki obliczeniowej. Dostajesz się do dziedzin, w których nie ma jednego najlepszego dokładnego algorytmu, ale ich różnorodność (estymatory statystyczne, w rzeczywistości), które mają różne właściwości, założenia i wydajność. To eksperymentalna matematyka. Są prawdopodobnie setki do tysięcy artykułów na ten temat.

Ostatnie pytanie brzmi, czy naprawdę potrzebujesz skewness i kurtosis samodzielnie, czy raczej jakieś inne parametry, które mogą być bardziej wiarygodny w charakteryzowaniu rozkładu prawdopodobieństwa (zakładając, że masz rozkład prawdopodobieństwa!). Spodziewasz się Gaussa?

Czy masz sposoby czyszczenia / wstępnego przetwarzania danych, aby były głównie Gaussiańskie? (na przykład kwoty transakcji finansowych są często nieco gaussowskie po wzięciu logarytmów). Czy spodziewasz się skończonych odchyleń standardowych? Spodziewasz się grubych ogonów? Czy ilości, na których Ci zależy, są w ogonach czy luzem?

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

Algorytm wyszukiwania większości używa jednego ze swoich rejestrów do tymczasowego przechowywanie pojedynczej pozycji ze strumienia; pozycja ta jest bieżącą kandydat na element większościowy. Drugi rejestr jest licznikiem zainicjalizowana na 0. Dla każdego elementu strumienia pytamy algorytm aby wykonać następującą rutynę. Jeśli licznik odczyta 0, zainstaluj aktualnego elementu strumienia jako nowego kandydata większościowego (wypierającego dowolne inny element, który może już być w rejestrze). Następnie, jeśli bieżący element pasuje do większości kandydatów, zwiększ licznik; w przeciwnym razie zmniejsz licznik. W tym momencie cyklu, jeśli część strumienia widzianego do tej pory ma większość element, ten element jest w rejestrze kandydatów, a licznik posiada wartość większą niż 0. Co jeśli nie ma elementu większościowego? Bez drugiego przepuszczania danych-co nie jest możliwe w środowisku strumieniowym - algorytm nie zawsze może dać jednoznaczną odpowiedź w tym okoliczności. Jedynie obiecuje prawidłowo zidentyfikować większość element, jeśli istnieje.

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

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

To powiedziane, jeśli nie masz do czynienia z jakąś patologiczną sytuacją, remedian (Rousseuw and Bassett 1990) może być wystarczająco dobry dla Twoich celów.

Bardzo prosto polega na obliczeniu mediany partii medianów.

Mediana i tryb nie mogą być obliczane online przy użyciu tylko stałej dostępnej przestrzeni. Ponieważ jednak mediana i tryb są i tak bardziej "opisowe" niż "ilościowe", można je oszacować np. poprzez pobranie zestawu danych.

Jeśli dane są normalne dystrybuowane w dłuższej perspektywie, możesz po prostu użyć swojej średniej do oszacowania mediany.

Można również oszacować medianę za pomocą następującej techniki: ustalić medianę estymacji M [i] dla każdego, powiedzmy, 1,000,000 wpisów w danych strumień tak, że M[0] jest mediana pierwszego miliona wpisów, M [1] mediana drugiego miliona wpisów itp. Następnie należy użyć mediany M[0]...M [k] jako Estymator 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 koleś 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;
}

Gdzie można już obliczyć wariancję próbki (svar) i średnią (avg) kieruj je do swoich funkcji, aby to zrobić.

Dla trybu floats nie ma znaczenia. zazwyczaj wsadza się je do pojemników o wielkości sginificant (jak 1/100 * (max - min)).

Ten problem został rozwiązany przez Pebay et al:

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

Zwykle używam wiadra, które mogą być adaptacyjne. Rozmiar łyżki powinien być dokładnością, której potrzebujesz. Następnie, gdy pojawi się każdy punkt danych, dodajesz go do licznika odpowiedniego zasobnika. Powinny one dać proste przybliżenia do mediany i kurtozy, licząc każdy kubeł jako jego wartość ważoną przez jego liczbę.

Jedynym problemem może być utrata rozdzielczości w zmiennoprzecinkowej po miliardach operacji, tzn. dodanie jednego nie zmienia już wartości! Aby obejść to, jeśli maksymalny rozmiar łyżki przekracza pewien limit, który można zdjąć dużą liczbę ze wszystkich liczeń.

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)