Rozdział 22 Inne testy statystyczne

22.1 Sprawdzanie zgodności z rozkładem Gaussa

Do popularnych „testów normalności” należą test Kołmogorowa-Smirnowa, test Shapiro-Wilka-Roystona, test Jarque'a-Bery.

22.1.1 Test Kołmogorowa-Smirnowa

22.1.2 Test Shapiro-Wilka-Roystona

22.1.3 Test Jarque'a-Bery

22.2 Test i przedział ufności dla współczynnika korelacji Pearsona

Współczynnik korelacji bada powiązanie pomiędzy dwoma zmiennymi ilościowymi. Może być dodatni i ujemny, przyjmuje wartości między -1 a 1. Współczynnik korelacji z próby często oznaczamy symbolem \(r\), zaś współczynnik z populacji oznaczamy grecką literą \(\rho\) („ro”, ang. „rho”).

Popularnym testem współczynnika korelacji jest test t. Test ten zakłada, że wspólny rozkład prawdopodobieństwa dwóch zmiennych to rozkład dwuwymiarowy normalny.

Hipoteza zerowa w tym teście stwierdza, że zmienne w populacji nie są skorelowane liniowo:

\[ H_0: \rho = 0, \]

a hipoteza alternatywna (dwustronna), że są:

\[ H_A: \rho \ne 0 \]

Możliwe jest również zastosowanie hipotez jednostronnych: \(H_A: \rho > 0\) lub \(H_A: \rho < 0\).

Statystyka testowa jest oparta na rozkładzie \(t\) o \(n-2\) stopniach swobody i przyjmuje postać:

\[ t = r {\sqrt {\frac {n-2}{1-r^{2}}}}, \] gdzie \(n\) to liczebność próbki, a \(r\) to współczynnik korelacji w próbie.

Wzór na przedział ufności współczynnika Pearsona jest dość złożony. Dolny kraniec uzyskujemy w następujący sposób:

\[ \frac{\text{exp}\left(ln\frac{1+r}{1-r}-\frac{2}{\sqrt{n-3}}z_{\alpha/2}\right)-1}{\text{exp}\left(ln\frac{1+r}{1-r}-\frac{2}{\sqrt{n-3}}z_{\alpha/2}\right)+1} \]

zaś górny (prawy) kraniec w następujący:

\[ \frac{\text{exp}\left(ln\frac{1+r}{1-r}+\frac{2}{\sqrt{n-3}}z_{\alpha/2}\right)-1} {\text{exp}\left(ln\frac{1+r}{1-r}+\frac{2}{\sqrt{n-3}}z_{\alpha/2}\right)+1} \]

22.4 Szablony

Sprawdzanie normalności — arkusz Google

Sprawdzanie normalności — szablon w Excelu

Współczynnik korelacji — szablon — arkusz Google

Współczynnik korelacji — szablon — szablon w Excelu

Kod w R

# Sprawdzenie normalności
# Dane
data<-c(76, 62, 55, 62, 56, 64, 44, 56, 94, 37, 72, 85, 72, 71, 60, 61, 64, 71, 69, 70, 64, 59, 71, 65, 84)
# Średnia
m<-mean(data)
# Odchylenie standardowe
s<-sd(data)
# Skośność
skew<-e1071::skewness(data, type=2)
# Kurtoza
kurt<-e1071::kurtosis(data, type=2)
# iloraz IQR do odchylenia standardowego
IQR_to_s<-IQR(data)/s
#udział obserwacji odległych o mniej niż 1 odchylenie standardowe od średniej
within1s<-mean(abs(data-m)<s)
#...o mniej niż 2 odchylenia...
within2s<-mean(abs(data-m)<2*s)
#...o mniej niż 3 odchylenia...
within3s<-mean(abs(data-m)<3*s)
# Test Shapiro-Wilka
SW_res<-shapiro.test(data)
# Test Jarque'a-Bery
JB_res<-tseries::jarque.bera.test(data)
# Wyniki
print(c('Średnia' = m, 
        'Odchylenie st.' = s,
        'Liczebność' = length(data),
        'Skośność (~=0?)' = skew,
        'Kurtoza (~=0?)' = kurt,
        'IQR/s (~=1,3?)' = IQR_to_s,
        'Reguła 68%' = within1s*100,
        'Reguła 95%' = within2s*100,
        'Reguła 100%' = within3s*100,
        'Stat. test. - test SW' = SW_res$statistic,
        'Wartość p - test SW' = SW_res$p.value,
        'Stat. test. - test JB' = JB_res$statistic,
        'Wartość p - test JB' = JB_res$p.value
))
##                         Średnia                  Odchylenie st. 
##                    65.760000000                    12.145918382 
##                      Liczebność                 Skośność (~=0?) 
##                    25.000000000                    -0.002892879 
##                  Kurtoza (~=0?)                  IQR/s (~=1,3?) 
##                     1.121315636                     0.905654036 
##                      Reguła 68%                      Reguła 95% 
##                    80.000000000                    92.000000000 
##                     Reguła 100%         Stat. test. - test SW.W 
##                   100.000000000                     0.964892812 
##             Wartość p - test SW Stat. test. - test JB.X-squared 
##                     0.520206266                     0.479578632 
##             Wartość p - test JB 
##                     0.786793609
# Wykresy
hist(data, prob=TRUE)
curve(dnorm(x, m, s), col="darkblue", lwd=2, add=TRUE)

qqnorm(data, pch=16)
qqline(data, col='darkblue')

#dane:
x <- c(4, 8, 6, 5, 3, 8, 6, 5, 5, 6, 3, 3, 3, 7, 3, 6, 4, 4, 8, 7)
y <- c(4, 6, 8, 3, 4, 6, 5, 8, 1, 4, 4, 3, 1, 6, 7, 5, 4, 4, 9, 5)
#test i przedział ufności:
cor.test(x, y, alternative = "two.sided", conf.level = 0.95)
## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = 2.56, df = 18, p-value = 0.01968
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.09607185 0.78067292
## sample estimates:
##       cor 
## 0.5166288

Kod w Pythonie

import numpy as np
from scipy import stats
from scipy.stats import skew, kurtosis

data = [76, 62, 55, 62, 56, 64, 44, 56, 94, 37, 72, 85, 72, 71, 60, 61, 64, 71, 69, 70, 64, 59, 71, 65, 84]

m = np.mean(data)

s = np.std(data, ddof=1)

skew = skew(data, bias=False)

kurt = kurtosis(data, bias=False)

IQR_to_s = stats.iqr(data) / s

within1s = np.mean(np.abs(data - m) < s)

within2s = np.mean(np.abs(data - m) < 2 * s)

within3s = np.mean(np.abs(data - m) < 3 * s)

SW_res = stats.shapiro(data)

JB_res = stats.jarque_bera(data)

results = {
    'Średnia': m, 
    'Odchylenie st.': s,
    'Liczebność': len(data),
    'Skośność (~=0?)': skew,
    'Kurtoza nadwyżkowa (~=0?)': kurt,
    'IQR/s (~=1,3?)': IQR_to_s,
    'Reguła 68%': within1s * 100,
    'Reguła 95%': within2s * 100,
    'Reguła 100%': within3s * 100,
    'Stat. test. - test SW': SW_res[0],
    'Wartość p - test SW': SW_res[1],
    'Stat. test. - test JB': JB_res[0],
    'Wartość p - test JB': JB_res[1]
}


for key, value in results.items():
    print(f"{key}: {value}")
## Średnia: 65.76
## Odchylenie st.: 12.145918381634768
## Liczebność: 25
## Skośność (~=0?): -0.002892878582795321
## Kurtoza nadwyżkowa (~=0?): 1.121315636006976
## IQR/s (~=1,3?): 0.9056540357320815
## Reguła 68%: 80.0
## Reguła 95%: 92.0
## Reguła 100%: 100.0
## Stat. test. - test SW: 0.9648927450180054
## Wartość p - test SW: 0.5202044248580933
## Stat. test. - test JB: 0.47957863171421605
## Wartość p - test JB: 0.7867936085428064
import matplotlib.pyplot as plt

# Histogram
plt.hist(data, density=True)
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = stats.norm.pdf(x, m, s)
plt.plot(x, p, 'k', linewidth=2, color="darkblue")
plt.show()

# Q-Q plot
qq_data = stats.probplot(data, dist="norm", plot=plt)
plt.show()

import scipy.stats as stats
import numpy as np

x = [4, 8, 6, 5, 3, 8, 6, 5, 5, 6, 3, 3, 3, 7, 3, 6, 4, 4, 8, 7]
y = [4, 6, 8, 3, 4, 6, 5, 8, 1, 4, 4, 3, 1, 6, 7, 5, 4, 4, 9, 5]

correlation_test = stats.pearsonr(x, y)

# Przedział ufności:
ci_tmp = stats.norm.interval(0.95, 
  loc=0.5*math.log((1+correlation_test[0])/(1-correlation_test[0])), 
  scale=1/((len(x)-3)**0.5))
confidence_interval = (np.exp(2*np.array(ci_tmp))-1)/(np.exp(2*np.array(ci_tmp))+1)

print('Test:\n\n', correlation_test, '\n\n Przedział ufności: \n\n', confidence_interval)
## Test:
## 
##  PearsonRResult(statistic=0.5166287822129882, pvalue=0.01968490591016886) 
## 
##  Przedział ufności: 
## 
##  [0.09607185 0.78067292]

22.5 Zadania

Zadanie 22.1 30 studentów próbowało odgadnąć wiek pewnej osoby na podstawie fotografii. Każdy student podawał swoje oszacowanie. Potraktujmy tych 30 studentów jak próbę losową z populacji, na temat której chcemy wnioskować. Czy należy odrzucić hipotezę zerową stwierdzającą, że rozkład liczb podawanych przez studentów jako oszacowanie wieku osoby na zdjęciu ma rozkład Gaussa? Oto liczby podane przez 30 studentów: 35, 35, 44, 36, 25, 27, 37, 35, 35, 39, 25, 34, 28, 37, 29, 39, 40, 38, 30, 32, 30, 40, 34, 27, 36, 27, 29, 34, 19, 39.