WnioskowanieStatystyczne/ Testowanie hipotez

Testowanie hipotez dotyczących jednej lub dwóch populacji

Wstęp

Schemat weryfikowania hipotez omówiony jest w wykładzie Weryfikacja hipotez statystycznych. Tu przypomnimy tylko krótko podstawowe pojęcia i decyzje, które trzeba pojąć w procedurze weryfikacji.

Hipoteza zerowa i alternatywna

Podstawą sukcesu w statystycznym testowaniu hipotez jest prawidłowe ich sformułowanie. Hipotezy muszą być rozłączne. Najczęściej jako hipotezę zerową przyjmujemy zdanie, które chcemy odrzucić, gdyż błąd takiej decyzji można kontrolować. Logika testowania jest następująca: tworzymy funkcję od zmiennych losowych, dla której przy spełnieniu przez owe zmienne hipotezy zerowej potrafimy podać prawdopodobieństwa z jakimi przyjmuje ona różne wartości. Ta funkcja nazywana jest statystyką. Następnie obliczamy wartość tej funkcji dla badanej próby. Jeśli prawdopodobieństwo osiągnięcia otrzymanej bądź jeszcze bardziej ekstremalnej wartości statystyki jest niskie to wątpimy, że nasze dane są zgodne z hipotezą zerową i jesteśmy skłonni przyjąć hipotezę alternatywną.

Wybór statystyki

Wybierając statystykę można posłużyć się następującym schematem:

• Jeżeli znamy rozkład prawdopodobieństwa, z którego pochodzą nasze dane, lub umiemy je przetransformować do znanego rozkładu, to wybierzemy klasyczny test parametryczny np. test t (ttest_rel, ttest_ind), [math]\chi^2[/math], [math]F[/math] itp.
• Jeżeli nie znamy rozkładu prawdopodobieństwa naszych danych albo nie chcemy nic o nim zakładać to mamy dwie możliwości:
  • korzystamy z klasycznego testu nieparametrycznego np.:

    test Wilcoxona dla obserwacji sparowanych

    testuje hipotezę zerową, że dwie próby [math]X[/math] i [math]Y[/math], które ze sobą porównujemy pochodzą z tej samej populacji ciągłej (przekłada się to na równość dystrybuant). Próby [math]X[/math] i [math]Y[/math] są sparowane. W pythonie mamy ten test zaimplementowany jako: scipy.stats.wilcoxon(x, y=None). Ta implementacja stosuje przybliżenie dużych prób i zalecana jest dla [math]n\gt 20[/math].

    test Manna-Whitney'a 

    testuje hipotezę zerową, że dwie próby [math]X[/math] i [math]Y[/math], które ze sobą porównujemy pochodzą z tej samej populacji ciągłej. Próby [math]X[/math] i[math]Y[/math] nie są sparowane. Implementacja w pyhtonie: scipy.stats.mannwhitneyu(x, y, use_continuity=True) stosuje przybliżenia i zalecana jest dla [math]n\gt 20[/math] w każdej z prób. Funkcja zwraca wartość p dla hipotezy jednostronnej. Jeśli testujemy hipotezę dwustronną trzeba otrzymane p pomnożyć przez 2.

  • wytwarzamy rozkład statystyki na podstawie naszych danych przez repróbkowanie (bootstrap lub permutacje). W podejściu repróbkowania tworzymy statystyczny model badanego procesu zgodny z hipotezą zerową i następnie badamy w drodze symulacji prawdopodobieństwa generowania przez ten model interesujących nas sytuacji. Największą uwagę musimy tu poświęcić na prawidłowe sformułowanie modelu, a następnie precyzyjne określenie prawdopodobieństwo jakiego zdarzenia nas naprawdę interesuje.

Poziom istotności [math]\alpha[/math] i wartość [math]p[/math]

Błąd pierwszego rodzaju 

błąd polegający na odrzuceniu hipotezy zerowej, która w rzeczywistości jest prawdziwa. Oszacowanie prawdopodobieństwa popełnienia błędu pierwszego rodzaju oznaczamy symbolem [math]\alpha[/math] i nazywamy poziomem istotności testu.

Błąd drugiego rodzaju 

polegaja na nieodrzuceniu hipotezy zerowej, która jest w rzeczywistości fałszywa. Oszacowanie prawdopodobieństwo popełnienia błędu drugiego rodzaju oznaczamy symbolem [math] \beta[/math], a jego dopełnienie do jedności nazywane jest mocą testu.

Wartość [math]p[/math] jest to wartość prawdopodobieństwa, że wobec posiadanych danych hipoteza zerowa jest prawdziwa. Jest ono obliczane jako prawdopodobieństwo zaobserwowania wartości statystyki takiej jak dla badanej próby, lub bardziej ekstremalnej, przy prawdziwej hipotezie zerowej. Najczęściej porównujemy p z wcześniej ustalonym poziomem istotności [math]\alpha[/math]. Poziom istotności [math]\alpha[/math] to wartość krytyczna prawdopodobieństwa, taka że jeżeli [math] p\lt \alpha[/math] to odrzucamy hipotezę zerową.

Formułowanie hipotez

Przykład: mutacje muszek owocowych

Załóżmy, że badamy muszki owocowe. W standardowej populacji proporcja samic do samców jest 1:1. Opracowaliśmy metodę powodującą taką mutację muszek owocowych, że potomstwo ich nie będzie miało jednakowej szansy na bycie samcem lub samiczką. W pierwszych 20 zbadanych przypadkach uzyskujemy 14 samców i 6 samiczek.

Pytanie naukowe

Czy wyniki eksperymentu potwierdzają, że nasza metoda zaburza proporcję płci?

Najpierw musimy przetłumaczyć pytanie naukowe na pytanie statystyczne.

Pytanie statystyczne

Jakie jest prawdopodobieństwo uzyskania zaobserwowanej proporcji (14:6), lub bardziej ekstremalnej w próbie 20 osobników, jeśli rzeczywista proporcja płci jest 1:1?

Z tego pytania wynikają dwie możliwe hipotezy:

• Hipoteza zerowa: Nowa metoda nie zaburza proporcji płci 1:1. Zaobserwowana próbka pochodzi z populacji, w której proporcja płci jest 1:1
• Hipoteza alternatywna: Zaobserwowana próbka pochodzi z populacji, w której proporcja płci nie jest 1:1.

Prawdopodobieństwo, które musimy oszacować

Jakie jest prawdopodobieństwo uzyskania 14 lub więcej jedynek w serii 20 prób, jeśli prawdopodobieństwo jedynki jest [math]1/2[/math]?

1. Oznaczmy 1: samiec 0: samiczka.
2. Zróbmy wektor 20-elementowy zawierający 10 zer i 10 jedynek.
3. Wylosujmy ze zwracaniem nowy wektor 20-elementowy. (Jest to nasz model uzyskiwania 20 elementowej próbki z populacji o proporcji 1:1.) Zapamiętajmy liczbę jedynek.
4. Powtórzmy poprzedni krok 1000 razy
5. Zróbmy histogram ilości jedynek.
6. Policzmy ile razy zdarzyło sie 14 lub więcej jedynek (to odpowiada 14 lub więcej samców) i dodajmy do tego ilość przypadków gdy mieliśmy 6 lub mniej jedynek (to odpowiada 14 lub więcej samiczek). Wynik podzielmy przez ilość losowań (1000).

Powyższa procedura opisuje test dwustronny. Testu dwustronnego musimy użyć jeśli nie mamy istotnych powodów, żeby wierzyć, że nowa metoda działa jedynie na zwiększenie szansy pojawienia się samca.

Jak powyższy problem rozwiązać w sposób parametryczny? Próby podlegają rozkładowi dwumianowemu, można więc obliczyć szukane prawdopodobieństwo korzystając z jego dystrybuanty.

Testowanie hipotez na temat średniej

Firma rozwożąca przesyłki po mieście deklaruje średni czas dostarczenia 28 minut. Przeprowadźmy test tej hipotezy na poziomie istotności 5%.

[math]H_0: \mu = 28[/math] [math]H_1: \mu \ne 28[/math]

Wybieramy losową próbę 100 przesyłek, mierzymy czas dostarczenia, liczymy średnią z próby [math]\bar x = 31,5[/math] minut i odchylenie standardowe próby [math]s = 5[/math] minut.

Test ten można przeprowadzić z użyciem przedziałów ufności:

• Konstruujemy przedziały ufności 95% dla średniej [math]\mu[/math].Formalnie, ponieważ wyliczyliśmy odchylenie standardowe z próby powinniśmy zastosować wartości krytyczne z rozkładu t. Rozmiar próby wynosi 100 więc rozkład t ma 99 stopni swobody. Wartość krytyczna w tym rozkladzie to [math]t_{97.5} = 1.984 [/math]. Mamy zatem:

[math]\bar x \pm 1,984 \frac{s}{\sqrt{n-1}} = 31,5 \pm 1,984*\frac{5}{\sqrt{99}} = 31,5 \pm 0,997 = [30,50 \quad 32,49][/math]

• Wnioskowanie: Możemy więc być na 95% pewni, że nieznana średnia leży pomiędzy 30,50 a 32,49 a więc na 95% nie leży poza tym przedziałem.

Skoro [math]H_0[/math] podaje [math]\mu = 28[/math] (poza przedziałem), możemy odrzucić tą hipotezę. Jeśli [math]H_0[/math] jest prawdziwe, istnieje prawdopodobieństwo 0,05, że skonstruowany przedział nie będzie zawierał [math]\mu[/math]. Istnieje zatem prawdopodobieństwo 0,05 popełnienia błądu I-go rodzaju. Mówimy, że przeprowadziliśmy test na poziomie istotności 0,05.

Test można też przeprowadzić następująco:

• Jako statystykę wybieramy [math]t = \frac{\bar x - \mu}{s/\sqrt{n-1}} [/math]
• Obliczmy jej wartość dla danych z próby: [math]t = \frac{\bar x - \mu}{s/\sqrt{n-1}} = (31,5 -28)/(5/\sqrt{99}) = 6.96 [/math]
• Porównujemy otrzymaną wartość z wartością krytyczna dla przyjętego poziomu istotności [math]\alpha[/math]. Konstuowany przez nas test ma być testem dwustronnym więc musimy wziąć do porównania wartość krytyczna dla [math]\alpha/2[/math]: [math] t_{\alpha/2} = t_{0,025}= -1,984[/math]
• Wnioskowanie: Dla naszej próby otrzymaliśmy wartość statystyki równą 6.96. Dla przyjętego poziomu istotności wartość statystyki wynosi 1,984. Prawdopodobieństwo zaobserwowania statystyki o wartości 6.96 lub bardziej ekstremalnej (tu: większej) jest mniejsze niż założony poziom istotności więc odrzucamy [math]H_0[/math] na poziomie [math]\alpha = 0,05[/math].

Testowanie hipotez na temat wariancji

Do testowania hipotez na temat wariancji używamy statystyki chi-kwadrat o [math]n-1[/math] stopniach swobody:

[math]\chi^2 = \frac{(n-1)s^2}{\sigma_0^2} [/math]

gdzie [math]\sigma_0^2[/math] jest wartością wariancji podaną w [math]H_0[/math].

Przykład

Do produkcji baterii używane są metalowe płytki o średniej średnicy 5mm. Jeśli wariancja średnicy płytki jest nie większa niż 1mm[math]^2[/math], produkcja jest kontynuowana. Jeśli wariancja przekracza 1mm[math]^2[/math] proces produkcji trzeba przerwać. Kontroler przeprowadza test na poziomie istotności [math]\alpha = 5[/math]%:

[math]H_0:[/math] [math]\sigma^2 \le 1[/math] i [math]H_1:[/math] [math]\sigma^2 \gt 1[/math].

Wybiera losową próbę 31 płytek i znajduje wariancję próby [math]s^2 = 1,62[/math]. Czy daje to podstawy do przerwania produkcji ?

[math]\chi^2 = \frac{(31-1) \cdot 1,62}{1} = 48,6[/math]. Znajdujemy poziom [math]p[/math] dla tej wartości [math]\chi^2[/math] z 30 stopniami swobody.

import scipy.stats as st
import pylab as py
p = 1-st.chi2.cdf(48.6,30)
print p

#ilustracja
x = py.arange(0,60,1)  
py.plot(x, st.chi2.pdf(x,30)) #rysujemy funkcję gęstości prawdopodobieństwa chi2 o 30 st. swobody
x2 = py.arange(48.8,60,1)
py.fill_between(x2,st.chi2.pdf(x2,30)) #cieniujemy obszar pod funkcją geęstości prawdopodobieństwa odpowiadający obliczonemu p

py.show()

Otrzymaliśmy [math] p = 0,017 [/math]. Wniosek: Odrzucamy [math]H_0[/math]

Błąd drugiego rodzaju. Moc testu.

Błąd II-go rodzaju popełniamy wtedy gdy przyjmujemy [math]H_0[/math] a prawdziwe jest [math]H_1[/math].

Przykład

Załóżmy następujący test:

[math]H_0: \mu = 60[/math]

[math]H_1: \mu = 65[/math]

Niech rozmiar próby wynosi [math]n = 100[/math] a odchylenie standardowe w populacji [math]\sigma = 20[/math].

Powinniśmy tu zastosować test jednostronny (mamy tylko dwie możliwości: [math]\mu = 60[/math] lub [math]65[/math]).

Znajdźmy punkt krytyczny [math]C[/math] dla poziomu istotności [math]\alpha = 0,05[/math]:

[math]C = \mu_0 + 1,645\frac{\sigma}{\sqrt{n}} = 60 + 1,645(20/10) = 63,29[/math]

Błąd pierwszego rodzaju popełnimy gdy [math]\bar x \gt C[/math] i prawdziwe będzie [math]H_0[/math]. Prawdopodobieństwo błędu pierwszego rodzaju ustaliliśmy z góry na poziomie [math]\alpha = 0,05[/math].

[math]\alpha = P(\bar x \gt C| \mu = \mu_0)[/math]

Błąd drugiego rodzaju popełnimy gdy [math]\bar x \lt C[/math] a prawdziwe będzie [math]H_1[/math]. Prawdopodobieństwo popełnienia tego błędu wynosi: [math]\beta = P(\bar x \lt C| \mu = \mu_1) = P\left( \frac{\bar x - \mu_1}{\sigma/\sqrt{n}} \lt \frac{C - \mu_1}{\sigma/\sqrt{n}}\right)= [/math]

[math]= P\left( Z \lt \frac{63,29 - 65}{20/10}\right) = P(Z \lt -0,855) = 0,1963[/math]

Moc testu jest dopełnieniem błędu II rodzaju. A moc testu czyli prawdopodobieństwo odrzucenia hipotezy zerowej podczas gdy jest ona fałszywa wynosi: [math]1 - \beta = 0,8037[/math].

Porównanie dwóch populacji

Dla przypomnienia: Jeśli mamy dwie próbki danych:

[math]x_1[/math] o liczebności [math]n_1[/math] i estymowanej wariancji [math]s_1^2[/math] i

[math]x_2[/math] o liczebności [math]n_2[/math] i estymowanej wariancji [math]s_2^2[/math]

gdzie [math]s_{i}^2= \frac{1}{(n_i-1)}\sum\limits_{j=1}^{n_i} (x_{j}-\overline{x}_{i})^{2}[/math]

pochodzących z rozkładu normalnego o takiej samej wariancji [math]\sigma[/math] to:

• wspólna wariancja może być estymowana jako: [math]s^2=\frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2} [/math]
• wariancja różnicy średnich może być estymowana jako: [math]s_{\Delta}^2=s_{\bar{x}_1}^2+s_{\bar{x}_2}^2=\frac{1}{n_1}s_1^2+\frac{1}{n_2}s_2^2 = \frac{n_1+n_2}{n_1n_2}s^2[/math]

Jeśli postawimy hipotezę zerową: [math]H_0:\; \mu_1 = \mu_2[/math] to

[math] t= \frac{\bar{x}_1-\bar{x}_2}{s_{\Delta}} [/math]

pochodzi z rozkładu [math]t[/math] o [math]n_1+n_2-2[/math] stopniach swobody.

Przykład: promocja

Producent odtwarzaczy CD chce sprawdzić czy małe obniżenie ceny produktu wpłynie korzystnie na sprzedaż. Losowa próba 15 tygodni sprzedaży przed obniżką dała średni dochód 6598 zł i standardowe odchylenie 844 zł. Losowa próba 12 tygodni sprzedaży w trakcie promocji dała średnią 6870 i odchylenie standardowe 669 zł. Czy dane te wykazują poprawę sprzedaży w trakcie promocji (przyjąć poziom istotności 5%) ?

[math]\bar{x}_1 = 6598[/math]

[math]\bar{x}_2 = 6870[/math]

Treść zadania odpowiada testowi jednostronnemu z poniższymi hipotezami:

[math]H_0: \mu_1 = \mu_2 [/math]

[math]H_1: \mu_1 \lt \mu_2 [/math]

Różnica średnich wynosi: [math]\bar{x}_1-\bar{x}_2 = 6598 - 6870 = -272 [/math]

Estymowana wariancja różnicy średnich: [math]s_{\Delta}^2 = \frac{1}{n_1}s_1^2+\frac{1}{n_2}s_2^2 = \frac{844^2}{15} + \frac{669^2}{12} = 84785,82[/math]

zatem:

[math]s_{\Delta} = 291,18[/math]

[math] t = -272/ 291,18 = - 0,934[/math]

Ilość stopni swobody: [math]n_1+n_2-2 = 25[/math]

[math]P(t_{25} \le -0,934) = 0.18 [/math]

Wniosek: Nie możemy odrzucić [math]H_0[/math]. Nie mamy podstaw by uznać że mała obniżka cen poprawiła sprzedaż.

Do testowania równości wariancji w dwóch populacjach stosuje się test [math]F[/math]: [math]F_{(n1 -1, n2 - 1)} = \frac{s_1^2}{s_2^2}[/math]

W przykładzie powyżej założyliśmy równość wariancji. Korzystając z testu [math]F[/math] sprawdzić czy założenie było uzasadnione.

Badanie założenia o normalności rozkładu

We wszystkich wspomnianych powyżej klasycznych testach statystycznych [math]t[/math], [math]z[/math], [math]F[/math], [math]\chi^2[/math] istotnym założeniem jest to, że dane wejściowe w próbie mają rozkład normalny. W powyższych zadaniach po prostu to zakładaliśmy, ale w praktyce, kiedy dostajemy próbę do analizy, musimy sami sprawdzić, czy możemy uznać ją za pochodzącą z rozkładu normalnego. Do weryfikacji takiej hipotezy służą narzędzia graficzne:

Histfit: histogram z naniesionym fitem rozkładu normalnego

def histfit(x,N_bins):  
''' 
x - dane
N_bins -ilość binów w histogramie

Funkcja rysuje histogram i na jego tle dorysowuje wykres 
funkcji gęstości prawdopodobieństwa rozkładu normalnego 
o średniej i wariancji estymowanych z x.

Funkcja wymaga zaimportowania modułów pylab as py i scipy.stats as st''' 
 
	n, bins, patches = py.hist(x, N_bins, normed=True, facecolor='green', alpha=0.75)
	# Rysujemy histogram i w jawny sposób odbieramy zwracane przez p.hist obiekty
	#   - normujemy histogram do jedności
	#   - ustalamy kolor prostokątów na zielony
	#   - ustawiamy przezroczystość prostokątów na 0.75


	bincenters = 0.5*(bins[1:]+bins[:-1])
	# wytwarzamy tablicę z centrami binów korzystając z granic binów
	# zwróconych przez py.hist w macierzy bins

 
	y = st.norm.pdf( bincenters, loc = np.mean(x), scale = np.std(x))
	# obliczamy momenty rozkładu x: średnią i wariancję (tak naprawdę to jej pierwiastek czyli standardowe odchylenie)
	# obliczamy wartości w normalnym rozkładzie gęstości prawdopodobieństwa
	# o średniej np.mean(x) i standardowym odchyleniu np.std(x) dla wartości bincenters
 
	l = py.plot(bincenters, y, 'r--', linewidth=1)
	# do histogramu dorysowujemy linię

Normplot

Jest to wykres wartości w próbie, wzg. prawdopodobieństwa uzyskania takiej wartości w rozkładzie normalnym. Wykres ten jest szczególnym przypadkiem kwantylowego wykresu prawdopodobieństwa [1]. Konstruuje się go następująco:

• Na osi pionowej odkładamy uporządkowane zaobserwowane wartości [math]x_1 \le x_2 \le \dots \le x_n [/math]. Estymują one położenie kwantyli w populacji.
• Na osi poziomej odkładamy położenia kwantyli w rozkładzie normalnym. Obliczamy je: [math]F^{-1}\left(\frac{i-r_{pop}}{n+n_{pop}} \right)[/math] gdzie [math]F^{-1}[/math] jest funkcją odwrotną do dystrybuanty rozkładu teoretycznego (w tym szczególnym przypadku normalnego) [math]i[/math] jest numerem obserwacji a [math]r_{pop}, n_{pop}[/math] są czynnikami ([math]\le 0.5 [/math]) zapewniającymi, że argument funkcji odwrotnej do dystrybuanty nie przyjmuje wartości 0 ani 1.

Jeśli dane x podlegają rozkładowi normalnemu to ich położenia kwantyli powinny być zgodne z położeniami kwantyli w rozkładzie normalnym, zatem wykres powinien przedstawiać linię prostą. Odstępstwa od prostej świadczą o odstępstwie od rozkładu normalnego. W pythonie możemy ten typ wykresu zaimplementować np. tak (współczynniki zaczerpnięte z [2]):

# -*- coding: utf-8 -*-
import scipy.stats as st
import pylab as py
import numpy as np

def normplot(x):
	'''normplot: x dane do testowania'''
	x_ord = sorted(x)
	N = len(x)
	y = np.zeros(N)
	y[0]=st.norm.ppf(1- 0.5**(1.0/N) )
	y[N-1] = st.norm.ppf(0.5**(1.0/N) )
	for i in range(1,N-1):
		arg = float(i-0.3175)/(N+0.365)
		y[i] = st.norm.ppf(arg)
	py.plot(y,x_ord,'.')

Testy nieparametryczne:

Test Shapiro-Wilka

W,p = st.shapiro(x)

Funkcja ta zwraca wartość statystyki W i prawdopodobieństwo p zaobserwowania takiej lub bardziej ekstremalnej wartości statystyki W dla danych podlegających rozkładowi normalnemu. Jeśli p ma wysoką wartość (większą niż przyjęta [math]\alpha[/math]) to nie możemy odrzucić hipotezy, że dane podlegają rozkładowi normalnemu.

Test Kolmogorova-Smirnova

Jest on oparty na badaniu maksymalnej różnicy pomiędzy dystrybuantą empiryczną (z próby) a teoretyczną dystrybuantą rozkładu. Dla testowania normalności próby powinniśmy wywołać

D , p =st.kstest(x, 'norm')

Funkcja zwraca wartość statystyki D i prawdopodobieństwo zaobserwowania takiej bądź bardziej ekstremalnej wartości statystyki jeśli testowane dane pochodzą z tego samego rozkładu.

Przykład

Proszę wygenerować 1000 liczb z rozkładu normalnego. Liczby te obejrzyjmy na wykresach histfit oraz normplot i zbadamy ich normalność testem Shapiro-Wilka i Kołmogorova-Smirnova.

# -*- coding: utf-8 -*-
import scipy.stats as st
import pylab as py
import numpy as np

def normplot(x):
	'''normplot: x dane do testowania'''
	x_ord = sorted(x)
	N = len(x)
	y = np.zeros(N)
	y[0]=st.norm.ppf(1- 0.5**(1.0/N) )
	y[N-1] = st.norm.ppf(0.5**(1.0/N) )
	for i in range(1,N-1):
		arg = float(i-0.3175)/(N+0.365)
		y[i] = st.norm.ppf(arg)
	py.plot(y,x_ord,'.')
	
def histfit(x,N_bins):  
	''' funkcja rysuje histogram i na jego tle dorysowuje wykres 
	funkcji gęstości prawdopodobieństwa rozkładu normalnego 
	o średniej i wariancji estymowanych z x 
	Funkcja wymaga zaimportowania modułów pylab as py i scipy.stats as st'''  
	n, bins, patches = py.hist(x, N_bins, normed=True, facecolor='green', alpha=0.75)
	bincenters = 0.5*(bins[1:]+bins[:-1]) 
	y = st.norm.pdf( bincenters, loc = np.mean(x), scale = np.std(x)) 
	l = py.plot(bincenters, y, 'r--', linewidth=1)
def porownania(x):
	py.subplot(2,2,1)
	normplot(x);
	py.subplot(2,2,2)
	histfit(x,15)
	W,p_sw = st.shapiro(x);
	D,p_ks = st.kstest(x,'norm',args=(np.mean(x),np.std(x) ))
	title = 'SW: %(sw).2f KS: %(ks).2f'%{'sw':p_sw, 'ks':p_ks}
	py.title(title)
	# wybieramy pierwsze dziesięć punktów  
	y=x[0:10]
	py.subplot(2,2,3)
	normplot(y);
	py.subplot(2,2,4)
	histfit(y,15)
	W,p_sw = st.shapiro(y);
	D,p_ks = st.kstest(y,'norm');
	title = 'SW: %(sw).2f KS: %(ks).2f'%{'sw':p_sw, 'ks':p_ks}
	py.title(title)

x = st.norm.rvs(size=1000,loc=0,scale=10);
py.figure(1)
porownania(x)
# A teraz zbadajmy dane z rozkładów innych niż normalny:

x = st.t.rvs(2,size=1000,loc=0,scale=1);
py.figure(2)
porownania(x)

x = st.expon.rvs(2,size=1000,loc=0,scale=1);
py.figure(3)
porownania(x)


py.show()

Proszę zapuścić skrypt kilka razy i zwrócić uwagę na to, jak trudno jest ocenić normalność danych przy małych próbach.

Przykład: transformacja Boxa-Coxa

Często normalność danych można poprawić przez zastosowanie odpowiedniej transformacji. Ogólną rodzinę transformacji, które często prowadzą do normalizacji danych można zapisać tak (trasformacja Box-Cox'a [3]):

[math] y^{(\lambda)}=\left\{\begin{matrix} \frac{y^\lambda-1} {\lambda} & \mbox{dla }\lambda \ne 0 \\ \ln(y) & \mbox{dla }\lambda = 0\end{matrix}\right. [/math]

W module scipy.stats mamy tę transformację zaimplementowaną jako boxcox().

Zobaczmy jej działanie na następującym przykładzie:

Proszę pobrać i zapisać w pliku tekstowym następujące dane. Zawierają one 8 kolumn charakterystyk samochodów:

• zużycie paliwa
• cylindry
• pojemność skokowa
• moc w koniach mechanicznych
• masa
• przyspieszenie
• rocznik
• pochodzenie

Proszę narysować histfit i normplot oraz policzyć testy Shapiro-Wilka i Kolmogorowa-Smirnowa dla masy pojazdów, a następnie zastosować do niej transformację Boxa-Coxa i zbadać normalność przetransformowanych danych.

# -*- coding: utf-8 -*-
import scipy.stats as st
import pylab as py
import numpy as np

def normplot(x):
	'''normplot: x dane do testowania'''
	x_ord = sorted(list(x))
	N = len(x)
	y = np.zeros(N)
	y[0]=st.norm.ppf(1- 0.5**(1.0/N) )
	y[N-1] = st.norm.ppf(0.5**(1.0/N) )
	for i in range(1,N-1):
		arg = float(i-0.3175)/(N+0.365)
		y[i] = st.norm.ppf(arg)
	py.plot(y,x_ord,'.')
	
def histfit(x,N_bins):  
	''' funkcja rysuje histogram i na jego tle dorysowuje wykres 
	funkcji gęstości prawdopodobieństwa rozkładu normalnego 
	o średniej i wariancji estymowanych z x 
	Funkcja wymaga zaimportowania modułów pylab as py i scipy.stats as st'''  
	n, bins, patches = py.hist(x, N_bins, normed=True, facecolor='green', alpha=0.75)
	bincenters = 0.5*(bins[1:]+bins[:-1]) 
	y = st.norm.pdf( bincenters, loc = np.mean(x), scale = np.std(x)) 
	l = py.plot(bincenters, y, 'r--', linewidth=1)
def porownania(x):
	py.subplot(2,1,1)
	normplot(x);
	py.subplot(2,1,2)
	histfit(x,15)
	W,p_sw = st.shapiro(x);
	D,p_ks = st.kstest(x,'norm',args=(np.mean(x),np.std(x) ))
	title = 'SW: %(sw).2f KS: %(ks).2f'%{'sw':p_sw, 'ks':p_ks}
	py.title(title)
	

dane = np.loadtxt('Samochody.dat')
# Badamy przyspieszenia
w = dane[:,5]
py.figure(1)
porownania(w)
# A teraz stosujemy transformację Box-Coxa
wt,lam = st.boxcox(w)
py.figure(2)
porownania(wt)

# badamy masy
w = dane[:,4]
py.figure(3)
porownania(w)
# A teraz stosujemy transformację Box-Coxa
wt,lam = st.boxcox(w)
py.figure(4)
porownania(wt)

py.show()

W pierwszym przypadku - asymetrię rozkładu przyspieszeń udało się transformacją B-C poprawić, ale w drugim przypadku, masa, asymetrię dało się skorygować (widać to na histfit) ale nie da się poprawić ciężkich ogonów - widać to zarówno na normplocie jak i na wynikach testów. Ogólnie: zanim zaczniemy analizować dane dobrze jest je pooglądać na różnych wykresach i chwilę pomyśleć.

Przykład (zastosowanie różnych testów do tych samych danych): karma

Badamy dwie nowe karmy A i B. Mamy dwie grupy po 12 zwierząt. Uzyskane przyrosty masy są następujące:

A: 31 34 29 26 32 35 38 34 31 29 32 31

B: 26 24 28 29 30 29 31 29 32 26 28 32

Pytanie: Czy któraś z karm daje istotnie większe przyrosty masy?

Poniżej rozwiążemy to zadanie stopniowo różnymi metodami. Kolejne kawałki kodu można dopisywać do tego samego pliku.

ROZWIĄZANIE: Przyjmujemy poziom istotności, na którym przeprowadzamy testy [math]\alpha = 0.05[/math].

Badamy rozkłady danych:

import scipy.stats as st
import pylab as py
import numpy as np

A=[ 31, 34, 29, 26, 32, 35, 38, 34, 31, 29, 32, 31];
B=[ 26, 24, 28, 29, 30, 29, 31, 29, 32, 26, 28, 32];
W, p_A = st.shapiro(A)
print 'Dla grupy A:', p_A
W, p_B = st.shapiro(B)
print 'Dla grupy B:', p_B

Dla obu grup test Shapiro-Wilka nie daje nam podstaw do odrzucenia założenia o normalności rozkładów.

Test parametryczny

Nie odrzuciliśmy hipotezy o normalnym rozkładzie danych zatem możemy zastosować test t dla różnicy średnich.

Formułujemy hipotezy:

• [math]H_0[/math]: średni przyrost masy w grupie A [math]=[/math] średni przyrost masy w grupie B
• [math]H_1[/math]: średni przyrost masy w grupie A [math]\ne[/math] średni przyrost masy w grupie B

Przeprowadzamy test:

t, p = st.ttest_ind(A,B)

Otrzymujemy p = 0.01.
[math]p \lt \alpha [/math], zatem na przyjętym poziomie istotności odrzucamy hipotezę zerową i stwierdzamy, że grupa A ma inną średnią niż grupa B.

Test nieparametryczny

Nie zakładajac postaci rozkładu danych mozemy zastosować test ze statystykami opartymi na rangach. Formułujemy hipotezy:

• [math]H_0[/math]: mediana przyrostu masy w grupie A [math]=[/math] mediana przyrostu masy w grupie B
• [math]H_1[/math]: mediana przyrostu masy w grupie A[math]\ne[/math] mediana przyrostu masy w grupie B

Przeprowadzamy test:

U, p = st.mannwhitneyu(A, B)
p_dwustronne = 2*p

W pythonie zaimplementowana jest wersja jednostronna tego testu. Aby otrzymać prawdopodobieństwo p dla testu dwustronnego musimy pomnożyć je przez 2.

Testy bootstrapowe

Teraz to samo sprawdzimy za pomocą testu repróbkowanego. Przyda nam się tu funkcja do pobierania losowej próbki z powtórzeniami z danych:

def randsample(x, N):
	'''zwraca wektor o dłougości N z losowo wybranymi elementami wektora x. Losowanie odbywa się z powtórzeniami''' 
	n=len(x)
	ind = np.random.randint(n, size = N)
	y = x[ind]
	return y

W testech repróbkowanych statystykę możemy wybrać dość dowolnie, ale jak pokażemy poniżej nie każda jest równie dobra. Zgodnie z hipotezą zerową próbka A i B pochodza z tej samej populacji. Nasza najlepsza wiedza o owej populacji to połączone próbki A i B:

POP=np.concatenate((A, B))

N=len(POP)
NA=len(A)
NB=len(B)
# Zasymulujemy N_rep razy wyciagniecie z POP prob o rozmiarach NA i NB i
# zobaczymy jak czesto zdarzają się wartości statystyki  sie roznica srednich taka jak w oryginalnym
# pomiarze lub jeszcze wieksza. 
N_rep=10000
# oryginalna roznica srednich i median:
mi_0 = np.abs(np.mean(A) - np.mean(B))
T_0  = np.abs(np.mean(A) - np.mean(B))/np.std(POP)
me_0 = np.abs(np.median(A) - np.median(B))

mi = np.zeros(N_rep)
T  = np.zeros(N_rep)
me = np.zeros(N_rep)

for i in range(N_rep):
	AA = randsample(POP,NA)
	BB = randsample(POP,NB)
	R_POP = np.concatenate((AA,BB))
	mi[i] = np.abs(np.mean(AA)-np.mean(BB)) # abs bo test dwustronny
	T[i]  = np.abs(np.mean(AA)-np.mean(BB))/np.std(R_POP)
	me[i] = np.abs(np.median(AA)-np.median(BB))

p_mi = np.sum(mi>=mi_0)/float(N_rep)
p_T  = np.sum( T>=T_0 )/float(N_rep)
p_me = np.sum(me>=me_0)/float(N_rep)

print 'testy repróbkowane: '
print 'rożnica średnich: ', p_mi
print 'pseudo T: ', p_T
print 'różnica median: ', p_me

Widzimy, że dla testu ze statystyką różnicy średnich i pseudo T dostajemy podobne wyniki, z tym, że pseudo T jest nieco silniejszy. Test ze statystyką różnicy średnich jest na tyle słaby, że nie pozwala na odrzucenie hipotezy zerowej.

Wnioskowanie w oparciu o przedziały ufności

Wnioskowanie o równości średnich dwóch grup można też przerowadzić w oparciu o przedziały ufności. W naszym przykładzie przedziały można skonstruować zarówno parametrycznie jak i nieparametrycznie.

Wersja parametryczna

Konstruujemy 95% przedział ufności wokół oryginalnej różnicy średnich. Dla dwóch grup wariancję różnicy średnich znajdujemy sumując wariancje śrenich estymowane dla każdej z grup:

[math] \mathrm{var}_{\Delta} = \frac{1}{N_1}\mathrm{var}(x_1) + \frac{1}{N_2} \mathrm{var}(x_2) [/math]

ilość stopni swobody: [math] df = N_1+N_2-2 [/math] co prowadzi do estymatora

[math] s^2_\Delta = \frac{N_1s_1^2 + N_2s_2^2}{N_1+N_2-2} \cdot \frac{N_1+N_2}{N_1N_2} [/math]

roznica_oryginalna = np.mean(A) - np.mean(B)
f = NA+NB-2;
v_A = np.var(A)
v_B = np.var(B)
sig = np.sqrt( ((NA*v_A +NB*v_B))/f * (NA + NB)/float(NA*NB) )
t_2_5 = st.t.ppf(0.025,f);
t_97_5 = st.t.ppf(0.975,f);
print 'przedział ufności dla różnicy średnich przy założeniu normalności %(d).2f %(g).2f'%{'d':sig*t_2_5+ roznica_oryginalna,'g':sig*t_97_5+roznica_oryginalna}

Otrzymujemy w wyniku przedział ufności dla różnicy średnich przy założeniu normalności 0.77 5.56. Oznacza to, że w 95% analogicznych badań powinniśmy otrzymać różnicę średnich zawartą w tym przedziale. Przedził ten nie zawiera wartości 0, zatem na przyjętym poziomie istotności średnie grupy A i B są różne.

Wersja nieparametryczna

# POP zawiera świat zgodny z H0
roznica_oryginalna = np.mean(A) - np.mean(B)
alfa = 0.05
N_rep = 10000
r = np.zeros(N_rep)
for i in xrange(N_rep):
	gA = randsample(POP, NA)
	gB = randsample(POP, NB)
	r[i] = np.mean(gA) - np.mean(gB)
ci_d = st.scoreatpercentile(r, per = alfa/2*100)
ci_g = st.scoreatpercentile(r, per = (1-alfa/2)*100)
# print 'przedział ufności: %(d).2f %(g).2f'%{'m':np.mean(r),'d':ci_d+roznica_oryginalna, 'g':ci_g+roznica_oryginalna}
print 'przedział ufności: %(d).2f %(g).2f'%{'m':np.mean(r),'d':ci_d, 'g':ci_g}
print 'oryginalna różnica średnich: %(ro).2f'%{'ro':roznica_oryginalna}

Wynik: przedział ufności: -2.50 2.50 oryginalna różnica średnich: 3.17. Przedział ufności 95% na różnicę skonstruowany zgodnie z [math]H_0[/math] nie zawiera oryginalnej różnicy średnich, zatem różnica 3.17 w świecie zgodnym z [math]H_0[/math] zdarza się nie częściej niż w 5% przypadków. Wniosek: średnie grup A i B są różne na przyjętym poziomie ufności.

Zadanie: Przeżywalność myszy

Mamy 7 myszy, którym podano środek, który miał poprawić ich przeżywalność po operacji oraz 9 myszy kontrolnych, którym owego środka nie podano. Myszy traktowane specjalnie przeżyły
94 38 23 197 99 16 141 dni
a myszy traktowane standardowo:
52 10 40 104 51 27 146 30 46 dni
Średnia różnica wynosi 30,63 dni dłużej dla myszy traktowanych po nowemu. Pytanie, na które chcielibyśmy znać odpowiedź to: Czy nowy środek faktycznie poprawia przeżywalność. Zadanie proszę rozwiązać wszystkimi możliwymi sposobami, analogicznie do powyższego przykładu.

Odp:

przedział ufności: [−51,64 52,89]

oryginalna różnica średnich: 30,63

przedział ufności dla różnicy średnich przy założeniu normalności [−27,99 89,26]

*

Przykład: Linie lotnicze

Linie lotnicze, projektując nowy samolot chcą sprawdzić czy średnia waga bagażu ręcznego zabieranego przez pasażerów nie zmieniła się od czasu poprzednich badań i wynosi wciąż 12 kg. Analiza ma być przeprowadzona na poziomie istotności [math]\alpha = 0,05[/math]. Analityk pobrał próbę bagażu ręcznego 144 pasażerów i obliczył wartość średnią z próby [math]\bar x = 14,6[/math] kg i odchylenie standardowe z próby [math]s = 7,8[/math]. Przeprowadź test hipotezy, że [math]\mu = 12[/math].

Zadanie: Agencja nieruchomości

Agencja nieruchomości w Japonii podała, że ceny gruntu w centrum Tokio wzrosły o 49% w ciągu ostatniego roku. Inwestor chcąc przetestować te dane, znajduje próbę 18 nieruchomości w centrum Tokio, dla których zna cenę obecna i sprzed roku. Zakłada, że ceny podlegają rozkładowi normalnemu. Dla każdej nieruchomości oblicza procentowy wzrost wartości a następnie znajduje średnią i odchylenie standardowe z próby. Estymatory dla próby wynoszą [math]\bar x = 38[/math]% i [math]s = 14[/math]%. Przeprowadź test na poziomie istotności [math]\alpha = 0,01[/math].

*

Odp. Odrzucamy [math]H_0:[/math] [math]\mu_0 = 49[/math], na poziomie istotności 0,01.

Zadanie: Zabiegi bio-inżynieryjne

Załóżmy, że krowy są bardziej wartościowe od byków. Bio-inżynier twierdzi, że przy pomocy pewnych zabiegów jest w stanie spowodować zwiększenie szansy na urodzenie się krowy powyżej 50%. W jego eksperymencie na 10 urodzonych zwierząt 9 było krowami, a tylko 1 bykiem. Czy powinnniśmy wierzyć temu bio-inżynierowi? Jakia jest szansa na uzyskanie takiego, bądź bardziej ekstremalnego wyniku przy założeniu, że procedura stosowana przez naszego inżyniera nia ma żadnych efektów? W tym problemie dla odmiany założymy, że w normalnych warunkach 100 spośród 206 cieląt to krowy. Zadanie rozwiązać metodą parametryczną i przez repróbkowanie. Wskazówka

*

Odp: p = 0,008. Odrzucamy H₀ o braku efektów.

Zadanie: Porównanie lekarstwa i placebo

Badamy skuteczność leku na raka. Mamy grupę 12 chorych: 6 osobom podajemy lek — poprawa wystąpiła u 5 osób, pozostałym sześciu osobom podajemy placebo — poprawa wystąpiła u 2 osób. Czy te wyniki upoważniają do stwierdzenia, że lek istotnie zwiększa szansę poprawy? Test przeprowadzić na poziomie istotności 5%.

Wskazówka: jako statystykę przyjąć różnicę w ilości popraw między grupą z lekiem a grupą z placebo. Interesuje nas prawdopodobieństwo zaobserwowania takiej (3) bądź większej różnicy.

* 

Odp: Prawdopodobieństwo wylosowania takiej samej bądź większej różnicy w ilości popraw estymowane z symulacji: 0,0702

Prawdopodobieństwo estymowane parametrycznie: 0,0699. Wniosek: nie mamy podstaw do odrzucenia hipotezy zerowej.

Przykład: Średnie grup sparowanych: Lek przeciwdepresyjny

Poniższa tabela prezentuje wyniki 9 pacjentów wykonujących pewien test diagnostyczny przed podaniem leku i po podaniu leku.

przed	po
1,83	0,878
0,50	0,647
1,62	0,598
2,48	2,05
1,68	1,06
1,88	1,29
1,55	1,06
3,06	3,14
1,3	1,29

Skonstruować test, który pozwoli stwierdzić czy lek jest skuteczny. Porównać różne wersje testu:

• bootstrapową (losowanie z powtórzeniami),
• permutacyjną,
• test parametryczny
• test nieparametryczny.

Jakie założenia przyjmujemy przy każdej z wersji testu?

Rozwiązanie

W tym zadaniu mamy dwie grupy przed i po ale oprócz tego istnieje ścisły porządek w parach, bez sensu jest porównywanie przed od jednego pacjenta z po drugiego pacjenta. Musimy stosować testy, które biorą ten porządek pod uwagę (testy pairwise).

Zadania

Zanieczyszczenie środowiska

Agencja ochrony środowiska ustaliła limit na koncentrację zanieczyszczeń emitowanych przez fabryki. Załóżmy, że dopuszczalny poziom zanieczyszczeń wynosi 55 cząstek na milion (cz/m) w promieniu dwóch kilometrów od fabryki. Kontrola przeprowadza 100 pomiarów o różnej porze dnia i roku w promieniu dwóch km od pewnej fabryki. Średnia z próby wyniosła 60 cz/m a odchylenie standardowe [math]s = 20[/math] cz/m. Czy dane te są wystarczające by na poziomie istotności [math]\alpha = 0,01[/math] uznać, że fabryka łamie prawo?

Fabryka łamie prawo jeśli emituje zanieczyszczenia na poziomie wyższym niż dopuszczalny więc należy przeprowadzić test jednostronny (w tym przypadku prawostronny). Czy moglibyśmy odrzucić [math]H_0[/math] na tym samym poziomie stosując test dwustronny? Jest ważne aby w zależności od problemu wybrać odpowiedni test: jedno- lub dwustronny.

Odp. Możemy odrzucić hipotezę [math]H_0[/math] (głoszącą, że fabryka nie łamie prawa) na poziomie [math]\alpha[/math] = 0,01.

* 

Wzrost mocy turbin

Turbina hydroelektryczna generuje moc średnią 25,2 kW. Po unowocześnieniu maszyny chcemy przetestować czy średnia moc generowana zmieniła się (na + lub −). Przeprowadzono 115 pomiarów, które dały średnią 26,1 kW i odch. std. 3,2 kW. Przeprowadzić test statystyczny na poziomie istotności 1%, zinterpretować wynik. Wnioskowanie przeprowadzić także w oparciu oprzedział ufności.

*

Odp.: Poziom p = 0,0026. Przedział: 24,43-25,97.

Sonda

Władze miasta chciałyby wiedzieć, czy przyznać koncesję operatorowi sieci kablowej. W tym celu zleciły nam przeprowadzenie sondy wśród mieszkańców. Zapytaliśmy o zdanie 50 przypadkowo wybranych osób. 30 osób powiedziało „tak” a 20 „nie”. Na ile pewnie otrzymane wyniki wskazują, że mieszkańcy chcą tej kablówki?

Celem naszych badań jest uniknięcie błędu polegającego na tym, że powiemy iż większość mieszkańców chce kablówki podczas gdy tak na prawdę to nie chce.

Wskazówka: Granicznym przypadkiem popełnienia tego błędu jest proporcja 1:1 zwolenników i przeciwników kablówki. Jeśli przeciwników kablówki byłoby jeszcze więcej to uzyskanie naszych wyników byłoby jeszcze mniej prawdopodobne.

*

Wybory prezydenckie

W ankiecie uzyskaliśmy 840 głosów popierających kandydaturę A i 660 kandydaturę B. Jaka jest szansa, że tak naprawdę kandydat B ma poparcie 50% lub większe? Jakie jest prawdopodobieństwo pojawienia sie zaobserwowanej próbki lub próbki wskazującej na jeszcze większe poparcie dla kandydata A, jeśli w rzeczywistości poparcie kandydata A byłoby 50% lub mniej.

*

Czy stosunek do marihuany się zmienił?

Rozważmy dwie ankiety przeprowadzone w USA, pytano 1500 respondentów o stosunek do legalizacji marihuany. Pierwszą ankietę przeprowadzono w 1980, wówczas za legalizacją opowiadało się 52% a drugą w 1985 i za legalizacją było 46%. Czy wyniki tych dwóch ankiet są istotnie różne?

Z jaką proporcją powinniśmy porównywać te wyniki? Jaka jest hipoteza zerowa?

*

Zawały serca i cholesterol

Badano grupę 605 osób. 135 osób z tej grupy miało wysoki poziom cholesterolu a 470 niski. W grupie z wysokim poziomem cholesterolu odnotowano 10 przypadków zawału serca a w grupie z niskim poziomem 21, w czasie 16 lat obserwacji. Nasze pytanie brzmi: Czy możemy uznać, że wysoki poziom cholesterolu zwiększa ryzyko zawału serca?

*

Czy gęstości planet się różnią?

Rozważmy pięć planet znanych w antycznym świecie. Chcemy zbadać, czy planety wewnętrzne Merkury (0,68) i Wenus (0,94) mają istotnie większe gęstości niż planety zewnętrzne Mars (0,71) Jowisz (0,24) i Saturn (0,12)?

Wskazówki:

• Unikalność zestawu planet wskazuje na zastosowanie testu permutacyjnego.
• Moduł implementujący funkcje kombinatoryczne w pythonie to itertools . Zawiera on funkcję permutations. Wywołanie itertools.permutations(sekwencja[, r]) zwraca obiekt permutacji. Obiekt ten zwraca kolejne permutacje o długości r elementów w sekwencji (np. w wektorze). Przykład:

import itertools
for kolejny in itertools.permutations(range(3), 2):
    print kolejny

*