Ten wpis jest wprowadzeniem / rozpoznaniem gruntu dla pewnej ciekawej inicjatywy wydawniczej, o której napiszę w czwartek.

Tekst ten jest bardzo techniczny, przeznaczony dla osób znających postawy teorii testowania hipotez.

Pod koniec stycznia miałem przyjemność prowadzić szkolenie z testów statystycznych ,,na poważnie”. Odbiorcami była wymagająca grupa dociekliwych, pracujących w KNFie. Część szkolenia dotyczyła testów normalności, czyli, moim zdaniem, metodologicznie bardzo ciekawego zagadnienia. Podczas przerwy zimowej uzupełniłem tę część materiałów, tak by można je było udostępnić na blogu. W uzupełnianiu pomógł mi bardzo Krzysztof Trajkowski sugerując różne rozszerzenia, za co bardzo dziękuję.

Aktualna wersja w postaci PDF znajduje się tutaj [[Wersja PDF]], a po konwersji do html i drobnych zmianach jest też przedstawiona poniżej. Wersja html była generowana z użyciem narzędzia 'latex2html’ co samo w sobie było pewnym wyzwaniem, lepiej złożona moim zdaniem jest wersja pdf.

Miłej lektury!

Wybrane testy normalności

Założenie o rozkładzie normalnym zmiennych jest podstawą wielu klasycznych technik modelowania (przykładowo analiza wariancji, regresja gaussowska, modele mieszane). Stąd też popularność technik testowania tego założenia.

W tym dokumencie przedstawiamy czternaście przykładowych testów normalności podzielonych na pięć grup, różniących się sposobem konstrukcji statystyki testowej. Poniżej przedstawiono testy, które mnie osobiście z różnych powodów interesowały, ale to nie są wszystkie testy do badania normalności. Moim celem było pokazanie mnogości podejść do tematu badania zgodności z rodziną rozkładów normalnych a nie wykonanie pełnego przeglądu testów normalności.

Ten dokument był częścią materiałów przygotowanych w ramach szkolenia o testach dla KNFu, a po uzupełnieniu i za zgodą został przekazany fundacji SmarterPoland.pl. Za pomoc w uzupełnieniu i rozbudowaniu pierwszej wersji dziękuję Krzysztofowi Trajkowskiemu, który jak zwykle wskazał wiele punktów, które warto poprawić i testów o których warto napisać.

Mam nadzieję, że dokument okaże się przydatnym. Wszelkie komentarze, sugestie propozycje do tego dokumentu proszę zgłaszać na adres Przemyslaw.Biecek@gmail.com. Z dokumentu można korzystać na zasadzie uznania autorstwa, cytując:

Przemysław Biecek, Wybrane testy normalności, 2013. Materiały Fundacji SmarterPoland.pl.

Wprowadzenie

Interesuje nas testowanie zgodności z rodziną rozkładów normalnych. Obserwujemy próbę prostą obserwacji , , niezależnych realizacji zmiennej losowej o rozkładzie . Interesuje nas weryfikacja hipotezy

$\begin{displaymath} H_0: \exists_{\mu, \sigma^2} F = \mathcal N(\mu, \sigma^2), \end{displaymath}$

przeciwko hipotezie

$\begin{displaymath} H_1: \neg\exists_{\mu, \sigma^2} F = \mathcal N(\mu, \sigma^2). \end{displaymath}$

Ponieważ parametrów $\mu$ i $\sigma^2$ najczęściej nie znamy, więc testowana hipoteza zerowa jest hipotezą złożoną, co nastręcza pewnych trudności. W poniżej opisanych testach najczęściej te parametry są estymowane z próby, z użyciem być może niestandardowych estymatorów, a statystyka testowa odpowiada badaniu zgodności z hipotezą prostą

$\begin{displaymath} H_0: F = \mathcal N(\hat \mu, \hat\sigma^2). \end{displaymath}$

Taka zamiana hipotezy złożonej na prostą ma wpływ na rozkłady estymatorów statystyk testowych.

Poniżej przedstawiamy różne testy do weryfikacji hipotezy o zgodności z rodziną rozkładów normalnych, ale wiele z prezentowanych metod może być zaadoptowanych do zagadnienia testowania zgodności z innymi rodzinami rozkładów.

Dla większości testów przedstawiamy też przykładowe alternatywy, dla których dany test ma największą moc. Wybrane scenariusze alternatyw to jednak kropla w zbiorze możliwych odstępstw od normalności.

Naturalnym pytaniem jest ,,którego testu używać”, poniższe przykłady wykazują, że nie ma testu dobrego w każdym przypadku, ale o dwóch warto wspomnieć. Testem często stosowanym jest test Shapiro Wilka, popularny i zaimplementowany w większości pakietów statystycznych. Testem o szerokim spektrum wykrywania możliwych odstępstw jest test gładki test Neymana zaimplementowany w pakiecie R.

Testy zgodności oparte o dystrybuantę empiryczną

Testy wymienione w tym rozdziale badają zgodność z rozkładem normalnym przez ocenę odległości pomiędzy dystrybuantą empiryczną a dystrybuantą rozkładu normalnego.

Opisane w tym rozdziale testy określa się terminem ,,distribution free”, co oznacza, że rozkład statystyki testowej przy prawdziwej hipotezie zerowej nie zależy od rozkładu, z którym badana jest zgodność.

Test Cramera von Misesa

Ten test jest oparty o odległość Cramera von Misesa pomiędzy dystrybuantami empiryczną i teoretyczną

$\begin{displaymath}d(F_n, F) = n \int\limits_{-\infty}^\infty (F_n(x) - F(x))^2 \, dF(x).\end{displaymath}$

(1)

Statystykę testową dla próby zbudowaną na powyższej odległości można zapisać jako

$\begin{displaymath}\begin{array}{lll}T &=& \frac{1}{12n} + \sum_{i=1}^n \left[ \frac{2i-1}{2n}-F(x_i) \right]^2.\end{array}\end{displaymath}$

(2)

W przypadku badania zgodności z rozkładem normalnym dystrybuanta jest zastępowana przez $\Phi(\hat \mu, \hat\sigma^2)$ dla oszacowanych parametrów średniej i wariancji. Rozkład statystyki przybliża się po uprzedniej modyfikacji [Stephens1986]

$\begin{displaymath}T^{*}=T(1.0 + 0.5/n).\end{displaymath}$

(3)

Przyglądając się odległości widzimy, że ,,zauważane” są różnice w centrum rozkładu, tam gdzie jest duża. Spodziewamy się ,,czułości” testu dla alternatyw zniekształconych w środku rozkładu.

W programie R ten test jest zaimplementowany w funkcji cvm.test{nortest}.

set.seed(1313); x <- rnorm(100)
library(nortest)
cvm.test(x)
##	Cramer-von Mises normality test
##data:  x 
#W = 0.0607, p-value = 0.3665
#cvm.test(x)$p.value
# 0.3664784

**Figure:**Moc vs. alternatywa Legendre’a rzędu 4 (opis tej klasy alternatyw jest w rozdziale 5), n=100, $\alpha=0.05$ . Na górnym panelu gęstość i dystrybuantę rozkładu dla którego badano moc (na czerwono) porównano z gęstością i dystrybuantą rozkładu normalnego (na czarno). Na dolnym panelu przedstawiono moc wybranych testów. Na czerwono zaznaczono wyniki nie gorsze niż 80% mocy najlepszego z testów (patrz pionowa czerwona linia). Mniej klasyczna alternatywa, ale też przykład w którym najlepiej wypadł test Cramera von Misesa. Moc oszacowana na bazie 10 000 powtórzeń.

Test Andersona-Darlinga

Ten test jest oparty o ważoną odległość Cramera von Misesa pomiędzy dystrybuantami empiryczną

i teoretyczną

z wagami odpowiadającymi odwrotności wariancji dystrybuanty empirycznej [Anderson1954] (pamiętamy, że

ma rozkład dwumianowy).

$\begin{displaymath}d(F_n, F) = n \int\limits_{-\infty}^\infty \frac{(F_n(x) - F(x))^2 }{[F(x)\; (1-F(x))]} \, dF(x).\end{displaymath}$

(4)

Statystykę testowę opartą o powyższą odległość dla próby prostej można przedstawić jako

$\begin{displaymath}\begin{array}{lll}S &=& \sum_{k=1}^n \frac{2k-1}{n}\left[\l......1-F(x_{n+1-k})\right)\right], \\A^2 &=& -n-S. \\\end{array}\end{displaymath}$

(5)

Podobnie jak w teście Cramera von Misesa, testując zgodność z rozkładem normalnym przyjmuje się $F = \Phi(\hat \mu, \hat \sigma^2)$ , a rozkład statystyki testowej przybliża się po po uprzedniej korekcie [Stephens1986] (w literaturze spotyka się też inne korekty, poniżej prezentowana jest używana w programie R)

$\begin{displaymath}A^{*,2}=A^2(1.0 + 0.75/n +2.25/n^{2}).\end{displaymath}$

(6)

W programie R ten test jest zaimplementowany w funkcji ad.test{nortest}.

set.seed(1313); x <- rnorm(100)
library(nortest)
ad.test(x)
##	Anderson-Darling normality test
##data:  x 
#A = 0.4141, p-value = 0.3301
#ad.test(x)$p.value
# [1] 0.3300575

**Figure:**Moc vs. alternatywa Legendre’a rzędu 6 (opis tej klasy alternatyw jest w rozdziale 5), n=250, $\alpha=0.05$ .Na górnym panelu gęstość i dystrybuantę rozkładu dla którego badano moc (na czerwono) porównano z gęstością i dystrybuantą rozkładu normalnego (na czarno). Na dolnym panelu przedstawiono moc wybranych testów. Na czerwono zaznaczono wyniki nie gorsze niż 80% mocy najlepszego z testów (patrz pionowa czerwona linia). Z wykryciem tej alternatywy najlepiej poradził sobie test Andersona-Darlinga, zapewnie przez cięższy prawy i lżejszy lewy ogon. Moc oszacowana na bazie 10 000 powtórzeń.

Test Kołmogorowa-Smirnowa

Test oparty o odległość supremum pomiędzy dystrybuantami empiryczną i teoretyczną

$\begin{displaymath}d(F_n, F) = \sup_x\vert F_n (x) - F(x)\vert.\end{displaymath}$

(7)

Statystyka testowa oparta o powyższą odległość sprowadza się do liczenia maksimum modułu różnicy dystrybuant w punktach skoku dystrybuanty empirycznej

$\begin{displaymath}D = \max_{x_i}\vert F_n (x_i) - F(x_i)\vert.\end{displaymath}$

(8)

Rozważa się też ważoną wersją testu Kołmogorowa Smirnowa, bardziej czułą na różnice w ogonach

$\begin{displaymath}D' = \max_{x_i}\frac{\vert F_n (x_i) - F (x_i)\vert}{\sqrt{ F(x_i) (1 - F(x_i)) }}.\end{displaymath}$

Rozkład statystyki testowej można wyznaczyć w sposób dokładny (rekurencyjnie) dla prostej hipotezy zerowej, a więc dla porównania z jednym określonym rozkładem. Asymptotycznie, ta statystyka przemnożona przez $\sqrt{n}$ ma rozkład Kołmogorowa.Test ten pomimo łatwego opisu probabilistycznego nie jest stosowany z uwagi na moc niższą niż konkurencja.

W programie R ten test jest zaimplementowany w funkcji ks.test{stats}.

set.seed(1313)
x <- rnorm(100)
ks.test(scale(x), "pnorm")
##	One-sample Kolmogorov-Smirnov test
##data:  scale(x) 
#D = 0.0543, p-value = 0.9297
#alternative hypothesis: two-sided 
#ks.test(scale(x), "pnorm")$p.value
# [1] 0.9297479

Test Lillieforsa

Ten test to modyfikacja testu Kołmogorowa-Smirnowa zaproponowana przez Huberta Lillieforsa [Lilliefors1967], pozwalająca na testowanie zgodności z całą rodziną rozkładów normalnych, bez znajomości parametrów średniej i odchylenia standardowego (test Kołmogorowa-Smirnowa pozwala na zbadanie zgodności z jednym określonym rozkładem).

Statystyka testowa w przypadku testu Lillieforsa wygląda tak samo jak w przypadku testu Kołmogorowa Smirnowa. Różnica polega na zastosowaniu innego rozkładu dla statystyki testowej (przybliżenie rozkładu dokładnego), uwzględniającego to, że hipoteza zerowa jest hipotezą złożoną.

W programie R ten test jest zaimplementowany w funkcji lillie.test{nortest} i lillieTest{fBasics}.

set.seed(1313); x <- rnorm(100)
lillie.test(x)$p.value
# [1] 0.6643097
lillieTest(x)@test$p.value
# 0.6643097

**Figure:**Moc vs. alternatywa w której część rozkładu skupiono w punkcie, n=100, $\alpha=0.05$ .Na górnym panelu gęstość i dystrybuantę rozkładu dla którego badano moc (na czerwono) porównano z gęstością i dystrybuantą rozkładu normalnego (na czarno). Na dolnym panelu przedstawiono moc wybranych testów. Na czerwono zaznaczono wyniki nie gorsze niż 80% mocy najlepszego z testów (patrz pionowa czerwona linia). Z wykryciem tej alternatywy najlepiej poradził sobie test Lillieforsa. Moc oszacowana na bazie 10 000 powtórzeń.

Testy oparte o momenty z próby

Testy oparte o dystrybuantę empiryczną można było stosować do badania zgodności z dowolnym rozkładem. Poniższe testy wykorzystują specyficzną właściwość rozkładu normalnego, tzn. zerową skośność i kurtozę równą 3 (licząc zgodnie z definicją ) lub kurtozę równą 0 (licząc zgodnie z definicją ).

Aby łatwiej nam było opisać poniższe testy, przyjmijmy następujące oznaczenia na to k-ty centralny moment z próby, skośność i kurtozę.

$\begin{displaymath}\begin{array}{lll}m_k &= &\frac 1n \sum_i (x_i-\bar x)^k,\\......{b_1} &=& m_3/m_2^{3/2}, \\b_2 &=& m_4/m_2^2. \\\end{array}\end{displaymath}$

(9)

Test D’Agostino oparty o skośność

Test D’Agostino (1970) za statystykę testową wykorzystuje przekształconą próbkową skośność.W celu ,,normalizacji” rozkładu statystyki testowej stosuje się następujące przekształcenie

$\begin{displaymath}\begin{array}{lll}Y &=& \sqrt{b_1}\left(\frac{(n+1)(n+3)}{6......pha + \left((Y/\alpha)^2+1\right)^{1/2}\right). \\\end{array}\end{displaymath}$

(10)

Statystyką testową jest $Z(\sqrt{b_1})$ . Przy założeniu prawdziwej hipotezy zerowej ta statystyka ma asymptotyczny rozkład normalny [Agostino1970].

W programie R ten test jest zaimplementowany w funkcji dagoTest{fBasics}.

set.seed(1313); x <- rnorm(100)
dagoTest(x)@test$p.value[2]
# Skewness Test 
#     0.9287743

**Figure:**Moc vs. rozkład Gumbela (wartości ekstremalnej), n=60, $\alpha=0.05$ . Na górnym panelu gęstość i dystrybuantę rozkładu dla którego badano moc (na czerwono) porównano z gęstością i dystrybuantą rozkładu normalnego (na czarno). Na dolnym panelu przedstawiono moc wybranych testów. Na czerwono zaznaczono wyniki nie gorsze niż 80% mocy najlepszego z testów (patrz pionowa czerwona linia). Ciężki prawy ogon rozkładu wartości ekstremalnej najlepiej wykrywa test skośności D’Agostino. Moc oszacowana na bazie 10 000 powtórzeń.

Test D’Agostino oparty o kurtozę

Test D’Agostino za statystykę testową wykorzystuje przekształconą próbkową kurtozę.

Bazując na wynikach [Anscombe1983] w celu ,,normalizacji” rozkładu oceny kurtozy stosuje się następujące przekształcenie

$\begin{displaymath}\begin{array}{lll}\mathsf{E}\{b_2\} &=& 3(n-1)/(n+1), \\\......{2/(A-4)}}\right)^{1/3}\right)/\sqrt{\frac{2}{9A}},\end{array}\end{displaymath}$

(11)

gdzie $\mathsf{E}\{\cdot\}$ , $\mathsf{var}\{\cdot\}$ to obciążenie i wariancja oceny kurtozy dla rozkładu normalnego.

Statystyką testową jest $Z(\sqrt{b_2})$ . Przy założeniu prawdziwej hipotezy zerowej ta statystyka ma asymptotyczny rozkład normalny.

W programie R ten test jest zaimplementowany w funkcji dagoTest{fBasics}.

set.seed(1313); x <- rnorm(100)
dagoTest(x)@test$p.value[3]
#Kurtosis Test 
#    0.4417257

**Figure:**Moc vs. rozkład mieszaniny , , n=500, $\alpha=0.05$ . Na górnym panelu gęstość i dystrybuantę rozkładu dla którego badano moc (na czerwono) porównano z gęstością i dystrybuantą rozkładu normalnego (na czarno). Na dolnym panelu przedstawiono moc wybranych testów. Na czerwono zaznaczono wyniki nie gorsze niż 80% mocy najlepszego z testów (patrz pionowa czerwona linia). Z wykryciem spłaszczenia rozkładu najlepiej poradził sobie test kurtozy D’Agostino. Moc oszacowana na bazie 10 000 powtórzeń.

Test typu omnibus D’Agostino-Pearsona oparty o kurtozę i skośność

Łącząc dwa powyższe testy otrzymuje się test czuły na odstępstwa od normalności zarówno w postaci niezerowej skośności jak i kurtozy istotnie różniej od 3 [Agostino1973].

Statystyką testową jest

$\begin{displaymath}K^2 = \left(Z\left(\sqrt{b_1}\right)\right)^2 + \left(Z\left(b_2\right)\right)^2,\end{displaymath}$

(12)

gdzie $Z\left(\sqrt{b_1}\right)$ to statystyka testowa testu opartego o skośność a $Z\left(b_2\right)$ to statystyka testowa testu opartego o kurtozę.

Asymptotyczny rozkład tej statystyki to oczywiście rozkład $\chi^2_2$ .

W programie R ten test jest zaimplementowany w funkcji dagoTest{fBasics}.

set.seed(1313); x <- rnorm(100)
dagoTest(x)@test$p.value[1]
#Omnibus  Test 
#    0.7408977

Rozszerzeniem testu D’Agostino-Pearsona na przypadek testowania zgodności z wielowymiarowym rozkładem normalnym, dla statystyki testowej opartej o kurtozę i skośność, jest test Doornika-Hansena [Doornik1994]. Jest on zaimplementowany w funkcji normality.test1{normwhn.test}, a ta implementacja pozwala też na testowanie kierunkowych alternatyw dla skośności i kurtozy, czyli że rozkład jest bardziej/mniej lewoskośny/spłaszczony niż rozkład normalny. W funkcji normality.test2{normwhn.test} zaimplementowana jest modyfikacja pozwalająca na testowanie zgodności z rozkładem normalnym przy osłabionym założeniu dotyczącym niezależności. Tzn. jest to test normalności dopuszczający słabą strukturę zależności pomiędzy obserwacjami, patrz [Velasco2004].

Więcej o teście Doornika-Hansena i innych wielowymiarowych testach normalności (np. teście Mardii, wielowymiarowym Shapiro-Francia) przeczytać można w [BiecekTrajkowski2011].

Test typu omnibus Jarque-Bera oparty o kurtozę i skośność

Innym testem opartym o kurtozę i skośność jest test Jarque-Bera [Jarque1987]. Statystyka testowa w przypadku tego testu ma łatwiejszą postać niż dla testu D’Agostino-Pearsona. Traci się jednak na niedokładnym oszacowaniu wartości krytycznych przy niewielkich wielkościach próby. Asymptotycznie ten test jest tak samo mocny jak test D’Agostino-Pearsona, ale na asymptotykę można liczyć jedynie w przypadku dużych prób.

Statystyka testowa ma postać

$\begin{displaymath}\text{\textsc{JB}}= \frac n6 \left(\left(\sqrt{b_1}\right)^2 + \frac 14 \left(b_2 - 3\right)^2\right).\end{displaymath}$

(13)

Asymptotyczny rozkład statystyki testowej $\text{\textsc{JB}}$ to rozkład $\chi^2_2$ . W funkcji jbTest{fBasics} zaimplementowane są też dwie inne metody liczenia p-wartości i szacowania kwantyli rozkładu statystyki JB oparte o mnożniki Lagrange, przy czym raportowane przez tę funkcję p-wartości są zaokrąglone do trzeciego miejsca po przecinku dziesiętnym. P-wartość wyznaczane przez funkcję jarqueberaTest{fBasics} nie jest zaokrąglana, jest więc podawana z większą dokładnością, ale jedynie dla przybliżenia rozkładu statystyki testowej rozkładem asymptotycznym. Zaokrąglanie p-wartości do tysięcznych uniemożliwia testowanie na poziomie istotności rzędu $\alpha=0.001$ lub mniejszym.

set.seed(1313); x <- rnorm(100)
jbTest(x)@test$p.value
# LM p-value ALM p-value  Asymptotic 
#      0.894       0.810       0.900

Testy oparte o statystyki pozycyjne

Kolejne dwa testy oparte są o statystyki pozycyjne z próby. Ideę ich działania najłatwiej uchwycić przyglądając się wykresowi kwantylowemu QQ-plot.

Test Shapiro-Wilka

Jednym z bardziej popularnych testów klasycznych jest test Shapiro-Wilka [Shapiro1965]. Statystyka testowa jest oparta o iloczyn skalarny unormowanych uporządkowanych statystyk porządkowych z unormowanych wartościami oczekiwanymi statystyk porządkowych. Innymi słowy

$\begin{displaymath}\begin{array}{lll}W &=& \frac{\left(\sum_{i=1}^n a_i x_{(i)...... &=& \frac{m^T V^{-1}}{(m^T V^{-1}V^{-1} m)^{1/2}},\end{array}\end{displaymath}$

(14)

gdzie to wektor wartości oczekiwanych uporządkowanych statystyk pozycyjnych w rozkładzie normalnym, a to macierz kowariancji tych statystyk testowych.

Wartości wektora i macierzy oraz rozkładu statystyki testowej pracowicie stablicowano [Sarhan1956]. W kolejnych latach proponowano różne aproksymacje tych wartości [Davis1978,Royston1982,Royston1992]. Wykorzystując tę ostatnią aproksymację wzór na statystykę testową można wyznaczyć następującym algorytmem.

W pierwszym kroku wyznaczane są wartości oczekiwane statystyk porządkowych i wektor wag

$\begin{displaymath}\begin{array}{lll}\tilde m_i &=& \Phi^{-1}\{\frac{i-3/8}{n+...... \sqrt{(\mathbf{\tilde m}^T\mathbf{\tilde m})}. \\\end{array}\end{displaymath}$

(15)

W zależności od rozmiaru próby stosuje się dwa różne przybliżenia.

Dla $4 \leq n < 12$ używa się przybliżenia

$\begin{displaymath}\begin{array}{lll}\tilde a_n &=& - 2.706056u^5 + 4.434685u^......f{\tilde m} - 2\tilde m_n^2)/(1-2\tilde a_n^2). \\\end{array}\end{displaymath}$

(16)

Dla $n \geq 12$ używa się przybliżenia

$\begin{displaymath}\begin{array}{lll}\tilde a_n &=& - 2.706056u^5 + 4.434685u^......m_{n-1}^2)/(1-2\tilde a_n^2-2\tilde a_{n-1}^2). \\\end{array}\end{displaymath}$

(17)

Wyznaczywszy te współczynniki, statystykę testową wyznaczamy ze wzoru

$\begin{displaymath}W=\left(\sum_{i=1}^n \tilde a_i x_i \right)^2 / \sum_{i=1}^n (x_i-\bar x)^2.\end{displaymath}$

(18)

Statystyka po przekształceniu przez funkcję ma rozkład bliski rozkładowi normalnemu o średniej $\mu$ i odchyleniu standardowym $\sigma$ [Royston1993].

Dla $4 \leq n < 12$ używa się przybliżenia

$\begin{displaymath}\begin{array}{lll}g(W) &=& -\ln\left(\gamma-\ln(1-W)\right)...... 0.062767n^2 -0.77857n + 1.3822\},\\u &=& \ln(n).\end{array}\end{displaymath}$

(19)

Dla $n \geq 12$ używa się przybliżenia

$\begin{displaymath}\begin{array}{lll}g(W) &=& \ln(1-W), \\\mu &=& 0.0038915u......0.0030302u^2 -0.082676u -0.4803\},\\u &=& \ln(n).\end{array}\end{displaymath}$

(20)

Statystyka $Z=\left(g(W)-\mu\right)/\sigma$ ma w przybliżeniu rozkład $\mathcal N(0,1)$ .

W programie R ten test jest zaimplementowany w funkcji shapiro.test{stats}. Ta implementacja pozwala na testowanie normalności wektorów o długości od 3 do 5000 obserwacji.

1
2
3

set.seed(1313); x <- rnorm(100)
shapiro.test(x)$p.value
# [1] 0.4879372

**Figure:**Moc vs. rozkład lognormalny, n=20, $\alpha=0.01$ .Na górnym panelu gęstość i dystrybuantę rozkładu dla którego badano moc (na czerwono) porównano z gęstością i dystrybuantą rozkładu normalnego (na czarno). Na dolnym panelu przedstawiono moc wybranych testów. Na czerwono zaznaczono wyniki nie gorsze niż 80% mocy najlepszego z testów (patrz pionowa czerwona linia). Skośność i ciężkie ogony rozkładu log normlanego najlepiej wykrywa test Shapiro-Wilka. Moc oszacowana na bazie 10 000 powtórzeń.

Test Shapiro-Francia

Test Shapiro-Francia jest uproszczoną wersją testu Shapiro-Wilka, w którym macierz kowariancji statystyk pozycyjnych zastąpiono przez macierz identycznościową [Shapiro1972].Dla dużych prób ten test ma podobne zachowanie jak test Shapiro Wilka.

Statystyką testową jest

$\begin{displaymath}\begin{array}{lll}W' &=& \frac{\left(\sum_{i=1}^n a_i x_{(i...... \\(a_1, ..., a_n) &=& \frac{m^T}{(m^T m)^{1/2}}.\end{array}\end{displaymath}$

(21)

Warto zauważyć, że statystyka to korelacja punktów na wykresie kwantylowym QQ-plot.Do obliczeń wykorzystuje się wzory

$\begin{displaymath}\begin{array}{lll}\tilde m_i &=& \Phi^{-1}\{(i-3/8)/(n+1/4)......ilde a_i x_i \right)^2 / \sum_i (x_i-\bar x)^2. \\\end{array}\end{displaymath}$

(22)

Podobnie jak w teście Shapiro-Wilka statystyka $g(W')=\ln(1-W')$ jest przybliżana rozkładem normalnym o parametrach

$\begin{displaymath}\begin{array}{lll}\mu &=& 1.0521u - 1.2725, \\u &=& \ln\l......030, \\v &=& \ln\left(\ln(n)\right)+2/\ln(n), \\\end{array}\end{displaymath}$

(23)

dla prób z przedziału 5 do 5000 obserwacji (w R nie można testować tym testem większych prób).

Statystyka $Z = \left(g(W^T)-\mu\right)/\sigma$ ma w przybliżeniu rozkład normalny $\mathcal N(0,1)$ .

W programie R ten test jest zaimplementowany w funkcji sf.test{nortest}. Ta implementacja pozwala na testowanie normalności wektorów o długości od 5 do 5000 obserwacji.

1
2
3

set.seed(1313); x <- rnorm(100)
sf.test(x)$p.value
# [1] 0.3414609

**Figure:**Moc vs. rozkład Cauchy’ego, n=20, $\alpha=0.01$ .Na górnym panelu gęstość i dystrybuantę rozkładu dla którego badano moc (na czerwono) porównano z gęstością i dystrybuantą rozkładu normalnego (na czarno). Na dolnym panelu przedstawiono moc wybranych testów. Na czerwono zaznaczono wyniki nie gorsze niż 80% mocy najlepszego z testów (patrz pionowa czerwona linia). Ciężkie ogony rozkładu Cauchego łatwo wykryć, ale najlepiej sobie z tym zadaniem poradził test Shapiro-Francia. Moc oszacowana na bazie 10 000 powtórzeń.

Testy oparte o statystykę $\chi^2$

Dla rozkładów o dyskretnym nośniku do badania zgodności można wykorzystać test zgodności $\chi^2$ . W postaci ogólnej jest on zaimplementowany w programie R w funkcji chisq.test(stats).

W przypadku testowania zgodności z rozkładem ciągłym, można też wykorzystać ten test, dzieląc uprzednio nośnik rozkładu na przedziałów, domyślnie w R liczbę przedziałów wyznacza się z wzoru $k = 2 n^{2/5}$ . Następnie bada się zgodność z takim zdyskretyzowanym rozkładem.

Test Pearsona

Test Pearsona jest oparty o statystykę

$\begin{displaymath}X^2 = \sum_{i=1}^{k} \frac{(O_i - E_i)^2}{E_i},\end{displaymath}$

która dla prawdziwej hipotezy zerowej ma asymptotyczny rozkład $\chi^2$ . Symbol

oznacza oczekiwaną liczebność obserwacji w klasie

, a

oznacza obserwowaną liczebność obserwacji w klasie

. Przybliżenie rozkładem asymptotycznym ma sens gdy liczebności

nie są bardzo małe, zaleca się przynajmniej 10 obserwacji w klasie, choć niektóre źródła podają za wymóg 5 lub 20 obserwacji.

Do przybliżenia rozkładu statystyki testowej używa się rozkładu $\chi^2_{k-1}$ lub $\chi^2_{k-3}$ . W przypadku nieznanych parametrów $\mu$ i $\sigma$ żaden z tych rozkładów nie jest dokładnym rozkładem statystyki testowej ale zaleca się wyznaczanie p-wartości z użyciem rozkładu $\chi^2_{k-3}$ .

W programie R ten test jest zaimplementowany w funkcji pearson.test{nortest} oraz pchiTest{fBasics}. W pierwszej z tych funkcji za wybór rozkładu asymptotycznego odpowiada opcja adjust=TRUE ( $\chi^2_{k-3}$ , domyślnie) lub adjust=FALSE ( $\chi^2_{k-1}$ ). W drugiej z tych funkcji liczone są dwie p-wartości, dla rozkładu $\chi^2_{k-1}$ i dla $\chi^2_{k-3}$ .

set.seed(1313); x <- rnorm(100)
pearson.test(x)$p.value
# [1] 0.5878833pchiTest(x)@test$p.value
#    Adjusted Not adjusted 
#   0.5878833    0.7515082

**Figure:**Moc vs. alternatywa Legendre’a rzędu 7 (opis tej klasy alternatyw jest w rozdziale 5), n=100, $\alpha=0.05$ .Na górnym panelu gęstość i dystrybuantę rozkładu dla którego badano moc (na czerwono) porównano z gęstością i dystrybuantą rozkładu normalnego (na czarno). Na dolnym panelu przedstawiono moc wybranych testów. Na czerwono zaznaczono wyniki nie gorsze niż 80% mocy najlepszego z testów (patrz pionowa czerwona linia). Z wykryciem wielomodalności najlepiej poradził sobie test Pearsona. Moc oszacowana na bazie 10 000 powtórzeń. Test ddst był wywołany z domyślnym ograniczeniem na liczbę statstyk , stąd też niska moc w wykrywaniu takich złożonych odstępstw, gdyby pozwolić temu testowi na eksplorację większych to prawdopodobnie miałby zdecydowanie większą moc.

Testy ,,data-driven”

Poza wymienionymi klasycznymi testami zgodności, w literaturze znaleźć można też ,,nowoczesne” testy w znacznym stopniu oparte na możliwościach obliczeniowych współczesnych komputerów. Poniżej przedstawiona będzie modyfikacja gładkiego testu Neymana, ciekawa metodologicznie i łatwo dostępna w programie R.

Gładki test Neymana zgodności z rozkładem normalnym

Poniższy opis statystyki testowej dotyczy badania zgodności z rodzinami rozkładów z parametrem położenia i skali, a w szczególności z rodziną rozkładów normalnych [Janic2008]. Więcej o gładkich testach znaleźć można w [Kallenberg1997].

Statystyką testową jest

$\begin{displaymath}\begin{array}{lll}W_k &=& \left[\frac 1{\sqrt n} \sum_{i=1}......t[\frac 1{\sqrt n} \sum_{i=1}^n \ell(x_i)\right]^T,\end{array}\end{displaymath}$

(24)

gdzie $\ell(x_i)$ to

-wymiarowy wektor statystyk pomocnicznych a

to macierz kowariancji tych statystyk dla prawdziwej hipotezy zerowej, ${I}=Cov_{\theta_0}[\ell(x_1)]^T[\ell(x_1)]$ . Aby pracować z jednostkową macierzą kowariancji, za statystyki wybierane są $\ell(x_i)=(\phi_1(F(x_i)),...,\phi_k(F(x_i)))$ , gdzie

jest dystrybuantą rozkładu z którym badamy zgodność a $\phi_i()$ to ortonormalne wycentrowane wielomiany na odcinku

. Dla hipotez złożonych i badania zgodności z rodzinami rozkładów, bada się zachowanie statystyk $\ell(x_i; \hat\gamma)$ wyznaczonych dla rozkładu o parametrach $\hat \gamma$ .

Wybór statystyk i ich liczba wpływa na to, na jakie alternatywy ten ten będzie ,,wyczulony”. W poniżej przedstawionej implementacji, liczba jest wybierana z użyciem pewnego kryterium wyboru modelu, a $\phi$ to pierwsze funkcji z bazy funkcji ortonormalnych, np. bazy kosinusów lub bazy wielomianów Legendre’a. Użycie kryterium wyboru modelu do wyboru liczby powoduje, że im większa próba, tym na większą liczbę alternatyw ten test będzie czuły.

Parametry rozkładu normalnego ocenia się estymatorami

$\begin{displaymath}\begin{array}{lll}\tilde \mu &=& \frac 1n \sum_{i=1}^n x_i,......sum_{i=1}^{n-1}(x_{(i+1)}-x_{(i)})(H_{i+1}-H_i),\\\end{array}\end{displaymath}$

(25)

gdzie $x_{(i)}$ to uporządkowane rosnąco statystyki porządkowe próby, a $H_i=\Phi^{-1}((i-3/8)(n+1/4))$ .

Statystyka testowa $W_k(\tilde \mu, \tilde \sigma)$ opisana jest szczegółowo w [Janic2008]. Wartość statystyki testowej można wyznaczyć dzięki stablicowanym macierzom $[I^*(\tilde \mu, \tilde \sigma)]^{-1}$ . Do wyznaczania p-wartości stosuje się aproksymacje funkcji kwantylowej.

W programie R ten test jest zaimplementowany w funkcji ddst.norm.test{ddst}.

set.seed(1313); x <- rnorm(100)
ddst.norm.test(x, compute.p=TRUE)
#	Data Driven Smooth Test for Normality
##data:  x,   base: ddst.base.legendre,   c: 100 
#WT* = 0.0318, n. coord = 1, p-value = 0.872

**Figure:**Moc vs. alternatywa Legendre’a rzędu 2, n=100, $\alpha=0.05$ . Na górnym panelu gęstość i dystrybuantę rozkładu dla którego badano moc (na czerwono) porównano z gęstością i dystrybuantą rozkładu normalnego (na czarno). Na dolnym panelu przedstawiono moc wybranych testów. Na czerwono zaznaczono wyniki nie gorsze niż 80% mocy najlepszego z testów (patrz pionowa czerwona linia). Z rozpoznaniem zaburzenia masy w środku rozkładu, najlepiej poradził sobie gładki test Neymana. Moc oszacowana na bazie 10 000 powtórzeń.

Test oparty o statystykę E (energii)

Do studium symulacyjnego dodano test zgodności z wielowymiarowym rozkładem normalnym (w tym przypadku jednowymiarowym), test oparty o statystykę E (patrz [Szekely2005]), czyli średnią odległość pomiędzy każdą parą obserwacji (to nie jest prawdziwa odległość w matematycznym sensie, chodzi o $h(x,y) = E\vert\vert x-Z'\vert\vert + E\vert\vert y-Z'\vert\vert-E\vert\vert Z'-Z''\vert\vert-\vert\vert x-y\vert\vert$ )

$\begin{displaymath}E = -\frac 1n \sum_{j,k = 1}^n (y_j - y_k)^2 + 2 \sum_{j=1}^n \left(2 y_j \Phi(y_j) + \phi(y_j) - y_j\right) - c,\end{displaymath}$

gdzie

to standaryzowana próba

, $\Phi$ i $\phi$ to dystrybuanta i gęstość rozkładu normalnego $\mathcal N(0,1)$ a

to stała równa

gdzie $Z', Z'' \sim_{iid} \mathcal N(0,1)$ .

W programie R ten test jest zaimplementowany w funkcji mvnorm.etest{energy}.

set.seed(1313); x <- rnorm(100)
mvnorm.etest(x, R = 999)
##        Energy test of multivariate normality: estimated parameters
##data:  x, sample size 100, dimension 1, replicates 999 
#E-statistic = 0.4521, p-value = 0.3373

Test oparty o statystykę SJ

W symulacjach uwzględniono również test zgodności z rozkładem normalnym oparty o statystykę SJ (S od ,,standard deviation” i J od ,,average absolute deviation from the median”) [Gel2007].Jest to test kierunkowy badający ciężko-ogonowość rozkładu, stąd też dla pewnych alternatyw moc tego testu będzie znacznie niższa niż poziom istotności.

$\begin{displaymath}\begin{array}{lll}SJ &=& \frac{S_n}{J_n}, \\S_n &=& \frac......qrt{\pi/2}}{n} \sum_{i=1}^n \vert x_i - M\vert, \\\end{array}\end{displaymath}$

(26)

to mediana próby

a $\bar x$ to średnia z próby. Statystyka $\sqrt{n}(SJ - 1)$ ma asymptotyczny rozkład normalny.

W programie R ten test jest zaimplementowany w funkcji sj.test{lawstat}.

set.seed(1313); x <- rnorm(100)
sj.test(x)
##        Test of Normality - SJ Test
##data:  x 
#Standardized SJ Statistic = 0.0353, ratio of S to J = 1.001, p-value = 0.4865

**Figure:**Moc vs. rozkład logistyczny, n=350, $\alpha=0.05$ . Na górnym panelu gęstość i dystrybuantę rozkładu dla którego badano moc (na czerwono) porównano z gęstością i dystrybuantą rozkładu normalnego (na czarno). Na dolnym panelu przedstawiono moc wybranych testów. Na czerwono zaznaczono wyniki nie gorsze niż 80% mocy najlepszego z testów (patrz pionowa czerwona linia). Aby rozpoznać rozkład logistyczy z przyzwoitą mocą potrzeba było ponad 300 obserwacji. Najlepiej z tym zadaniem poradził sobie test oparty o statystykę SJ, który częściej wykrywał cięższe ogony niż w rozkładzie normalnym. Moc oszacowana na bazie 10 000 powtórzeń.

Bibliography

Anscombe1983

Anscombe FJ, Glynn WJ. Distribution of the Kurtosis Statistic b 2 for Normal Samples. Biometrika. 1983; 70(1):227-34. Blom G. Statistical Estimates and Transformed Beta-Variables. New York, NY. John Wiley & Sons, 1958.

Agostino1970

D’Agostino RB. Transformation to Normality of the Null Distribution of g1. Biometrika. 1970; 57(3):679-81.

Agostino1973

D’Agostino R, Pearson ES. Tests for Departure from Normality. Empirical Results for the Distributions of b 2 and b 1. Biometrika. 1973; 60(3):613-22.

Agostino1990

D’Agostino RB, Belanger A, D’Agostino Jr RB. A Suggestion for Using Powerful and Informative Tests of Normality. The American Statistician. 1990; 44(4):316-21.

Anderson1954

Anderson, T.W. and Darling, D.A. (1954). A Test of Goodness-of-Fit. Journal of the American Statistical Association 49: 765ľ769.

BiecekTrajkowski2011

Przemyslaw Biecek, Krzysztof Trajkowski (2011). Na przelaj przez Data Mining, URL http://www.biecek.pl/NaPrzelajPrzezDataMining.

Davis1978

Davis C, Stephens M. Approximating the Covariance Matrix of Normal Order Statistics. Applied Statistics. 1978; 27(2):206-212.

Dickey1979

Dickey, D.A. and W.A. Fuller (1979). Distribution of the Estimators for Autoregressive Time Series with a Unit Root. Journal of the American Statistical Association, 74, p. 427ľ431.

Doornik1994

Doornik, J.A., and H. Hansen (1994). An Omnibus Test for Univariate and Multivariate Normality, Working Paper, Nuffield College, Oxford University, U.K. Lobato, I.,

Gel2007

Gel, Y. R., Miao, W., and Gastwirth, J. L. (2007). Robust Directed Tests of Normality Against Heavy Tailed Alternatives. Computational Statistics and Data Analysis 51, 2734-2746.

Jarque1980

Jarque CM, Bera AK. Efficient tests for normality, homoscedasticity and serial independence of regression residuals. Economics Letters. 1980; 6(3):255-259.

Jarque1987

Jarque CM, Bera AK. A Test for Normality of Observations and Regression Residuals. International Statistical Review/Revue Internationale de Statistique. 1987; 55(2):163ľ172.

Janic2008

Janic, A. and Ledwina, T. (2008). Data-driven tests for a location-scale family revisited. J. Statist. Theory. Pract. Special issue on Modern Goodness of Fit Methods.

Kallenberg1997

Kallenberg, W.C.M., Ledwina, T. (1997 a). Data driven smooth tests for composite hypotheses: Comparison of powers. J. Statist. Comput. Simul. 59, 101ľ121.

Lilliefors1967

Lilliefors, H. (June 1967), On the KolmogorovľSmirnov test for normality with mean and variance unknown, Journal of the American Statistical Association, Vol. 62. pp. 399ľ402.

Royston1982

Royston JP. Algorithm AS 177: Expected normal order statistics (exact and approximate). Journal of the Royal Statistical Society. Series C (Applied Statistics). 1982; 31(2):161ľ165.

Royston1992

Royston P. Approximating the Shapiro-Wilk W-test for non-normality. Statistics and Computing. 1992; 2(3):117-119.

Royston1993

Royston P. A Toolkit for Testing for Non-Normality in Complete and Censored Samples. The Statistician. 1993; 42(1):37.

Sarhan1956

Sarhan AE, Greenberg BG. Estimation of location and scale parameters by order statistics from singly and doubly censored samples. The Annals of Mathematical Statistics. 1956; 27(2):427ľ451.

Shapiro1965

Shapiro SS, Wilk MB. An analysis of variance test for normality (complete samples). Biometrika. 1965; 52(3-4):591-611.

Shapiro1972

Shapiro SS, Francia R. An approximate analysis of variance test for normality. Journal of the American Statistical Association. 1972; 67(337):215ľ216.

Stephens1986

Stephens, M.A. (1986): Tests based on EDF statistics. In: D’Agostino, R.B. and Stephens, M.A., eds.: Goodness-of-Fit Techniques. Marcel Dekker, New York.

Szekely2005

Szekely, G. J. and Rizzo, M. L. (2005) A New Test for Multivariate Normality, Journal of Multivariate Analysis, 93/1, 58-80.

Velasco2004

C. Velasco (2004). A Simple Test of Normality of Time Series, Econometric Theory, 20, pp. 671-689, Cambridge University Press.

Dodatek, kody w programie R

Kod w programie R rysujący gęstości i dystrybuanty wybranych alternatyw

Rozkład Cauchego

x <- seq(-5,5,0.01)
plot(x, dnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="g&#281;sto&#347;&#263;")
lines(x, dcauchy(x,0,0.8), type="l", lwd=3, col="red3")
plot(x, pnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="dystrybuanta")
lines(x, pcauchy(x,0,0.8), type="l", lwd=3, col="red3")

Rozkład Gumbela (wartości ekstremalnej)

plot(x, dnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="g&#281;sto&#347;&#263;")
lines(x, dgumbel(x,-0.3), type="l", lwd=3, col="red3")
plot(x, pnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="dystrybuanta")
lines(x, pgumbel(x,-0.3), type="l", lwd=3, col="red3")

Rozkład mieszaniny dwóch rozkładów

plot(x, dnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="g&#281;sto&#347;&#263;")
lines(x, (dnorm(x,-1)+dnorm(x,1))/2, type="l", lwd=3, col="red3")
plot(x, pnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="dystrybuanta")
lines(x, (pnorm(x,1)+pnorm(x,-1))/2, type="l", lwd=3, col="red3")

Rozkład logistyczny

plot(x, dnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="g&#281;sto&#347;&#263;")
lines(x, dlogis(x,0,0.627), type="l", lwd=3, col="red3")
plot(x, pnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="dystrybuanta")
lines(x, plogis(x,0,0.627), type="l", lwd=3, col="red3")

Rozkład lognormalny

plot(x, dnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="g&#281;sto&#347;&#263;")
lines(x, dlnorm(x+2.2,0.83,0.5), type="l", lwd=3, col="red3")
plot(x, pnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="dystrybuanta")
lines(x, plnorm(x+2.2,0.83,0.5), type="l", lwd=3, col="red3")

Rozkład alternatywa Legendre’a rzędu 2

tx <- seq(0.0001,1-0.0001,0.0001)
ty <- (ddst:::ddst.polynomial.fun[[3]](tx)+2)/2
plot(x, dnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="g&#281;sto&#347;&#263;")
lines(qnorm(tx)[-1], ty[-1]/diff(qnorm(tx))/10000, type="l", lwd=3, col="red3")
plot(x, pnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="dystrybuanta")
lines(qnorm(tx), cumsum(ty)/10000, type="l", lwd=3, col="red3")

Rozkład alternatywa Legendre’a rzędu 7

tx <- seq(0.0001,1-0.0001,0.0001)
ty <- (ddst:::ddst.polynomial.fun[[7]](tx)+2)/2
plot(x, dnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="g&#281;sto&#347;&#263;",ylim=c(0,0.53))
lines(qnorm(tx)[-1], ty[-1]/diff(qnorm(tx))/10000, type="l", lwd=3, col="red3")
plot(x, pnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="dystrybuanta")
lines(qnorm(tx), cumsum(ty)/10000, type="l", lwd=3, col="red3")

Rozkład alternatywa Legendre’a rzędu 6

tx <- seq(0.0001,1-0.0001,0.0001)
ty <- (ddst:::ddst.polynomial.fun[[6]](tx)+4)/4
plot(x, dnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="g&#281;sto&#347;&#263;",ylim=c(0,0.5))
lines(qnorm(tx)[-1], ty[-1]/diff(qnorm(tx))/10000, type="l", lwd=3, col="red3")
plot(x, pnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="dystrybuanta")
lines(qnorm(tx), cumsum(ty)/10000, type="l", lwd=3, col="red3")

Rozkład alternatywa Legendre’a rzędu 4

tx <- seq(0.0001,1-0.0001,0.0001)
ty <- (ddst:::ddst.polynomial.fun[[4]](tx)+3 )/3
plot(x, dnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="g&#281;sto&#347;&#263;")
lines(qnorm(tx)[-1], ty[-1]/diff(qnorm(tx))/10000, type="l", lwd=3, col="red3")
plot(x, pnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="dystrybuanta")
lines(qnorm(tx), cumsum(ty)/10000, type="l", lwd=3, col="red3")

Rozkład alternatywa w której część rozkładu normalnego skupiono w punkcie

dd <- dnorm(x)
dd[x < 0 &amp;amp; x > -0.32] = 0
dd[abs(x) < 0.001] = 100
plot(x, dnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="g&#281;sto&#347;&#263;", ylim=c(0,0.5))
lines(x,  dd, type="l", lwd=3, col="red3")
pp <- pnorm(x)
pp[x < 0 &amp; x > -0.32] = pnorm(-0.32)
plot(x, pnorm(x), type="l", lwd=3, las=1, ylab="",xlab="", main="dystrybuanta")
lines(x, pp, type="l", lwd=3, col="red3")

Kod w programie R rysujący wykresy mocy

library(fBasics)
library(lawstat)
library(nortest)
library(energy)
library(ddst)
#
# funkcja liczy p-wartosci dla rozwa&#380;anych test&oacute;w
# za wyj&#261;tkiem testu E, dla kt&oacute;rego liczona jest wartosc statystyki testowej
zbiorTestow <- function(x) {
  c(ad=ad.test(x)$p.value, 
      cvm = cvm.test(x)$p.value, 
      lillie=lillie.test(x)$p.value,
      ks=ks.test(scale(x),"pnorm")$p.value,
      dagoTest(x)@test$p.value,
      jbTest(x)@test$p.value,
      shapiro=shapiro.test(x)$p.value, 
      francia=sf.test(x)$p.value,
      pearson=pearson.test(x)$p.value,
      ddst=ddst.norm.test(x, compute.p=TRUE)$p.value,
      normal.e(x),
      sj.test(x)$p.value)
}
#
# scenariusze alternatyw
getData <- function(k = 1) {
  switch(k,
         '1'= runif(50),
         '2'= rcauchy(20),
         '3'= rlnorm(20),
         '4'= rgumbel(60),
         '5'= ifelse(runif(500)<0.5,rnorm(500),rnorm(500,2)),
         '6'= ifelse(runif(300)<0.5,rnorm(300),rnorm(300,0,2)),
         '7'= rnorm(25),
         '8' = rchisq(45,3),
         '9' = rweibull(20,0.8),
         '10' = rweibull(40,1.5),
         '11' = rexp(25,0.7),
         '12' = ifelse(runif(1000)<0.08,rnorm(1000),rnorm(1000,0,10)),
         '13' = round(rnorm(500)*5,0),
         '14' = runif(500)+runif(500),
         '15' = runif(1000)+runif(1000)+runif(1000),
         '16' = rlogis(350),
         '17' = rt(100,3),
         '18' = qnorm(sample(tx,80,replace=T,prob=wx)),
         '19' = qnorm(sample(tx,100,replace=T,prob=wx2)),
         '20' = qnorm(sample(tx,100,replace=T,prob=wx3)),
         '21' = qnorm(sample(tx,250,replace=T,prob=wx4)),
         '22' = qnorm(sample(tx,1000,replace=T,prob=wx5)),
         '23' = qnorm(sample(tx,100,replace=T,prob=wx6)))
}
#
# gestosci pomocnicze
tx <- seq(0.0001,1-0.0001,0.0001)
wx <- ddst:::ddst.polynomial.fun[[5]](tx)+2
wx2 <- ddst:::ddst.polynomial.fun[[3]](tx)+2
wx3 <- ddst:::ddst.polynomial.fun[[7]](tx)+2
wx4 <- ddst:::ddst.polynomial.fun[[6]](tx)+4
wx5 <- ddst:::ddst.polynomial.fun[[2]](tx)+2
wx6 <- ddst:::ddst.polynomial.fun[[4]](tx)+3
#
# poziomy istotnosci dla kolejnych testow
alfy <- c(0.05,0.01,0.01,0.05,0.05,0.05,0.01,0.01,0.01,0.05,
          0.03,0.05,0.01,0.05,0.05,0.05,0.05,0.05,0.05,0.05,
          0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05)
#
# wyznacz moce dla wszystkich alternatyw
N <- 10000
for (alte in 1:23) {
  # policz p-wartosci
  res <- replicate(N, zbiorTestow(getData(alte)) ) 
  # dla testu E wyznacz empirycznie p-wartosci
  n <- length(getData(alte))
  res0 <- replicate(N, normal.e(rnorm(n)) ) 
  res["energy",] <- sapply(res["energy",], function(r) (sum(r < res0) + 1)/(N+1))
  # rysuj wykres do pliku pdf
  pdf(paste0("moc",alte,".pdf"),9,4)
  par(mar=c(3,10,2,2))
  moce <- rowMeans(res < alfy[alte])
  th <- max(moce,na.rm=T)*0.8
  bb <- barplot(moce[-(8:9)], horiz=T,las=1, xlim=c(0,1), main="", 
                 col=ifelse(moce[-(8:9)] > th, "red3", "grey"),
                 border=ifelse(moce[-(8:9)] > th, "red3", "grey"))
  abline(v = seq(0,1,0.1), lty = 3, col = "#77777777")
  abline(v = th, lty = 2, lwd = 1, col = "red4")
  rect(moce[-(8:9)]+0.005, bb - 0.35, moce[-(8:9)]+0.08, bb + 0.4, 
            col="white", border="white")
  text(moce[-(8:9)]+0.01, bb - 0.25, substr(as.character(round(moce[-(8:9)],3)),2,10), 
            adj=c(0,0), cex=0.8)
  axis(3,th,round(th,3))
  dev.off()
}

library(fBasics) library(lawstat) library(nortest) library(energy) library(ddst) # # funkcja liczy p-wartosci dla rozważanych testów # za wyjątkiem testu E, dla którego liczona jest wartosc statystyki testowej zbiorTestow <- function(x) { c(ad=ad.test(x)$p.value, cvm = cvm.test(x)$p.value, lillie=lillie.test(x)$p.value, ks=ks.test(scale(x),"pnorm")$p.value, dagoTest(x)@test$p.value, jbTest(x)@test$p.value, shapiro=shapiro.test(x)$p.value, francia=sf.test(x)$p.value, pearson=pearson.test(x)$p.value, ddst=ddst.norm.test(x, compute.p=TRUE)$p.value, normal.e(x), sj.test(x)$p.value) } # # scenariusze alternatyw getData <- function(k = 1) { switch(k, '1'= runif(50), '2'= rcauchy(20), '3'= rlnorm(20), '4'= rgumbel(60), '5'= ifelse(runif(500)<0.5,rnorm(500),rnorm(500,2)), '6'= ifelse(runif(300)<0.5,rnorm(300),rnorm(300,0,2)), '7'= rnorm(25), '8' = rchisq(45,3), '9' = rweibull(20,0.8), '10' = rweibull(40,1.5), '11' = rexp(25,0.7), '12' = ifelse(runif(1000)<0.08,rnorm(1000),rnorm(1000,0,10)), '13' = round(rnorm(500)*5,0), '14' = runif(500)+runif(500), '15' = runif(1000)+runif(1000)+runif(1000), '16' = rlogis(350), '17' = rt(100,3), '18' = qnorm(sample(tx,80,replace=T,prob=wx)), '19' = qnorm(sample(tx,100,replace=T,prob=wx2)), '20' = qnorm(sample(tx,100,replace=T,prob=wx3)), '21' = qnorm(sample(tx,250,replace=T,prob=wx4)), '22' = qnorm(sample(tx,1000,replace=T,prob=wx5)), '23' = qnorm(sample(tx,100,replace=T,prob=wx6))) } # # gestosci pomocnicze tx <- seq(0.0001,1-0.0001,0.0001) wx <- ddst:::ddst.polynomial.fun[[5]](tx)+2 wx2 <- ddst:::ddst.polynomial.fun[[3]](tx)+2 wx3 <- ddst:::ddst.polynomial.fun[[7]](tx)+2 wx4 <- ddst:::ddst.polynomial.fun[[6]](tx)+4 wx5 <- ddst:::ddst.polynomial.fun[[2]](tx)+2 wx6 <- ddst:::ddst.polynomial.fun[[4]](tx)+3 # # poziomy istotnosci dla kolejnych testow alfy <- c(0.05,0.01,0.01,0.05,0.05,0.05,0.01,0.01,0.01,0.05, 0.03,0.05,0.01,0.05,0.05,0.05,0.05,0.05,0.05,0.05, 0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05) # # wyznacz moce dla wszystkich alternatyw N <- 10000 for (alte in 1:23) { # policz p-wartosci res <- replicate(N, zbiorTestow(getData(alte)) ) # dla testu E wyznacz empirycznie p-wartosci n <- length(getData(alte)) res0 <- replicate(N, normal.e(rnorm(n)) ) res["energy",] <- sapply(res["energy",], function(r) (sum(r < res0) + 1)/(N+1)) # rysuj wykres do pliku pdf pdf(paste0("moc",alte,".pdf"),9,4) par(mar=c(3,10,2,2)) moce <- rowMeans(res < alfy[alte]) th <- max(moce,na.rm=T)*0.8 bb <- barplot(moce[-(8:9)], horiz=T,las=1, xlim=c(0,1), main="", col=ifelse(moce[-(8:9)] > th, "red3", "grey"), border=ifelse(moce[-(8:9)] > th, "red3", "grey")) abline(v = seq(0,1,0.1), lty = 3, col = "#77777777") abline(v = th, lty = 2, lwd = 1, col = "red4") rect(moce[-(8:9)]+0.005, bb - 0.35, moce[-(8:9)]+0.08, bb + 0.4, col="white", border="white") text(moce[-(8:9)]+0.01, bb - 0.25, substr(as.character(round(moce[-(8:9)],3)),2,10), adj=c(0,0), cex=0.8) axis(3,th,round(th,3)) dev.off() }

11 thoughts on “Czternaście wybranych testów normalności”

Paweł K pisze:

9 kwietnia 2013 o 11:31

Pytanie laika – praktyka, czyli osoby, która od czasu do czasu potrzebuje policzyć jakieś statystki ze swoich wyników, a nie pasjonuje się statystyką jako dziedziną nauki, który test jest najlepszy i do czego w życiu codziennym (przykład) się je stosuje (skoro wszystkie służą do określenia zgodności rozkładu zmiennej z rozkładem normalnym)?

Na co dzień korzystam zupełnie nie świadomie z testu Shapiro-Wilka i nie zastanawiałem się czy inne testy byłyby lepsze, czy gorsze, a teraz zacząłem mieć poważne wątpliwości…

Bardzo ciekawe zestawienie (drobna uwaga na grafie Test Kołmogorowa-Smirnova nakłada się z testem SJ, poza tym napis wychodzi poza koło – może warto by wydłużyć gałąź główną od tych testów opartych o dystrybuantę empiryczną (chyba że nakładanie się tych testów było celowe), a w nazwie „test Kołmogorowa-Smirnowa” zmniejszyć czcionkę lub odległość między znakami).

Odpowiedz
1. smarterpoland pisze:
  
  9 kwietnia 2013 o 12:19
  
  Uczciwą odpowiedzią będzie, że nie ma testu najlepszego dla wszystkich alternatyw.
  Wiekszość testów nie widzi pewnego rodzaju odstepstw, test shapiro-wilka ma tą zaletę, że jest popularny[?].
  
  Inną sprawą jest, że gdy nie jesteśmy pewni czy mamy do czynienia z rozkłądem normalnym to lepiej użyć metod odpornych, które są już dosyć dostępne.
  
  Odpowiedz
Cześnik pisze:

10 kwietnia 2013 o 21:08

Proponuję zamienić wszystkie statystki na statystyki. Jest to bardzo popularny (ciekawe, dlaczego?), ale jednak błąd.

Odpowiedz
1. smarterpoland pisze:
  
  10 kwietnia 2013 o 21:16
  
  Dobry pomysł, dziękuję.
  
  Odpowiedz
piotrek pisze:

11 kwietnia 2013 o 15:54

Dobry tekst.
Drobne literówki,:
str. 19 powinno być „stablicowana”, nie „ztablicowana”
str. 19 na dole jest „verage absolute deviation from the median” – chyba „average”

Odpowiedz
1. smarterpoland pisze:
  
  11 kwietnia 2013 o 16:14
  
  dziękuję, poprawione
  
  Odpowiedz
Zbigniew Kaminski pisze:

18 maja 2013 o 13:14

Dzień dobry,
To bardzo ciekawy materiał dla osób wykonujących od czasu do czasu obliczenia satystyczne. Brakuje mi tylko informacji o stosowalności poszczególnych testów ze wzgledu na liczebność próby. W zagadnieniach technicznych liczebność jest zazwyczaj niewielka – od kilku do kilkunastu elementów w danej próbie, zwykle nie przekracza 30.
Z powazaniem

Odpowiedz
1. smarterpoland pisze:
  
  18 maja 2013 o 20:30
  
  Prezentowane studium mocy dotyczyło większych prób, ale pokazywało, że w zależności od tego jakiego typu jakiego rodzaju odstępstwo od normalności występuje, inny test jest najsilniejszy. Podobnie będzie dla mniejszych prób, może więcej stracą bardziej złożone testy, a więcej zyskają te proste [np. Shapiro Wilka].
  
  Odpowiedz
mb pisze:

8 stycznia 2014 o 21:53

Studiuję matematykę i trochę z testami miałam do czynienia. Ubolewam, że nie ma tu informacji o teście Mardii (odsyłacz nie znalazł takiego artykułu). Jednak, (zwłaszcza), informacje o teście Kołmogorowa- Smirnowa były dla mnie wartościowe. Dziękuję.

Odpowiedz
1. smarterpoland pisze:
  
  8 stycznia 2014 o 23:30
  
  Cieszę się, że się przydało.
  Testy wielowymiarowe to oddzielny temat. Aby kompleksowo do niego podejść trzeba by ten opis zwiększyć kilkukrotnie. Krzysiek napisał o teście Mardii w dokumencie http://www.biecek.pl/NaPrzelajPrzezDataMining [dziękuję za informację o niedziałających linkach, poprawiłem].
  
  Swoją drogą, biorąc pod uwagę, że Kanti Mardia to mężczyzna, to czy powinno mówić się o teście Mardii czy teście Mardiego?
  
  Odpowiedz
Marcin pisze:

8 czerwca 2020 o 15:38

Czy badając normalność za pomocą kilku testów równocześnie np. stosując test Kołmogorowa-Smirnowa i test Shapiro-Wilka należy zastosować poprawkę np. Bonferroniego?

Odpowiedz

11 thoughts on “Czternaście wybranych testów normalności”

Dodaj komentarz Anuluj pisanie odpowiedzi