Testy normalności

Testy normalności są szczególnym przypadkiem testów zgodności, czyli sytuacji gdy testujemy, czy zmienne losowe pochodzą z określonego rozkładu, ale bez zakładania ustalonych parametrów.

Na razie skupimy się na testach normalności, czyli niech \(X_1, \dots, X_n\) będzie próbą prostą z rozkładu o nieznanej, ciągłej dystrybuancie \(F_X\). Testujemy

\(H_0: \exists (m,\sigma^2)\in (\mathbb{R}, \mathbb{R}_+)\quad F_X=F_{N(0,1)},\)

przeciwko

\(H_1: \forall (m,\sigma^2)\in (\mathbb{R}, \mathbb{R}_+)\quad F_X\neq F_{N(0,1)}.\)

W tym przypadku maksymalnym niezmiennikiem jest \[(Z_1, \dots, Z_n)=\left(\frac{X_1-\bar{X}}{\sqrt{\frac{1}{n-1}\sum_{i=1}^n(X_i-\bar{X})^2}}, \dots,\frac{X_n-\bar{X}}{\sqrt{\frac{1}{n-1}\sum_{i=1}^n(X_i-\bar{X})^2}} \right)\]

względem grupy przekształceń \(G=\left\{g: \mathbb{R}^n\rightarrow \mathbb{R}^n: g(\mathbf{x})=(ax_1+b, \dots, ax_n+b), a>0, b\in \mathbb{R}\right\}\).

Zauważyć warto, ze statystyka \((Z_1, \dots, Z_n)\) przy prawdziwości hipotezy zerowej nie zależy od \((m, \sigma^2)\).

Test Lilleforsa

Test ten porównuje dystrybuantę empiryczną zestandaryzowanej próby \((Z_1, \dots, Z_n)\) z dystrybuantą \(F_{N(0,1)}\) za pomocą statystyki testowej Kołmogorowa-Smirnowa:

\[L=\sup_{t\in \mathbb{R}}\left|\hat{F}_Z(t)-F_{N(0,1)}\right|.\] Wartość statystyki najłatwiej obliczyć korzystając ze wzoru : \[L=\max_{i=1, \dots, n}\left(\max\left\{\left|\frac{i}{n}-F_{N(0,1)}(Z_{(i)})\right|, \left|\frac{i-1}{n}-F_{N(0,1)}(Z_{(i)})\right|\right\}\right).\]

P-wartość obliczamy korzystając z metody MC:

Generujemy \(mc\) prób prostych \(X_1,\dots, X_n\sim N(01,)\) (może być dowolny rozkład normalny).
Dla każdej próby z punktu \(1\) wyznaczamy wektor \((Z_1, \dots, Z_n)\) i obliczamy statystykę \(L\), otrzymując \(L_1, \dots, L_{mc}\).
\(\hat{p}=\frac{1}{mc}\sum_{i=1}^{mc}\mathbb{I}(L_i\geq L)\).

Zadanie 1

\(30\) samochodów użyto do zbadania dwóch rodzajów paliwa pod kątem spalania. Samochody zatankowały \(10\) litrów paliwa typu \(1\) i jeździły po torze do momentu aż paliwo się skończyło. Zanotowano liczbę przejechanych kilometrów. Później te same samochody zatankowały paliwo typu \(2\) i powtórzono eksperyment. Czy możemy założyć że różnice w spalaniu są normalne? (dane w pliku Baza 2.1).

set.seed(321)
library(readxl)
Baza<- read_excel("Baza 2.1.xlsx")
X<-Baza$"Typ 1"- Baza$"Typ 2"
n<-length(X)
Z<-(X-mean(X))/sqrt(sum((X-mean(X))^2)/(n-1))
Z=sort(Z)
LT<-c()
for (i in 1:n){
  LT[i]=max(abs(i/n-pnorm(Z[i],0,1)),abs((i-1)/n-pnorm(Z[i],0,1)))
}
L<-max(LT)
mc<-10000
LMC<-c()
for (k in 1:mc){
  XMC <- rnorm(n)
  ZMC<-(XMC-mean(XMC))/sqrt(sum((XMC-mean(XMC))^2)/(n-1))
  ZMC=sort(ZMC)
  LTMC<-c()
  for (i in 1:n){
    LTMC[i]=max(abs(i/n-pnorm(ZMC[i],0,1)),abs((i-1)/n-pnorm(ZMC[i],0,1)))
  }
  LMC[k]<-max(LTMC)
}
pL <- mean(LMC > L)
pL

## [1] 0.9033

Można też użyć gotowej funkcji lillie.test z biblioteki nortest.

library(nortest)
lillie.test(X)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  X
## D = 0.079161, p-value = 0.9021

Test normalności Andersona - Darlinga

Podobnie, jak wcześniej porównujemy dystrybuantę empiryczną zestandaryzowanych danych \(Z_1,\dots, Z_n\) z dystrybuantą \(F_{N(0,1)}\), ale wykorzystując do tego statystykę testową Andersona - Darlinga.

Po pierwsze wykorzystamy fakt, że zmienna losowa o ciągłej dystrybuancie, obłożona swoją własną dystrybuantą ma rozkład \(U[0,1]\).

Niech \(X\) będzie zmienną losową o ciągłej dystrybuancie \(F\) oraz niech \(Y=F(X)\), wtedy

\(F_Y(t)=P(Y\leq t)=P(X\leq F^{-1}(T))=F(F^{-1}(t))=t\), gdzie \(t\in (0,1)\).

Zatem \(Y\sim U[0,1]\).

Wobec tego, możemy równoważnie zapisać \[\sup_{t\in \mathbb{R}}\left|\hat{F}_Z(t)-F_{N(0,1)}\right|=\sup_{t\in (0,1)}\left|\hat{F}_{F_{N(0,1)}(Z)}-t\right|.\] Wykorzystując statystykę Andersona-Darlinga można pokazać że:

Zadanie 1 - cd

Wykonaj analogiczny test jak w zadaniu \(1\), ale wykorzystaj statystykę testową Andersona-Darlinga.

set.seed(321)
A<-0
for (i in 1:n){
  A=A+(2*i-1)*(log(pnorm(Z[i],0,1))+log(1-pnorm(Z[n+1-i],0,1)))
}
A=-1-1/n^2*A
mc<-10000
AMC<-c()
for (k in 1:mc){
  XMC <- rnorm(n)
  ZMC<-(XMC-mean(XMC))/sqrt(sum((XMC-mean(XMC))^2)/(n-1))
  ZMC=sort(ZMC)
  AMC[k]<-0
  for (i in 1:n){
    AMC[k]=AMC[k]+(2*i-1)*(log(pnorm(ZMC[i],0,1))+log(1-pnorm(ZMC[n+1-i],0,1)))
  }
  AMC[k]=-1-1/n^2*AMC[k]
}
pA <- mean(AMC > A)
pA

## [1] 0.9503

Zatem nie odrzucamy hipotezy zerowej. Jeżeli przeprowadzimy dodatkowo test na jednorodność wariancji i nie odrzuci on hipotezy zerowej, to będziemy mogli zastosować test t-Studenta.

Funkcja kwantyl-kwantyl

Funkcja kwantyl-kwantyl (Q-Q plot) służy do porównania rozkładu empirycznego z rozkładem teoretycznym. Na osi poziomej przedstawione są teoretyczne kwantyle, a na osi poziomej empiryczne.

Niech \(X_1, \dots, X_n\) będzie próbą prostą z rozkładu \(N(m, \sigma^2)\).

Można pokazać ze wszystkie punkty \(\left(E_{N(0,1)}(X_{(i)}), X_{(i)}\right)\) leżą na prostej \(y=\sigma^2x+m\).

Funkcja \(y=\hat{\sigma}x+\hat{m}\), gdzie \(\hat{m}=\frac{1}{n}\sum_{i=1}^nX_i\) oraz \(\hat{\sigma}=\sqrt{\frac{1}{n-1}\sum_{i=1}^n(X_i-\hat{m})}\) nazywamy funkcją kwantyl-kwantyl normalności rozkładu.

Oczywiście, empiryczne dane nie będą leżały idealnie na prostej, bo \(X_{(i)}=\sigma E_{N(0,1)}(X_{(i)})+m+\epsilon_i\).

W celu obliczenia teoretycznych kwantyli można użyć poniższego wzoru (Metoda Bloma): \[E_{N(0,1)}(X_{(i)})=F^{-1}\left(\frac{i-3/8}{n+1/4}\right).\] Empiryczne kwantyle to statystyka pozycyjna: \(X_{(1)}, \dots, X_{(n)}\).

Zadanie 2

Narysować funkcję Kwantyl-Kwantyl dla danych z zadania \(1\).

KY<-sort(X)

KX <- qnorm((1:n - 3/8) / (n + 1/4), mean = 0, sd = 1)

m=mean(X)
s=sqrt(sum((X-m)^2)/(n-1))

KK<-function(t){
  s*t+m
}

plot(KK, -3, 3)
points(KX,KY)

Zadanie 3

Narysować funkcję Kwantyl-Kwantyl dla \(100\) obserwacji wygenerowanych z rozkładu wykładniczego z parametrem \(1\).

Z <- rexp(100,1)
KY <- sort(Z)

KZ <- qnorm((1:100 - 3/8) / (100 + 1/4), mean = 0, sd = 1)

mZ = mean(Z)
sZ = sqrt(sum((Z-mZ)^2)/(100-1))

KK <- function(t){
  sZ*t+mZ
}

plot(KK, -3, 3)
points(KZ,KY)

Test Shapiro - Wilka

Jest to kolejny test normalności. Wykorzystuje on funkcję kwantyl - kwantyl do porównania empirycznych kwantyli (statystyk pozycyjnych) z ich wartościami oczekiwanymi (przy założeniu, że obserwacje pochodzą z rozkładu normalnego o nieznanych parametrach).

Idea testu:

Wartości oczekiwane leżą na prostej : \(y=\hat{\sigma}x+\hat{m}\),

Statystyki pozycyjne (po uwzględnieniu szumu) leżą “koło” prostej: \(y=\tilde{\sigma}x+\tilde{m}\).

Statystyka testowa będzie się opierać na unormowanym ilorazie \(\frac{\tilde{\sigma}}{\hat{\sigma}}\). Jeżeli \(\frac{\tilde{\sigma}}{\hat{\sigma}}\) będzie bliskie \(1\), to będzie to świadczyć o prawdziwości hipotezy zerowej. Natomiast, jeżeli \(\frac{\tilde{\sigma}}{\hat{\sigma}}\) będzie mniejsze od \(1\), to będzie to świadczyło o odrzuceniu hipotezy zerowej.

Statystyka testowa Shapiro-Wilka:

\[W=\frac{(a^tX_{()})^2}{\sum_{i=1}^n\left(X_i-\frac{1}{n}\sum_{j=1}^nX_j\right)},\] gdzie \(a=\frac{R^{-1}E}{\sqrt{E^tR^{-1}R^{-1}E}}\), \(E^t=\left(E_{N(0,1)}(X_{(1)}), \dots, E_{N(0,1)}(X_{(n)})\right)\) i \(R\) to macierz wariancji i kowariancji wektora błędów \(\epsilon\).

Warto zwrócić uwagę, że obliczenie macierzy \(R^{-1}\) jest bardzo kosztowne obliczeniowo. Dlatego też współczynniki wektora \(a\) zostały stabilicowane dla \(n\leq 50\), natomiast dla \(n>50\) statystykę \(W\) zastępuje się równoważną statystyką asymptotyczną \[\tilde{W}=\frac{\left(\sum_{i=1}^n\left((X_{(i)}-1/n\sum_{i=1}^nX_{(i)})E_{N(0,1)}(X_{(i)})\right)\right)^2}{\sum_{i=1}^n(X_{(i)}-1/n\sum_{i=1}^n X_{(i)})^2\sum_{i=1}^n\left(E_{N(0,1)}(X_{(i)})\right)}.\]

P -wartość obliczamy korzystając z metody MC:

Generujemy \(mc\) prób prostych \(X_1, \dots, X_n\sim N(0,1)\) (może być dowolny rozkład normalny).
Dla każdej próby z punktu \(1\) obliczamy statystykę \(W\) i otrzymujemy próbę prostą \(W_1, \dots, W_{mc}\).
Ponieważ Shapiro-Wilk jest testem lewostronnym, to p-wartość obliczamy według poniższego wzoru \[\hat{p}=\frac{1}{mc}\sum_{i=1}^n\mathbb{I}(W_i\leq W).\]

Zadanie 4

Przeprowadź test Shapiro-Wilko dla danych z zadania \(1\) z użyciem asymptotycznej statystyki testowej Shapiro-Wilka.

SX<-sort(X)


W<-(sum((SX-mean(SX))*KX))^2/(sum((SX-mean(SX))^2)*sum(KX^2))

mc<-10000
WMC<-c()
for (k in 1:mc){
  XMC <- rnorm(n)
  SXMC<- sort(XMC)
  WMC[k]<-(sum((SXMC-mean(SXMC))*KX))^2/(sum((SXMC-mean(SXMC))^2)*sum(KX^2))
}

pW= mean(WMC<W)
pW

## [1] 0.9467

Test Shapiro-Wilka można również przeprowadzić z użyciem gotowej funkcji shapiro.test()

shapiro.test(X)

## 
##  Shapiro-Wilk normality test
## 
## data:  X
## W = 0.98322, p-value = 0.9032

shapiro.test(Z)

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.839, p-value = 4.75e-09

Zadanie (test jednostajności)

Niech \(X_1, \dots, X_{30}\) będą próbą prostą z rozkładu \(Beta(1, b)\), gdzie \(b\in\{0.2, 0.4, \dots, 2\}\). Przetestować, czy dane pochodzą z rozkładu jednostajnego, czyli \(H_0: \exists_{a,b\in \mathbb{R}}\quad F_X=F_{U[a,b]}\), przeciwko \(H_1: \forall_{a,b, \in \mathbb{R}} \quad F_X \neq F_{U[a,b]}\).

Wyznaczyć analogiczną statystykę testową do \(\tilde{W}\) i narysować krzywą mocy empirycznej. Warto zauważyć, że dla \(b=1\) mamy rozkład jednostajny, więc moc powinna być równa założonemu poziomowi istotności (np.\(\alpha=0.05\)).

ODPOWIEDZI

Posiadane obserwacje \(X_1, \dots, X_{30}\sim Beta(1,b)\), chcemy porównać z obserwacjami generowanymi przy prawdziwości hipotezy zerowej, czyli \(Y_1, \dots, Y_{30}\sim U[0,1]\).

Zauważmy, że statystyka pozycyjna \(Y_{(k)}\sim Beta(k, n-k+1)\) (u nas \(n=30\)). Dodatkowo \(E\left(Y_{(k)}\right)=\frac{k}{n+1}\) oraz \(Var\left(Y_{(i)}\right)=\frac{k(n-k+1)}{(n+1)^2(n+2)}\).

Zatem kwadrat empirycznej korelacji dla \(\left(E_{U[0,1]}(Y_{(1)}), Y_{(1)}\right), \dots, \left(E_{U[0,1]}(Y_{(n)}), Y_{(n)}\right)\) będzie w następującej postaci

\[\tilde{W}_U=\frac{\left[E(X-E(X))E(Y-E(Y))\right]^2}{Var(X)Var(Y)}=\] \[=\frac{\left[\sum\left(Y_{(i)}-\bar{Y}_{(i)}\right)\sum\left(E(Y_{(k)})-1/n\sum E(Y_{(k)})\right)\right]^2}{\sum(Y_{(i)}-\bar{Y_{(i)}})\sum VarY_{(k)}} \] \[=\frac{\left[\sum_{i=1}^n\left(X_{(i)}-1/n\sum_{i=1}^nX_{(i)}\right) \left(\frac{i}{n+1}-1/2\right)\right]^2}{\sum_{i=1}^n \left(X_{(i)}-1/n\sum_{i=1}^nX_{(i)}\right)^2 \sum_{i=1}^n\left(\frac{i}{n+1}-1/2\right)^2}. \]

Labolatorium 9

Karolina Marek

2024-08-06