Statystyka odpornościowa

Z Wikipedii, wolnej encyklopedii
Skocz do: nawigacji, wyszukiwania

Statystyka odpornościowa[1] lub odporne metody statystyczne[2] (ang. robust statistics) – gałąź statystyki, obejmująca metody projektowane pod kątem odporności na niewielkie odejście od założeń modelu (szczególnie występowanie obserwacji odstających) lub rezygnacji z niektórych założeń.

Wprowadzenie[edytuj | edytuj kod]

Klasyczne metody statystyczne w dużym stopniu polegają na założeniach, które nie zawsze są spełnione w praktyce. W szczególności często zakłada się, że dane mają rozkład normalny (choćby w przybliżeniu), lub próbka jest dostatecznie duża, aby dało się zastosować centralne twierdzenie graniczne do otrzymania rozkładu normalnego błędu estymatora.

Estymatory bazujące na rozkładzie normalnym lub wielowymiarowym rozkładzie normalnym (takie jak np. współczynniki regresji liniowej, współczynnik korelacji Pearsona, czy odchylenie standardowe) mają zwykle tę cechę, iż obserwacje bardziej oddalone od średniej są bardziej wpływowe, czyli wchodzą z większą wagą do wyniku. Ma to groźną konsekwencję: w przypadku obecności w danych obserwacji odstających, czyli np. pomyłek typu przesunięcie przecinka, nawet jedna taka obserwacja może zakłócić wynik analizy.

Ponadto normalność rozkładu jest zawsze tylko przybliżeniem – realne dane są zawsze skwantowane (dyskretne) i ograniczone, a tymczasem zgodnie z rozkładem normalnym każda liczba rzeczywista ma dodatnie prawdopodobieństwo wystąpienia, czyli rozkład ten przewiduje np. istnienie ludzi o wzroście ujemnym lub wielkości galaktyki. Naturalnie przybliżenie to jest często tak dobre, że mimo to opłaca się go wykorzystywać, czasem jednak dane wyraźnie różnią się od rozkładu normalnego i standardowe metody statystyczne mogą dawać wątpliwe wyniki.

Istnieją miary "odporności" (ang. robustness) danej metody na naruszenie założeń modelu. Najczęściej używane są punkt załamania (ang. breakdown point) i funkcja wpływu (ang. influence function), opisane dalej.

Spora część odpornych metod statystycznych zastępuje rozkład normalny z metod klasycznych, rozkładem Studenta z niewielką liczbą stopni swobody, czyli dużą skośnością (w praktyce wykorzystuje się rozkłady z 4-6 stopniami swobody), lub stosuje kombinację dwóch lub większej liczby rozkładów.

Przykład: prędkość światła[edytuj | edytuj kod]

Gelman z zespołem w Bayesian Data Analysis (2004) rozważali zbiór danych związany z pomiarami prędkości światła wykonanymi przez Simona Newcomba[3].

Choć rozkład większości obserwacji wygląda na mniej więcej normalny, są też dwie oczywiste obserwacje odstające. Te obserwacje mają duży wpływ na średnią przyciągając ją do siebie i odciągając od średniej pozostałych. W ten sposób jeśli użyjemy średniej jako miary położenia rozkładu, wynik będzie obciążony.

Rozkład średniej zgodnie z centralnym twierdzeniem granicznym zbiega do rozkładu normalnego wraz ze wzrostem liczebności próbki do nieskończoności. Obserwacje odstające mogą jednak sprawić, że rozkład średniej nie będzie normalny nawet dla bardzo dużych prób. Ponadto obserwacje odstające, jeśli stanowią pewien stały procent próby, mogą zmienić parametry granicznego rozkładu normalnego.

Estymacja położenia[edytuj | edytuj kod]

Wykresy poniżej przedstawiają gęstość prawdopodobieństwa dla rozkładu prędkości światła w danych (wykres (a)). Pokazany jest również wykres kwantylowo-kwantylowy. Obserwacje odstające są łatwe do zauważenia.

Wykresy (c) i (d) pokazują rozkład bootstrap średniej arytmetycznej (c) i średniej ucinanej w 10%. Średnia ucinana jest prostym odpornym estymatorem położenia, który bazuje na danych z odciętą częścią (tu 10%) skrajnych obserwacji, obliczając dla pozostałych zwykłą średnią arytmetyczną. Analiza była wykonana dla 10000-elementowej próby bootstrap, zarówno w przypadku średniej arytmetycznej, jak i uciętej.

Jak w oczywisty sposób wynika z wykresu, rozkład średniej jest szerzej rozrzucony (czyli ma większy błąd normalny) w przypadku średniej arytmetycznej. Co więcej rozkład średniej uciętej jest w przybliżeniu normalny, podczas gdy rozkład średniej arytmetycznej jest mocno skośny. W tej próbie 66 obserwacji tylko dwie obserwacje odstające sprawiły, że centralne twierdzenie graniczne nie daje się zastosować (wymagałoby o wiele większej liczby obserwacji).

SpeedOfLight.png

Odporne metody statystyczne, których przykładem jest średnia ucinana, mają dawać lepsze od metod klasycznej statystyki rezultaty w sytuacjach występowania obserwacji odstających, lub ogólniej, gdy założenia modeli parametrycznych nie do końca są spełnione.

Chociaż średnia ucinana lepiej spisuje się w tym przykładzie od średniej arytmetycznej, istnieją jeszcze lepsze odporne estymatory. Średnia, mediana i średnia ucinana są szczególnymi przypadkami M-estymatorów, omówionych dalej w artykule.

Estymacja skali[edytuj | edytuj kod]

Obserwacje odstające w danych na temat prędkości światła wpłynęły nie tylko na średnią. Na ogół w celu estymacji skali (rozrzutu danego rozkładu wokół średniej) używane jest odchylenie standardowe, które jest w znacznie większym stopniu podatne na obserwacje odstające, ze względu na drugie potęgi we wzorze, które sprawiają, że obserwacje odstające są bardziej wpływowe.

Wykresy poniżej pokazują rozkład bootstrap odchylenia standardowego, odchylenia medianowego (MAD) i kwantylowego estymatora skali Qn (Rousseeuw i Croux, 1993). Wykresy opierają się na próbie bootstrap liczącej 10000 obserwacji, do której dodano pewien szum o rozkładzie normalnym (tzw. wygładzany bootstrap, ang. smoothed bootstrap). Wykres (a) pokazuje rozkład odchylenia standardowego, wykres (b) rozkład MAD, a wykres (c) rozkład Qn.

SpeedOfLightScale.png

Rozkład odchylenia standardowego jest obciążony i rozproszony w wyniku obecności obserwacji odstających. MAD zachowuje się lepiej, a Qn jest nawet bardziej efektywne od MAD. Ten prosty przykład pokazuje, że w obecności obserwacji odstających odchylenie standardowe nie powinno być stosowane jako estymator skali.

Ręczne wyszukiwanie obserwacji odstających[edytuj | edytuj kod]

Tradycyjnie statystycy ręcznie wyszukiwali obserwacje odstające w danych i usuwali je, potwierdzając ewentualnie w źródle danych, czy faktycznie są to pomyłki. Faktycznie, w naszych danych o prędkości światła łatwo znaleźć dwie obserwacje odstające i usunąć je przed wykonaniem jakichkolwiek dalszych analiz. Jednakże w dzisiejszych czasach zbiory danych często składają się ze zbyt wielu zmiennych i obserwacji aby coś takiego było możliwe. Ponadto może pojawić się zarzut braku obiektywnych kryteriów takiego usuwania.

Jedne obserwacje odstające czasami maskują inne. Niech w zbiorze danych będzie jedna obserwacja odstająca, bardzo daleka od średniej i druga znacznie bliższa. Ponieważ obserwacja daleka zwiększy wartość odchylenia standardowego, więc bazując na odchyleniu nie zauważymy że obserwacja bliższa także jest odstająca. Jeśli jednak usuniemy dalszą obserwację odstającą, wówczas odchylenie zmniejszy się i bliższa obserwacja odstająca stanie się łatwa do odróżnienia. Problem maskowania jest coraz bardziej dotkliwy, gdy zwiększa się złożoność danych.

Odporne metody statystyczne pozwalają na automatyczną detekcję obserwacji odstających, a następnie zmniejszenie ich wagi i oflagowanie lub usunięcie, w znaczący sposób upraszczając proces analizy.

Miary odporności[edytuj | edytuj kod]

Podstawowymi narzędziami używanymi do opisu i pomiaru odporności na obserwacje odstające, są punkt załamania (ang. breakdown point), funkcja wpływu (ang. influence function) i krzywa wrażliwości (ang. sensitivity curve).

Punkt załamania[edytuj | edytuj kod]

Intuicyjnie, punkt załamania estymatora jest proporcją niepoprawnych obserwacji (np. dowolnie dużych obserwacji odstających), które estymator jest w stanie wytrzymać, zanim zacznie dawać dowolnie duże wyniki. Na przykład mając próbę losową, czyli n niezależnych zmiennych losowych (X_1,\dots,X_n)\sim\mathcal{N}(0,1) i odpowiednie ich realizacje x_1,\dots,x_n, możemy użyć wzoru \overline{X_n}:=\tfrac{X_1+\cdots+X_n}{n} do estymacji średniej w populacji. Taki estymator ma punkt załamania zero, gdyż można uczynić średnią \overline{x} dowolnie dużą za pomocą zmiany jednej wartości x_1,\dots,x_n.

Im wyższy jest punkt załamania estymatora, tym bardziej estymator jest odporny. Punkt załamania nie może przekroczyć 50%, gdyż jeśli ponad połowa obserwacji jest wybranych z innej populacji, nie ma powodu aby traktować tę "inną" populację jako odstającą (obserwacje odstające z definicji stanowią mniejszość, w przeciwnym wypadku to reszta odstaje od nich). Istnieją estymatory, które osiągają ten maksymalny pułap (np. mediana ma punkt załamania równy 0.5). Średnia ucinana na poziomie X% (odcięte X% obserwacji z każdej strony rozkładu) ma punkt załamania X%[4].

Przykład: dane o prędkości światła[edytuj | edytuj kod]

W przykładzie prędkości światła, usunięcie dwóch najniższych wyników powoduje, że średnia wzrasta z 26.2 do 27.75 (zmiana o 1.55). Estymacja skali stworzona metodą Qn wynosi 6.3. Dzieląc to przez pierwiastek z wielkości próby możemy uzyskać odporny estymator błędu standardowego, wynoszący tutaj 0.78. Tak więc zmiana średniej związana z dwiema obserwacjami odstającymi jest tu prawie dwa razy większa od błędu standardowego.

Średnia ucinana na poziomie 10% wynosi tu 27.43. Usuwając dwie najmniejsze obserwacje uzyskujemy za jej pomocą 27.67 (zmiana o 0.24). Ewidentnie średnia ucinana jest mniej wrażliwa na obserwacje odstające, co wiąże się z jej wyższym punktem załamania.

Empiryczna funkcja wpływu[edytuj | edytuj kod]

Funkcja Tukeya (ang. Tukey's biweight function)

Empiryczna funkcja wpływu pokazuje, w jaki sposób estymator zachowuje się, gdy zmieniamy jedną obserwację w próbie i nie bierzemy pod uwagę ewentualnego pogwałcenia założeń modelu. Ilustracja przedstawia dwuwagową funkcję Tukeya pewnego dobrze zachowującego się estymatora.

Oznaczenia:

  • (\Omega,\mathcal{A},P) to przestrzeń probabilistyczna,
  • (\mathcal{X},\Sigma) to przestrzeń mierzonych stanów,
  • \Theta to przestrzeń parametrów wymiaru p\in\mathbb{N}^*,
  • (\Gamma,S) to przestrzeń pomiarów,
  • \gamma:\Theta\rightarrow\Gamma to funkcja symbolizująca pomiar,
  • \mathcal{F}(\Sigma) to zbiór wszystkich możliwych rozkładów na \Sigma

Na przykład:

  1. (\Omega,\mathcal{A},P) to przestrzeń probabilistyczna,
  2. (\mathcal{X},\Sigma)=(\mathbb{R},\mathcal{B}),
  3. \Theta=\mathbb{R}\times\mathbb{R}^+
  4. (\Gamma,S)=(\mathbb{R},\mathcal{B}),
  5. \gamma:\mathbb{R}\times\mathbb{R}^+\rightarrow\mathbb{R} jest zdefiniowana jako \gamma(x,y)=x.

Definicja empirycznej funkcji wpływu:
Niech n\in\mathbb{N}^* i X_1,\cdots,X_n:(\Omega,  \mathcal{A})\rightarrow(\mathcal{X},\Sigma) są niezależnymi zmiennymi losowymi o identycznym rozkładzie, a (x_1,\cdots,x_n) jest próbą statystyczną z tych zmiennych. T_n:(\mathcal{X}^n,\Sigma^n)\rightarrow(\Gamma,S) jest estymatorem. Niech i\in\{1,\cdots,n\}. Empiryczna funkcja wpływu EIF_i obserwacji i jest zdefiniowana jako:

EIF_i:x\in\mathcal{X}\mapsto T_n(x_1,\cdots,x_{i-1},x,x_{i+1},\cdots,x_n)\in\Gamma

Funkcja ta pokazuje zatem zależność od x\; wyników estymatora po zastąpieniu i-tej obserwacji przez x\;.

Funkcja wpływu i krzywa wrażliwości[edytuj | edytuj kod]

W drugim podejściu zamiast sprawdzać wpływ zmiany pojedynczej obserwacji na wynik estymacji, badany będzie wpływ punktowego zaburzenia na rozkład.

Niech T będzie asymptotycznie nieobciążonym estymatorem pewnego parametru rozkładu F spełniającego założenia modelu.

Niech G będzie pewnym rozkładem który nie spełnia założeń modelu statystycznego, różni się jednak nieznacznie od F

Wpływ na estymator T przejścia z F na G jest wyrażony wzorem:

dF_{G-F}(F) = \lim_{t\rightarrow 0^+}\frac{T(tG+(1-t)F) - T(F)}{t},

który jest pochodną kierunkową T w punkcie F w kierunku G.

Niech x\in\mathcal{X}. \Delta_x będzie miarą prawdopodobieństwa o masie 1 skupionej w punkcie x. Wybieramy G=\Delta_x. Funkcja wpływu jest zdefiniowana przez:

IF(x; T; F):=\lim_{t\rightarrow 0^+}\frac{T(t\Delta_x+(1-t)F) - T(F)}{t}.

Wzór ten opisuje wpływ infinitezymalnej zmiany w punkcie x na poszukiwaną estymatę, standaryzowany masą t zaburzenia (asymptotyczne obciążenie estymatora spowodowane przez zanieczyszczenie obserwacji).

Pożądane właściwości[edytuj | edytuj kod]

Właściwości funkcji wpływu, charakteryzujące najbardziej odporne estymatory:

  1. Skończony punkt odrzucenia (ang. rejection point) \rho^*,
  2. Mała maksymalna wrażliwość na błędy (ang. gross-error sensitivity) \gamma^*,
  3. Mała wrażliwość na lokalne przesunięcie (ang. local-shift sensitivity) \lambda^*.

Punkt odrzucenia[edytuj | edytuj kod]

\rho^*:=\inf_{r>0}\{r:IF(x;T;F)=0, |x|>r\}

Maksymalna wrażliwość[edytuj | edytuj kod]

\gamma^*(T;F) := \sup_{x\in\mathcal{X}}|IF(x; T ; F)|

Wrażliwość na lokalne przesunięcia[edytuj | edytuj kod]

\lambda^*(T;F) := \sup_{(x,y)\in\mathcal{X}^2, x\neq y}\left\|\frac{IF(y ; T; F) - IF(x; T ; F)}{y-x}\right\|

Ten wzór, podobny do definicji stałej Lipschitza, opisuje efekt przesunięcia obserwacji z punktu x do pobliskiego punktu y, czyli innymi słowy dodania obserwacji w punkcie y przy jednoczesnym usunięciu z x.

M-estymatory[edytuj | edytuj kod]

Powstało wiele rozmaitych podejść do konstrukcji odpornych estymatorów, w szczególności R-estymatory i L-estymatory. Obecnie tę dziedzinę zdominowały M-estymatory jako najbardziej ogólne, posiadające wysoki punkt załamania i wysoką efektywność (zobacz Huber (1981)).

M-estymatory są uogólnieniem estymatorów największej wiarygodności (MLE). Od MLE oczekujemy maksymalizacji prawdopodobieństwa \prod_{i=1}^n f(x_i) czyli (równoważnie) minimalizacji \sum_{i=1}^n-\log f(x_i). W 1964 Huber zaproponował uogólnienie tego podejścia poprzez minimalizację wyrażenia \sum_{i=1}^n \rho(x_i), gdzie \rho jest pewną funkcją. MLE są tym samym szczególnym przypadkiem M-estymatorów. Nazwa ta pochodzi od ang. Maximum likelihood type estimators, czyli estymatory typu największej wiarygodności.

Minimalizacja \sum_{i=1}^n \rho(x_i) często może być wykonana przez różniczkowanie \rho i rozwiązanie równania \sum_{i=1}^n \psi(x_i) = 0, gdzie \psi(x) = \frac{d\rho(x)}{dx}, oczywiście pod warunkiem, że \rho jest różniczkowalne.

Zostało zaproponowanych wiele funkcji \rho. Wykresy poniżej pokazują cztery funkcje \rho i odpowiadające im \psi.

RhoFunctions.png

Funkcja \rho(x) opisująca błąd kwadratowy rośnie coraz szybciej, podczas gdy dla błędów bezwzględnych rośnie liniowo. Stosując winsoryzację obserwujemy mieszaninę tych dwóch efektów: dla małych wartości x \rho rośnie z kwadratem x, a kiedy zostanie osiągnięty określony próg (w przykładzie 1,5), wzrost staje się liniowy.

Funkcja Tukey'a (ang. Tukey's biweight) zachowuje się początkowo w podobny sposób, jak funkcja błędu kwadratowego, ale dla większych błędów staje się płaska.

PsiFunctions.png

Właściwości M-estymatorów[edytuj | edytuj kod]

Zostało pokazane, że M-estymatory mają asymptotycznie rozkład normalny, więc o ile można obliczyć ich błąd standardowy, dla dużych prób da się stosować metody wnioskowania statystycznego oparte na rozkładzie normalnym.

Dla małych prób ich rozkład może znacząco różnić się od normalnego, lepiej więc stosować inne metody wnioskowania statystycznego, takie jak bootstrap. Należy jednak ostrożnie projektować schematy badań bootstrap, aby nie przekroczyć (jakkolwiek wysokiego) punktu załamania estymatora.

Funkcja wpływu M-estymatora[edytuj | edytuj kod]

Można pokazać, że funkcja wpływu M-estymatora T jest proporcjonalna do \psi (Huber, 1981 (i 2004), str. 45), co oznacza, że właściwości takiego estymatora (takie jak jego punkt załamania, sumaryczną wrażliwość i wrażliwość na lokalne przesunięcie) można wyprowadzić znając jego funkcję \psi:

IF(x;T,F) = M^{-1}\psi(x,T(F))\;

z p\times p danym przez:

M = -\int_{\mathcal{X}}\left(\frac{\partial \psi(x,\theta)}{\partial \theta}\right)_{T(F)}dF(x).

Wybór \psi i \rho[edytuj | edytuj kod]

W wielu praktycznych sytuacjach wybór funkcji \psi nie jest krytycznym warunkiem uzyskania dobrej odpornej estymaty, i wiele różnych funkcji da podobny odporny rezultat, znacząco lepszy w obecności obserwacji ostających od estymacji metodami klasycznymi (Huber, 1981).

Parametryczne estmatory odporne[edytuj | edytuj kod]

M-estymatory niekoniecznie są związane z funkcją gęstości i tym samym nie są rozwiązaniem w pełni parametrycznym. W pełni parametryczne podejście do odpornego modelowania i wnioskowania statystycznego, zarówno bayesowskiego, jak i największej wiarygodności, zwykle polega na zastosowaniu rozkładów takich jak rozkład t-Studenta z "ciężkimi ogonami".

Dla rozkładu t Studenta z \nu stopniami swobody, zachodzi:

\psi(x) = \frac{x}{x^2 + \nu}.

Dla \nu=1, rozkład t sprowadza się do rozkładu Cauchy'ego. Liczba stopni swobody jest czasem nazywana parametrem kurtozy. Ten właśnie parametr steruje wielkością "ogonów" rozkładu. Generalnie, \nu może być estymowany z próby jak każdy inny parametr. W praktyce jest to utrudnione, gdyż występuje wiele lokalnych maksimów gdy zezwolimy na zmienność \nu. Tak więc na ogół ustala się ten parametr sztywno na wartość w okolicach 4 lub 6. Wykres poniżej pokazuje, funkcję \psi dla 4 różnych wartości \nu.

TDistPsi.png

Przykład: prędkość światła[edytuj | edytuj kod]

Dla przykładu prędkości światła, uwalniając parametr kurtozy i maksymalizując wiarygodność, uzyskujemy:

 \hat\mu = 27.40, \hat\sigma = 3.81, \hat\nu = 2.13.

Ustalając \nu = 4 i maksymalizując wiarygodność, uzyskamy

\hat\mu = 27.49, \hat\sigma = 4.51.

Przypadki szczególne[edytuj | edytuj kod]

Jakkolwiek ten artykuł wiąże się z ogólnymi zasadami dotyczącymi jednowymiarowych metod statystycznych, dedykowane metody odporne istnieją dla problemów regresyjnych, GLM i estymacji parametrów różnych rozkładów.

Przypisy

  1. Polska nazwa za dr. hab. Antonim L. Dawidowiczem ([1], str. 4) oraz słownikiem ISI
  2. Elżbieta Pleszczyńska. Odporne metody statystyczne ("Robust Statistics"). „Wiadomości statystyczne”. 1-12, s. 20, 1976. Główny Urząd Statystyczny. 
  3. Dane te dostępne są wraz z materiałami uzupełniającymi książki Robust Regression and Outlier Detection (Rousseeuw i Leroy, 1986) na stronie Uniwersytetu w Kolonii.
  4. Huber (1981), Maronna et al (2006)

Bibliografia[edytuj | edytuj kod]

  • Frank R. Hampel, Elvezio M. Ronchetti, Peter J. Rousseeuw, Werner A. Stahel: Robust Statistics - The Approach Based on Influence Functions. Wiley, 1986. (ponownie opublikowane w 2005 w miękich okładkach)
  • Peter. J. Huber: Robust Statistics. Wiley, 1981. (ponownie opublikowane w 2004 w miękich okładkach)
  • Peter J. Rousseeuw, Annick M. Leroy: Robust Regression and Outlier Detection. Wiley, 1987.(ponownie opublikowane w 2003 w miękich okładkach)
  • Ricardo Maronna, Doug Martin, Victor Yohai: Robust Statistics - Theory and Methods. Wiley, 2006.
  • Andrew Gelman, John B. Carlin, Hal S. Stern, Donald B. Rubin: Bayesian Data Analysis. Chapman & Hall/CRC, 2004.
  • P. J. Rousseeuw, C. Croux. Alternatives to the Median Absolute Deviation. „Journal of the American Statistical Association”, s. 88, 1993. 

Zobacz też[edytuj | edytuj kod]

Linki zewnętrzne[edytuj | edytuj kod]