Algorytm do dynamicznego monitorowania kwantyli

24

Chcę oszacować kwantyl niektórych danych. Dane są tak ogromne, że nie można ich zapisać w pamięci. A dane nie są statyczne, wciąż pojawiają się nowe dane. Czy ktoś zna jakiś algorytm do monitorowania kwantyli danych obserwowanych do tej pory przy bardzo ograniczonej pamięci i obliczeniach? Uważam, że algorytm P2 jest użyteczny, ale nie działa zbyt dobrze w przypadku moich danych, które są bardzo rozproszone.

sinoTrinity
źródło
Niektóre pomysły (w kontekście szacowania median) można znaleźć w wątku na stronie stats.stackexchange.com/q/346/919 .
whuber
3
To pytanie jest crossposted na math.SE.
kardynał

Odpowiedzi:

16

Algorytm P2 to miłe znalezisko. Działa poprzez dokonanie kilku oszacowań kwantyla, okresową aktualizację i zastosowanie interpolacji kwadratowej (nieliniowej, nie sześciennej) do oszacowania kwantyla. Autorzy twierdzą, że interpolacja kwadratowa działa lepiej na ogonach niż interpolacja liniowa, a sześcienny stałby się zbyt wybredny i trudny.

Nie podajesz dokładnie, w jaki sposób to podejście zawodzi w przypadku twoich „grubościennych” danych, ale łatwo zgadnąć: szacunki ekstremalnych kwantyli dla gruboziarnistych dystrybucji będą niestabilne, dopóki nie zgromadzi się dużej ilości danych. Ale będzie to stanowić problem (w mniejszym stopniu), nawet jeśli miałbyś przechowywać wszystkie dane, więc nie oczekuj cudów!

W każdym razie, dlaczego nie ustawić pomocniczych znaczników - nazwijmy je i x 6 - w przypadku których masz pewność, że kwantyl będzie leżał, i zapisz wszystkie dane, które leżą między x 0 a x 6 ? Gdy bufor się zapełni, będziesz musiał zaktualizować te znaczniki, zawsze zachowując x 0x 6 . Prosty algorytm do tego celu można opracować na podstawie kombinacji (a) aktualnego oszacowania P2 kwantyla i (b) przechowywanych zliczeń liczby danych mniejszej niż x 0 i liczby danych większej niż x 6x0x6x0x6x0x6x0x6. W ten sposób możesz z dużą pewnością oszacować kwantyl tak samo dobrze, jak gdyby cały zestaw danych był zawsze dostępny, ale potrzebujesz tylko stosunkowo małego bufora.

W szczególności proponuję strukturę danych celu utrzymania częściowej informacji o sekwencji n wartości danych x 1 , x 2 , , x n . Tutaj y jest połączoną listą(k,y,n)nx1,x2,,xny

y=(x[k+1](n)x[k+2](n)x[k+m](n)).

W tej notacji oznacza i- najmniejszą z dotychczas odczytanych wartości n x . m jest stałą wielkością bufora y .x[i](n)ithn xmy

Algorytm zaczyna się od wypełnienia napotkanymi pierwszymi m wartościami danych i umieszczenia ich w posortowanej kolejności, od najmniejszej do największej. Niech q będzie kwantylem do oszacowania; np. q = 0,99. Po odczytaniu x n + 1 możliwe są trzy działania:ymqqxn+1

  • Jeżeli , przyrost k .xn+1<x[k+1](n)k

  • Jeśli , nic nie rób.xn+1>x[k+m](n)

  • W przeciwnym razie wstaw do y .xn+1y

W każdym razie przyrost .n

Procedura wstawiania umieszcza w y w posortowanej kolejności, a następnie eliminuje jedną z ekstremalnych wartości w y :xn+1yy

  • Jeśli , a następnie usunąć x ( n ) [ k + 1 ] z y a przyrost K ;k+m/2<nqx[k+1](n)yk

  • W przeciwnym razie usuń z y .x[k+m](n)y

Pod warunkiem, że jest wystarczająco duży, procedura ta będzie zawierała prawdziwe kwantyle rozkładu z dużym prawdopodobieństwem. Na dowolnym etapie n można to oszacować w zwykły sposób w kategoriach x ( n ) [ q n] i x ( n ) [ q n] , które prawdopodobnie będą znajdować się w y . (Uważam, że m musi być skalowane tylko jako pierwiastek kwadratowy maksymalnej ilości danych ( Nmnx[qn](n)x[qn](n)ymN), ale nie przeprowadziłem rygorystycznej analizy, aby to udowodnić.) W każdym razie algorytm wykryje, czy się udało (porównując i ( k + m ) / n do q ).k/n(k+m)/nq

Testowanie do 100 000 wartości przy użyciu iq=.5(najtrudniejszy przypadek) wskazuje, że ten algorytm ma 99,5% skuteczności w uzyskaniu prawidłowej wartości x ( n ) [ q n] . Dla strumieniaN=10 12 wartości wymagałoby to bufora tylko dwóch milionów (ale trzy lub cztery miliony byłoby lepszym wyborem). Użycie posortowanej podwójnie połączonej listy dla bufora wymagaO(log(m=2Nq=.5x[qn](n)N=1012=wysiłekO(log(N))podczas identyfikowania i usuwania maks. Lub min tooperacjeO(1). Stosunkowo drogie wstawienie zwykle wymaga wykonania tylkoO(O(log(N))O(log(N))O(1)razy. Zatem koszty obliczeniowe tego algorytmu wynosząO(N+O(N)w czasie iO(O(N+Nlog(N))=O(N)w magazynie.O(N)

Whuber
źródło
Jest to rozszerzona praca algorytmu P2. [link] sim.sagepub.com/content/49/4/159.abstract . Pamięć jest wciąż za duża dla mojej aplikacji, która działa na małych czujnikach o łącznej pojemności 10 KB. Mogę zużyć najwyżej kilkaset bajtów tylko do oszacowania kwantylowego.
sinoTrinity
@ whuber Właściwie implementuję rozszerzone P2 i testuję je z wygenerowanymi próbkami z różnych dystrybucji, takich jak jednolite i wykładnicze, gdzie działa świetnie. Ale kiedy stosuję je w stosunku do danych z mojej aplikacji, których rozkład jest nieznany, czasami nie zbiega się i daje błąd względny (abs (oszacowanie - rzeczywisty) / rzeczywisty) do 300%.
sinoTrinity,
2
@sino Jakość algorytmu w porównaniu do wykorzystania wszystkich danych nie powinna zależeć od ciężkości ogonów. Bardziej sprawiedliwym sposobem pomiaru błędu jest: niech będzie empirycznym cdf. Dla oszacowania q na q percentyla, jaka jest różnica między F ( q ) i F ( q ) ? Jeśli jest rzędu 1 / n , radzisz sobie bardzo dobrze. Innymi słowy, jaki percentyl zwraca algorytm P2 dla danych? Fq^qF(q^)F(q)1/n
whuber
Masz rację. Właśnie zmierzyłem F (qˆ) i F (q) w przypadku, o którym wspomniałem, z błędem względnym do 300%. Dla q 0,7 q7 wynosi prawie 0,7, co powoduje pomijalny błąd. Jednak dla q wynoszącego 0,9 qˆ wydaje się wynosić około 0,95. Myślę, że dlatego mam ogromny błąd do 300%. Wiesz, dlaczego to 0,95, a nie 0,9? BTW, czy mogę zamieścić tutaj rysunek i jak mogę opublikować wzór matematyczny tak jak Ty?
sinoTrinity
2
@ whuber Jestem całkiem pewien, że moja implementacja jest zgodna z rozszerzonym P2. 0,9 nadal idzie do 0,95 lub nawet więcej, gdy jednocześnie szacuję 0,8, 0,85, 0,9, 0,95 kwantyle. Jednak 0,9 zbliża się bardzo do 0,9, jeśli jednocześnie śledzone są kwantyle 0,8, 0,85, 0,9, 0,95 i 1,0 .
sinoTrinity
5

O(N)

Instead of tracking the quantiles at 0, p/2, p, (1+p)/2, and 1, as the original P2 algorithm suggests, you could simply keep track of more quantiles (but still a constant number). It looks like the algorithm allows for that in a very straightforward manner; all you need to do is compute the correct "bucket" for incoming points, and the right way to update the quantiles (quadratically using adjacent numbers).

Say you keep track of 25 points. You could try tracking the quantile at 0, p/12, , p11/12, p, p+(1p)/12, , p+11(1p)/12, 1 (picking the points equidistantly in between 0 and p, and between p and 1), or even using 22 Chebyshev nodes of the form p/2(1+cos(2i1)π22) and p+(1p)/2(1+cos(2i1)π22). If p is close to 0 or 1, you could try putting fewer points on the side where there is less probability mass and more on the other side.

If you decide to pursue this, I (and possibly others on this site) would be interested in knowing if it works...

Erik P.
źródło
+1 I think this is a great idea given the OP's constraints. All one can hope for is an approximation, so the trick is to pick bins that have a high likelihood of being narrow and containing the desired quantile.
whuber
3

Press et al., Numerical Recipes 8.5.2 "Single-pass estimation of arbitrary quantiles" p. 435, give a c++ class IQAgent which updates a piecewise-linear approximate cdf.

denis
źródło
books.google.com/… for a version that doesn't require Flash.
ZachB
2

This can be adapted from algorithms that determine the median of a dataset online. For more information, see this stackoverflow post - /programming/1387497/find-median-value-from-a-growing-set

benhamner
źródło
The computational resources required of the algorithm you link to are unnecessarily large and do not meet the requirements of this question.
whuber
2

I'd look at quantile regression. You can use it to determine a parametric estimate of whichever quantiles you want to look at. It make no assumption regarding normality, so it handles heteroskedasticity pretty well and can be used one a rolling window basis. It's basically an L1-Norm penalized regression, so it's not too numerically intensive and there's a pretty full featured R, SAS, and SPSS packages plus a few matlab implementations out there. Here's the main and the R package wikis for more info.

Edited:

Check out the math stack exchange crosslink: Someone sited a couple of papers that essentially lay out the very simple idea of just using a rolling window of order statistics to estimate quantiles. Literally all you have to do is sort the values from smallest to largest, select which quantile you want, and select the highest value within that quantile. You can obviously give more weight to the most recent observations if you believe they are more representative of actual current conditions. This will probably give rough estimates, but it's fairly simple to do and you don't have to go through the motions of quantitative heavy lifting. Just a thought.

Marc
źródło
1

It is possible to estimate (and track) quantiles on an on-line basis (the same applies to the parameters of a quantile regression). In essence, this boils down to stochastic gradient descent on the check-loss function which defines quantile-regression (quantiles being represented by a model containing only an intercept), e.g. updating the unknown parameters as and when observations arrive.

See the Bell Labs paper "Incremental Quantile Estimation for Massive Tracking" ( ftp://ftp.cse.buffalo.edu/users/azhang/disc/disc01/cd1/out/papers/kdd/p516-chen.pdf)

Ludo
źródło