Statystyka odpornościowa

Statystyka odpornościowa^[1] lub odporne metody statystyczne^[2] – 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

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 ujemnym wzroście lub ujemną wielkość galaktyki. Naturalnie przybliżenie to jest często tak dobre, że mimo to opłaca się je 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

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

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ą bootstrapowy rozkład ś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).

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

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.

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

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

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

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 ${\displaystyle n}$ niezależnych zmiennych losowych ${\displaystyle (X_{1},\dots ,X_{n})\sim {\mathcal {N}}(0,1)}$ i odpowiednie ich realizacje ${\displaystyle x_{1},\dots ,x_{n},}$ możemy użyć wzoru ${\displaystyle {\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ą ${\displaystyle {\overline {x}}}$ dowolnie dużą za pomocą zmiany jednej wartości ${\displaystyle 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

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

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:

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

Na przykład:

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

Definicja empirycznej funkcji wpływu:
Niech ${\displaystyle n\in \mathbb {N} ^{*}}$ i ${\displaystyle X_{1},\dots ,X_{n}:(\Omega ,{\mathcal {A}})\to ({\mathcal {X}},\Sigma )}$ są niezależnymi zmiennymi losowymi o identycznym rozkładzie, a ${\displaystyle (x_{1},\dots ,x_{n})}$ jest próbą statystyczną z tych zmiennych. ${\displaystyle T_{n}:({\mathcal {X}}^{n},\Sigma ^{n})\to (\Gamma ,S)}$ jest estymatorem. Niech ${\displaystyle i\in \{1,\dots ,n\}.}$ Empiryczna funkcja wpływu ${\displaystyle EIF_{i}}$ obserwacji ${\displaystyle i}$ jest zdefiniowana jako:

${\displaystyle EIF_{i}:x\in {\mathcal {X}}\mapsto T_{n}(x_{1},\dots ,x_{i-1},x,x_{i+1},\dots ,x_{n})\in \Gamma }$

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

Funkcja wpływu i krzywa wrażliwości

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 ${\displaystyle T}$ będzie asymptotycznie nieobciążonym estymatorem pewnego parametru rozkładu ${\displaystyle F}$ spełniającego założenia modelu.

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

Wpływ na estymator ${\displaystyle T}$ przejścia z ${\displaystyle F}$ na ${\displaystyle G}$ jest wyrażony wzorem:

${\displaystyle dF_{G-F}(F)=\lim _{t\to 0^{+}}{\frac {T(tG+(1-t)F)-T(F)}{t}}}$

który jest pochodną kierunkową ${\displaystyle T}$ w punkcie ${\displaystyle F}$ w kierunku ${\displaystyle G.}$

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

${\displaystyle IF(x;T;F):=\lim _{t\to 0^{+}}{\frac {T(t\Delta _{x}+(1-t)F)-T(F)}{t}}.}$

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

Pożądane właściwości

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

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

Punkt odrzucenia

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

Maksymalna wrażliwość

${\displaystyle \gamma ^{*}(T;F):=\sup _{x\in {\mathcal {X}}}|IF(x;T;F)|}$

Wrażliwość na lokalne przesunięcia

${\displaystyle \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 ${\displaystyle x}$ do pobliskiego punktu ${\displaystyle y,}$ czyli innymi słowy dodania obserwacji w punkcie ${\displaystyle y}$ przy jednoczesnym usunięciu z ${\displaystyle x.}$

M-estymatory

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 ${\displaystyle \prod _{i=1}^{n}f(x_{i}),}$ czyli (równoważnie) minimalizacji ${\displaystyle \sum _{i=1}^{n}-\log f(x_{i}).}$ W 1964 Huber zaproponował uogólnienie tego podejścia poprzez minimalizację wyrażenia ${\displaystyle \sum _{i=1}^{n}\rho (x_{i}),}$ gdzie ${\displaystyle \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 ${\displaystyle \sum _{i=1}^{n}\rho (x_{i})}$ często może być wykonana przez różniczkowanie ${\displaystyle \rho }$ i rozwiązanie równania ${\displaystyle \sum _{i=1}^{n}\psi (x_{i})=0,}$ gdzie ${\displaystyle \psi (x)={\frac {d\rho (x)}{dx}},}$ oczywiście pod warunkiem, że ${\displaystyle \rho }$ jest różniczkowalne.

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

Funkcja ${\displaystyle \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 ${\displaystyle \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 Tukeya (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.

Właściwości M-estymatorów

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, 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

Można pokazać, że funkcja wpływu M-estymatora ${\displaystyle T}$ jest proporcjonalna do ${\displaystyle \psi }$ (Huber, 1981 (i 2004), s. 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ę ${\displaystyle \psi {:}}$

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

z ${\displaystyle p\times p}$ danym przez:

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

Wybór ${\displaystyle \psi }$ i ${\displaystyle \rho }$

W wielu praktycznych sytuacjach wybór funkcji ${\displaystyle \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 estymatory odporne

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 ${\displaystyle \nu }$ stopniami swobody, zachodzi:

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

Dla ${\displaystyle \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, ${\displaystyle \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ść ${\displaystyle \nu .}$ Tak więc na ogół ustala się ten parametr sztywno na wartość w okolicach 4 lub 6. Wykres poniżej pokazuje, funkcję ${\displaystyle \psi }$ dla 4 różnych wartości ${\displaystyle \nu .}$

Przykład: prędkość światła

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

${\displaystyle {\hat {\mu }}=27{,}40,{\hat {\sigma }}=3{,}81,{\hat {\nu }}=2{,}13.}$

Ustalając ${\displaystyle \nu =4}$ i maksymalizując wiarygodność, uzyskamy

${\displaystyle {\hat {\mu }}=27{,}49,{\hat {\sigma }}=4{,}51.}$

Przypadki szczególne

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.

Zobacz też

• metody rangowe
• statystyka gradacyjna
• statystyka nieparametryczna

Przypisy

Bibliografia

• 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. 

Linki zewnętrzne

• Justyna Majewska: Estymatory odporne zmienności w modelu Blacka-Scholesa. [dostęp 2008-05-23]. [zarchiwizowane z tego adresu (2007-07-08)].
• Olfa Nasraoui: A Brief Overview of Robust Statistics (ang.). [dostęp 2008-05-23]. [zarchiwizowane z tego adresu (2008-04-30)].
• B.D. Ripley: Robust Statistics (ang.). [dostęp 2008-05-23]. [zarchiwizowane z tego adresu (2012-10-21)].