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}} = 31,5 \pm 1,984*\frac{5}{\sqrt{100}} = 31,5 \pm 0,992 = [30,51 \quad 32,49][/math]

• Wnioskowanie: Możemy więc być na 95% pewni, że nieznana średnia leży pomiędzy 30,51 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}} [/math]
• Obliczmy jej wartość dla danych z próby: [math]t = \frac{\bar x - \mu}{s/\sqrt{n}} = (31,5 -28)/(5/\sqrt{100}) = 7 [/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 7 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,'.')

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 Shapiro-Wilka nie jest wrażliwy na parametry rozkładu, z którego pochodzą dane (dane mogą pochodzić z rozkładu o dowolnej średniej i dowolnym odchyleniu standardowym).

Test Kołmogorowa-Smirnowa

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', args=(np.mean(x),np.std(x,ddof=1)) #sposób zalecany
#lub 
D , p = st.kstest((x-np.mean(x))/np.std(x,ddof=1), 'norm') #sposób niezalecany, ale również prawidłowy

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.

UWAGA! Test Kołmogorowa-Smirnowa jest wrażliwy na parametry rozkładu, z którego pochodzą dane. Wywołanie poniższego kodu jest błędem:

D , p = st.kstest(x, 'norm') #BŁĘDNE użycie testu Kołmogorowa-Smirnowa

Test dla innych postaci rozkładu

Za pomocą testu Kołmogorowa-Smirnowa możemy zbadać również założenie o pochodzeniu danych z populacji podlegającej dowolnemu innemu rozkładowi ciągłemu. W tym celu należy podać zamiast 'norm' odpowiednią nazwę rozkładu z modułu scipy.stats oraz parametry tegoż rozkładu w odpowiedniej kolejności (należy w tym celu zajrzeć do dokumentacji). Przykład jak to należy zrobić dla rozkładu wykładniczego:

D , p = st.kstest(x, 'expon', args=(0, 1/np.mean(x)) #test, czy dane podlegają rozkładowi wykładniczemu

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,ddof=1)) 
	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,ddof=1)))
	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', args=(np.mean(x),np.std(x,ddof=1)))
	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(df = 2, size=1000, loc=0, scale=1)
py.figure(2)
porownania(x)
 
x = st.expon.rvs(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.txt')
# 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)/N_rep
p_T  = np.sum( T>=T_0 )/N_rep
p_me = np.sum(me>=me_0)/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)/(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 range(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]

import scipy.stats as st
import numpy as np

'''Przedział ufności dla różnicy dwóch średnich

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ść.

Skonstruujmy przedział ufności 95% dla średniej różnicy w przeżywalności.

Uwaga: przy tym problemie każdą z grup traktujemy jako reprezentantów bardzo 
dużych populacji. '''
def randsample(x,ile):
	ind = st.randint.rvs(0,len(x),size = ile)
	y = x[ind]
	return y

m_sp = np.array([94, 38, 23, 197, 99, 16, 141])
N_sp = len(m_sp)
m_st = np.array([52, 10, 40, 104, 51, 27, 146, 30, 46])
N_st = len(m_st)
alfa = 0.05
# zgodnie z hipotezą zerową (H0) nie ma różnicy między grupami

# symulacja
# tworzymy świat zgodny z H0
m = np.concatenate((m_sp,m_st))

roznica_oryginalna = np.mean(m_sp) - np.mean(m_st)

N_rep = 10000
r = np.zeros(N_rep)
for i in range(N_rep):
	g1 = randsample(m, N_sp)
	g2 = randsample(m, N_st)
	r[i] = np.mean(g1) - np.mean(g2)
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, 'g':ci_g})
print('oryginalna różnica średnich: %(ro).2f'%{'ro':roznica_oryginalna})
# przedział ufności na różnicę skonstruowany zgodnie z H0 zawiera oryginalną różnicę średnich, 
# zatem różnica taka jest na przyjętym poziomie ufności możliwa 
# do zaobserwowania w przypadku braku różnicy między grupami

# zakładając normalność
f = N_sp+N_st-2;
v_1 = np.var(m_sp)
v_2 = np.var(m_st)
sig = np.sqrt( ((N_sp*v_1 +N_st*v_2))/f * (N_sp + N_st)/float(N_sp*N_st) )
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})


print('przedział ufności dla różnicy średnich przy założeniu normalności %(d).2f %(g).2f'%{'d':s_delta*t_2_5+ roznica_oryginalna,'g':s_delta*t_97_5+roznica_oryginalna})
# skonstruowany wokoł oryginalnej różnicy średnich przedział ufności zawiera 0, zatem nie możemy odrzucić możliwości,
# że nie ma różnicy między grupami

Zadanie: 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].

Rozwiązanie: Standardowe importy modułów

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

Zmienne występujące w treści przykładu

mu_0=12
a=0.05
N=144
x=14.6
s=7.8

Test dotyczy średniej wiec obliczamy odchylenie standardowe średniej:

std_mu=s/np.sqrt(N)

Wyjaśnienie |title= Obliczamy wartość statystyki:

 
t=(mu_0-x)/std_mu

Odchylenie standardowe estymowaliśmy z próby. Wartości powyższej statystyki podlegają zatem rozkładowi t. Obliczamy wartości krytyczne odpowiadające poziomowi istotności a. Test jest dwustronny mamy wiec dwie wartości krytyczne: jedna odcina obszar pod funkcja gęstości prawdopodobieństwa na lewo, a druga na prawo od siebie. Pole każdego z tych obszarów wynosi a/2:

 
  
t_kryt_lewy = st.t.ppf(a/2,N-1)
t_kryt_prawy = st.t.ppf( 1-a/2, N-1)

print('obliczona wartość statystyki t: ', t)
print('wartości krytyczne t: %(tl).2f %(tp).2f '%{'tl':t_kryt_lewy, 'tp':t_kryt_prawy})

Możemy też obliczyć prawdopodobieństwo zaobserwowania wartości t takiej jak w naszym zadaniu lub bardziej ekstremalnej:

 
p =  (st.t.cdf(-np.abs(t),N-1)) + (1-st.t.cdf(np.abs(t), N-1))  # sumujemy po obu ogonach  bo test jest dwustronny
print('Prawdopodobieństwo zaobserwowania bardziej ekstremalnych wartości t: %(p).4f'%{'p':p})

Wyniki te możemy zilustrować:

 
os_t = np.arange(-5, 5, 0.1)
py.plot(os_t, st.t.pdf(os_t,N-1)) #rysujemy funkcję gęstości prawdopodobieństwa t o N-1 st. swobody

#cieniujemy lewy obszary pod funkcją gęstości prawdopodobieństwa odpowiadający obliczonemu p
os_t2 = np.arange(-5, t_kryt_lewy, 0.1)
py.fill_between(os_t2,st.t.pdf(os_t2,N-1)) 
#cieniujemy prawy obszary pod funkcją gęstości prawdopodobieństwa odpowiadający obliczonemu p
os_t2 = np.arange(t_kryt_prawy, 5,  0.1)
py.fill_between(os_t2,st.t.pdf(os_t2,N-1)) 

# zaznaczamy obliczoną wartość statystyki:
py.plot((t,), (0,), 'ro')
py.show()

Odpowiedź: Wyliczone t leży poza obszarem akceptacji hipotezy zerowej, zatem odrzucamy hipotezę zerową i akceptujemy alternatywną.

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

Rozwiązanie:

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

mu_0=49;
x=38;
s=14;
N=18;
a=0.01;


# test dotyczy sredniej wiec jej std:
std_mu=s/np.sqrt(N)

# odchylenie std obliczylismy z proby => stosujemy test t
# obliczamy wartość statystyki 
t=(mu_0-x)/std_mu
# obliczamy wartości krytyczne odpowiadające poziomowi istotności a
# test jest dwustronny mamy wiec dwie wartosci krytyczne 
# jedna odcina obszar pod funkcja gęstości prawdopodobieństwa na lewo, 
# a druga na prawo od siebie. Pole każdego z tych obszarów wynosi a/2
  
t_kryt_lewy = st.t.ppf(a/2,N-1)
t_kryt_prawy = st.t.ppf( 1-a/2, N-1)

print('obliczona wartość statystyki t: ', t)
print('wartości krytyczne t: %(tl).2f %(tp).2f '%{'tl':t_kryt_lewy, 'tp':t_kryt_prawy})

# Możemy też obliczyć prawdopodobieństwo zaobserwowania wrtości t takiej jak w naszym zadaniu 
# lub bardziej ekstremalnej:
p = (1-st.t.cdf(t, N-1)) + (st.t.cdf(-t,N-1)) # sumujemy po obu ogonach  bo test jest dwustronny
		
print('Prawdopodobieństwo zaobserwowania bardziej ekstremalnych wrtości t: %(p).3f'%{'p':p})
# Wyniki te możemy zilustrować
os_t = np.arange(-5, 5, 0.1)
py.plot(os_t, st.t.pdf(os_t,N-1)) #rysujemy funkcję gęstości prawdopodobieństwa t o N-1 st. swobody

#cieniujemy lewy obszary pod funkcją gęstości prawdopodobieństwa odpowiadający obliczonemu p
os_t2 = np.arange(-5, t_kryt_lewy, 0.1)
py.fill_between(os_t2,st.t.pdf(os_t2,N-1)) 
#cieniujemy prawy obszary pod funkcją gęstości prawdopodobieństwa odpowiadający obliczonemu p
os_t2 = np.arange(t_kryt_prawy, 5,  0.1)
py.fill_between(os_t2,st.t.pdf(os_t2,N-1)) 

# zaznaczamy obliczoną wartość statystyki:
py.plot(t, 0, 'ro')
py.show()



# 'Odp: Wyliczone t lezy poza obszarem akceptacji hipotezy zerowej,
# zatej odrzucamy hipoteze zerowa i akceptujemy alternatywna.'

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

Rozwiązanie:

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

def randsample(x, N):
	'''zwraca wektor 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

# sposób pierwszy:
# zmienna urodzenie byka/krowy podlega rozkladowi dwumianowemu 
p = 100.0/206
N = 10
k = 9
p_bino = 1-st.binom.cdf(k-1,N,p) #prawdopodobieństwo wylosowania 9 lub 10 krów w 10 probach 
print('Prawdopodobieństwo wylosowania 9 lub 10 krów w 10 probach: %(p).4f'%{'p':p_bino})
# W jego mwetodzie chyba cos jest

# sposob drugi: repróbkowanie
# model swiata z ktorego pochodza byki(0)/krowy(1):
w = np.concatenate((np.ones(100), np.zeros(106)))

N_rep = 100000
wynik = np.zeros(N_rep)
for i in range(N_rep):
	wynik[i] = np.sum(randsample(w,10))

p_rep = float(np.sum(wynik>=k))/N_rep
print('Prawdopodobieństwo wylosowania 9 lub 10 krów w 10 probach estymowane z symulacji: %(p).4f'%{'p':p_rep})

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.

Rozwiązanie:

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

def randsample(x, N):
	'''zwraca wektor 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




# Formułujemy hipotezy 
# H0: lek nie daje poprawy 
# H1: lek daje poprawę



# zgodnie z H0 obie próby pochodzą ze świata:
# pierwszy sposób: repróbkowanie
w = np.concatenate((np.ones(7), np.zeros(5))) # jedynki -> wystąpiła poprawa

n_l = 5 # ilosc popraw w grupie leku
n_p = 2 #ilosc popraw w grupie placebo
# jako statystykę testową przyjmiemy różnicę w ilości popraw miedzy grupami
# w tym problemie istotne jest zwiększenie ilości popraw wiec stosujemy test
# jednostronny
N_rep = 100000

st_0= n_l - n_p
st_rep=np.zeros(N_rep)

for i in range(N_rep):
    n_l_rep   = np.sum(randsample(w,6))
    n_p_rep   = np.sum(randsample(w,6))
    st_rep[i] = n_l_rep - n_p_rep # wartość statystyki w i-tym repróbkowaniu
    


p_rep = float(np.sum(st_rep >= st_0))/N_rep
print('''Prawdopodobieństwo wylosowania takiej samej bądź większej różnicy w ilości popraw estymowane z symulacji: %(p).4f'''%{'p':p_rep})


# drugi sposób:
# zmienna uzyskanie poprawy podlega rozkładowi dwumianowemu 
p = 7.0/12
N = 6
k1 = 5
k2 = 2
p_param = 0
for d in range(N-st_0+1):  # w tej pętli sumujemy prawdopodobieństwa zdarzeń  sprzyjających zaobserwowaniu różnicy co najmniej st_0 popraw 
	p_bino1 = 1-st.binom.cdf(st_0 - 1 +d, N, p)  # prawdopodobieństwo uzyskania poprawy w co najmniej st_0 +d próbach
	p_bino2 = st.binom.pmf(d,N,p)   # prawdopodobieństwo uzyskania poprawy w d próbach
	p_param += p_bino1 * p_bino2 # prawdopodobieństwo zaobserwowania jednocześnie obu powyższych sytuacji

print('Prawdopodobieństwo estymowane parametrycznie: %(p).4f'%{'p':p_param})

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.

Zadanie: Pomiar masy cząstki elementarnej

W pomiarach wstępnych zbadano masę spoczynkową pewnej cząstki elementarnej. Otrzymano następujące wyniki [MeV/c²]:

139.20, 139.34, 140.22, 139.56, 139.42, 139.64, 139.22, 139.74, 139.38, 139.54, 139.38, 139.46, 140.09, 139.77, 139.52, 139.47, 139.89, 138.95, 139.99, 139.64, 139.37, 139.49, 139.15, 139.77, 140.10, 139.48, 139.84, 139.44, 140.13

Zbadaj na poziomie istotności 1%, czy cząstką tą mógł być naładowany pion.

Zbadaj na poziomie istotności 1%, czy cząstką tą mógł być neutralny pion.

Rozwiązanie:

import numpy as np
import scipy.stats as st

alfa = 0.01 # poziom istotności
dane = np.array([139.20, 139.34, 140.22, 139.56, 139.42, 139.64, 139.22, 139.74,
                 139.38, 139.54, 139.38, 139.46, 140.09, 139.77, 139.52, 139.47,
                 139.89, 138.95, 139.99, 139.64, 139.37, 139.49, 139.15, 139.77,
                 140.10, 139.48, 139.84, 139.44, 140.13])

Pi_plus_min = 139.57 # masa pionów naładowanych Pi+ i Pi-
Pi_neutral = 134.98  # masa pionu neutralnego Pi0

# test t hipotezy H0, że dane pochodzą z rozkładu normalnego o wartości oczekiwanej Pi_plus_min
t1, p1 = st.ttest_1samp(dane,Pi_plus_min)
print("Poziom p dla hipotezy, że był do pion naładowany",p1)
if p1>=alfa: print("Nie możemy odrzucić tej hipotezy na poziomie istotności {}%".format(100*alfa))
else: print("Możemy odrzucić tę hipotezę na poziomie istotności {}%".format(100*alfa))

# test t hipotezy H0, że dane pochodzą z rozkładu normalnego o średniej Pi_neutral
t2, p2 = st.ttest_1samp(dane,Pi_neutral)
print("\nPoziom p dla hipotezy, że był do pion neutralny",p2)
if p2>=alfa: print("Nie możemy odrzucić tej hipotezy na poziomie istotności {}%".format(100*alfa))
else: print("Możemy odrzucić tę hipotezę na poziomie istotności {}%".format(100*alfa))

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

PRZED = [1.83, 0.5, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.3] PO = [0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 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).

Jako miarę tego czy lek jest skuteczny przyjmiemy różnicę (po - przed). Każda z zaobserwowanych różnic to zmienna losowa. Ich średnia to nasza STATYSTYKA. Będziemy wierzyli, że lek działa jeśli ta różnica jest istotnie mniejsza od zera. Mamy stąd hipotezy:
[math]H_0: \bar r \geq 0[/math]
[math]H_1: \bar r \lt 0 [/math]

Wersja bootstrapowa

Założenie, które czynimy w wersji bootstrapowej testu jest następujące: Zaobserwowana grupa pacjentów jest reprezentatywna, tzn. istnieje duża populacja pacjentów spośród której moglibyśmy pobrać wiele innych grup podobnych pod względem interesujących nas parametrów do grupy zaobserwowanej w tym konkretnym eksperymencie. Konsekwencją tego założenia jest to, że jeśli w naszej grupie mamy już pacjenta z wynikami: [1.83, 0.878], to szansa na wylosowanie kolejnego pacjenta o takich wynikach się nie zmienia i nadal wynosi 1/9. Prowadzi to do implementacji zawierającej losowanie z powtórzeniami.

Losowość występuje tu w dwóch miejscach:

• losujemy pacjentów z powtórzeniami - zakładamy, że badana grupa jest reprezentatywna dla bardzo dużej populacji.
• dla każdego z wybranych pacjentów losujemy jego wynik przed oraz wynik po z wyników, które uzyskał w rzeczywistości --- to jest zgodne z hipotezą zerową.

Dla każdego losowania obliczamy wartość statystyki w tym losowania uśredniając różnice indywidualne. Wartości statystyki otrzymane dla wszystkich losowań tworzą jej empiryczny rozkład, przy założeniu, że hipoteza zerowa jest prawdziwa. Jako estymator prawdopodobieństwa zaobserwowania oryginalnej wartości statystyki mr (średnia różnica) przy prawdziwej hipotezie zerowej przyjmujemy frakcję rozkładu empirycznego, w której wartości statystyki były nie większe niż mr. Wynik ilustrujemy histogramem rozkładu empirycznego z zaznaczoną ową frakcją.

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

def randsample(x, N):
	'''zwraca wektor 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

def hist_z_markerem(x, N_bins, marker):
	'''Rysuje histogram wartości w tablicy x, używając N_bins binów. 
	Na lewo od wartości wskazanej przez marker dorysowywany jest prostokąt'''
	
	r = np.max(x) - np.min(x)
	szer_binu = r/N_bins

	#konstruujemy biny
        # robimy biny od markera co szerokość binu aż do x minimalnego
	biny_na_lewo = np.arange( marker, np.min(x), -szer_binu)
        # odwracamy kolejność tej sekwencji żeby była rosnąca
	biny_na_lewo = biny_na_lewo[-1::-1] 
         # robimy biny od markera co szerokość binu aż do x maksymalnego
	biny_na_prawo = np.arange(marker,np.max(x), szer_binu)
        # sklejamy oba zakresy binów
	biny = np.concatenate((biny_na_lewo, biny_na_prawo))
	(n,xx,patch) = py.hist(x,bins = biny)
	py.fill([np.min(xx), np.min(xx), marker, marker] , [0, np.max(n), np.max(n), 0] ,'r' ,alpha = 0.2)

A=np.array([[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.30,	1.29]])

r = A[:,1] - A[:,0] # od drugiej kolumny odejmuję pierwszą
mr = np.mean(r)  # średnia z tych różnic to wartość statystyki zaobserwowana dla oryginalnych danych
print( 'średnia różnica: %(mr).2f'%{'mr':mr})

N = len(r)

N_rep = 100000
r_boot = np.zeros(N_rep)
przed = np.zeros(N)
po = np.zeros(N)

for i in range(N_rep):
    ix=randsample(np.arange(0,N,1),N) # wybieramy pacjentów z powtórzeniami
    B=np.array(A[ix,:])
    for j in range(N):   # mieszamy losowo przypisując wyniki do grupy przed i po
                         # zakładając, że pacjent może uzyskać dwukrotnie taki sam wynik
        s = np.random.rand()
        if s > 0.5:
            przed[j]=B[j,0]
        else:
            przed[j]=B[j,1]
        s = np.random.rand()
        if s > 0.5:
            po[j]=B[j,1]
        else:
            po[j]=B[j,0]
    rr=po-przed
    r_boot[i] = np.mean(rr)

hist_z_markerem(r_boot,30,mr)
p_h0 = np.sum( r_boot <=  mr)/N_rep
s_boot = u'dla repróbkowanego testu jednostronnego: %(p_h0).5f'%{'p_h0': p_h0}
print( s_boot)
py.title(s_boot)
py.show()

Wersja permutacyjna

W tym teście zakładamy, że grupa pacjentów jest unikalna, więc w procedurze wytwarzania empirycznego rozkładu statystyki korzystamy z danych wszystkich pacjentów w każdej iteracji.

Zgodnie z hipotezą zerową pomiary przed i po są równoważne można je zatem zamieniać. Wykonamy wszystkie możliwe zamiany przed i po. Możliwych zamian jest [math]2^N[/math]. Skorzystamy z faktu, że bity w reprezentacji binarnej liczb całkowitych od 0 do [math]2^{N-1}[/math] zawierają wszystkie możliwe permutacje ciągów zer i jedynek o długości N. Wartości 1 zamienimy na logiczne True a wartości 0 na False. Zinterpretujemy True jako zamianę i False jako brak zamiany.

Dla każdej permutacji obliczamy wartość statystyki uśredniając różnice indywidualne. Wartości statystyki otrzymane dla wszystkich permutacji tworzą jej empiryczny rozkład, przy założeniu, że hipoteza zerowa jest prawdziwa. Jako estymator prawdopodobieństwa zaobserwowania oryginalnej wartości statystyki mr przy prawdziwej hipotezie zerowej przyjmujemy frakcję rozkładu empirycznego, w której wartości statystyki były nie większe niż mr. Wynik ilustrujemy histogramem rozkładu empirycznego z zaznaczoną ową frakcją.

# -*- coding: utf-8 -*-
import scipy.stats as st
import pylab as py
import numpy as np
 
def dec2bin(n, l):
    '''konwertuje dziesiętną liczbę całkowitą na tablicę 
    przedstawiającą reprezentację binarną tej liczby
    n liczba do konwersji
    l długość reprezentacji binarnej 
    zwracana jest binarna reprezentacja liczby 
    skonwertowana do tablicy logicznej (0->False, 1-> True)
    '''
    b = np.zeros(l, dtype = bool)
    if n < 0:  raise ValueError("must be a positive integer")
    i = 1
    while n > 0:
        b[l-i] = bool( n % 2 ) 
        n = n >> 1
        i += 1
    return b
 
def hist_z_markerem(x, N_bins, marker):
	'''Rysuje histogram wartości w tablicy x, używając N_bins binów. 
	Na lewo od wartości wskazanej przez marker dorysowywany jest prostokąt'''
 
	r = np.max(x) - np.min(x)	
	szer_binu = r/N_bins
 
	#konstruujemy biny
        # robimy biny od markera co szerokość binu aż do x minimalnego
	biny_na_lewo = np.arange( marker, np.min(x), -szer_binu)
        # odwracamy kolejność tej sekwencji żeby była rosnąca
	biny_na_lewo = biny_na_lewo[-1::-1] 
         # robimy biny od markera co szerokość binu aż do x maksymalnego
	biny_na_prawo = np.arange(marker,np.max(x), szer_binu)
        # sklejamy oba zakresy binów
	biny = np.concatenate((biny_na_lewo, biny_na_prawo))
	(n,xx,patch) = py.hist(x,bins = biny )
	py.fill([np.min(xx), np.min(xx), marker, marker] , [0, np.max(n), np.max(n), 0] ,'r' ,alpha = 0.2)
 
A=np.array([[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.30, 1.29]])
 
r = A[:,1] - A[:,0] # od drugiej kolumny odejmuję pierwszą
mr = np.mean(r)  # średnia z tych różnic to wartość statystyki zaobserwowana dla oryginalnych danych
print( 'średnia różnica: %(mr).2f'%{'mr':mr} )
 
N = len(r)
 
 
N_perm = 2**N
r_perm = np.zeros(N_perm)
for i in range(2**N):
	B = np.array(A)   # B zawiera kopię tablicy A
	zamiana = dec2bin(i,N) # w których wierszach dokonać zamiany?
	# print( i,': ', zamiana )
	# wiersze tablicy B wskazane przez True w wektorze indeksów ind zamieniamy wartości  'przed' z wartościami 'po'
	B[zamiana, 0] = A[zamiana, 1] 
	B[zamiana, 1] = A[zamiana, 0]
	rr = B[:,1] - B[:,0] # Obliczam wartości zmiennych losowych dla tej zamiany
	r_perm[i] = np.mean(rr) # Obliczmy wartość statystyki dla tej zamiany
 
 
p_h0 = np.sum( r_perm <=  mr)/N_perm
s_perm = u'dla permutacyjnego testu jednostronnego: %(p_h0).3f'%{'p_h0': p_h0}
print( s_perm)
hist_z_markerem(r_perm,30,mr)
py.title(s_perm)
py.show()

Wersja parametryczna

Jeśli badane różnice przed i po podlegają rozkładowi normalnemu to do testowania czy średnia wartość różnicy jest równa 0 można zastosować test t dla prób zależnych st.ttest_rel(). Aby się upewnić, że możemy zastosować ten test badamy normalność różnic przy pomocy normplot() oraz testu Shapiro-Wilka st.shapiro().

# -*- 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 = (i-0.3175)/(N+0.365)
		y[i] = st.norm.ppf(arg)
	py.plot(y,x_ord,'.')
 
A=np.array([[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.30, 1.29]])

r = A[:,1] - A[:,0] # od drugiej kolumny odejmuję pierwszą

normplot(r)
py.title("Wyniki testów na normalność różnic\n Shapiro-Wilka:  W=%.3f, p=%.3f"%st.shapiro(r)
        +"\nKołmogorowa-Smirnowa: D=%.3f, p=%.3f"%st.kstest(r, 'norm', args=(np.mean(r),np.std(r,ddof=1))))
t, p = st.ttest_rel(A[:,1],A[:,0])
p_t = p/2 # aby test był jednostronny
s_t = u'dla parametrycznego testu jednostronnego: %(p_h0).3f'%{'p_h0': p_t}
print(s_t)
py.show()

Wersja nieparametryczna

W ogólności, nie zakładając normalności różnic można by przeprowadzić test Wilcoxona. W aktualnej implementacji tego testu w scipy.stats jest on dla naszych danych niedokładny, gdyż mamy małą liczebność grupy, a implementacja stosuje przybliżenia asymptotyczne.

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

A=np.array([[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.30, 1.29]])
z, p = st.wilcoxon(A[:,1],A[:,0] )
p_w = p/2 # aby test był jednostronny
s_w = u'dla nieparametrycznego testu jednostronnego: %(p_h0).3f'%{'p_h0': p_w}
print( s_w)

Podsumowanie

średnia różnica: -0.43
dla repróbkowanego testu jednostronnego: 0.0013
dla permutacyjnego testu jednostronnego: 0.014
Wynik testu Shapiro-Wilka na normalność różnic 
	 statystyka W: 0.921 
	 prawdopodobieństwo takiej wartości statystyki dla rozkładu normalnego: 0.404
dla parametrycznego testu jednostronnego: 0.008
dla nieparametrycznego testu jednostronnego: 0.019

Wszystkie testy wskazują prawdopodobieństwo zaobserwowania odpowiadających im statystyk poniżej przyjętego poziomu istotności [math]\alpha = 0.05[/math]. Zatem hipotezę zerową należy odrzucić i przyjąć hipotezę alternatywną.

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.

Rozwiązanie:

# -*- coding: utf-8 -*-
import scipy.stats as st
import numpy as np
 
x_c = 55 #cząstek na milion (cz/m) w promieniu dwóch kilometrów od fabryki. 
#Kontrola przeprowadza 
N = 100 # pomiarów o różnej porze dnia i roku w promieniu dwóch km. od pewnej fabryki. 
#średnia z próby wyniosła 
x_s = 60 # cz/m 
# a odchylenie standardowe 
s = 20 # cz/m. 
# Czy dane te są wystarczające by uznać, że fabryka łamie prawo ?
alpha = 0.01
 
# H0: mu <=x_c
# H1: mu > x_c
# wystarczy sprawdzić jakie jest p dla największego dopuszczalnego stężenia 
# mu = x_c zaobserwowania wartości średniej większej bądź równej zaobserwowanej x_s. 
# Odchylenie std. dane jest dla populacji i trzeba je przeliczyć na odchylenie std. średniej.
p = 1 - st.t.cdf(x_s, df=N-1, loc = x_c, scale = s/np.sqrt(N) )
print("poziom p",p)


############################
#inne (tożsame) rozwiązania:

##wykorzystaniem "Survival function" zamiast dystrybuanty
p = st.t.sf(x_s, df=N-1, loc = x_c, scale = s/np.sqrt(N) )
print("poziom p",p)

##obliczając statystykę t
t=(x_s-x_c)/(s/N**0.5)
p = 1-st.t.cdf(t,df=N-1)
print("poziom p",p)

Odp. [math]p=0,007\lt \alpha[/math], zatem 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.

Rozwiązanie:

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


N=115
m=25.2	#hipoteza zerowa mówi, że tyle wynosi moc
sr=26.1	#w pomiarach uzyskaliśmy taką średnią moc
s=3.2	#przy takim odchyleniu standardowym

alfa=0.01

#obliczenie poziomu p
t=(m-sr)/(s/N**0.5)
p = st.t.cdf(t,df=N-1) * 2 #mnozymy przez 2, bo chcemy testu dwustronnego
print("poziom p",p)

#alternatywnie - obliczenie przedziału ufności
lo,hi = st.t.ppf([alfa/2, 1-alfa/2],df=N-1,loc=sr,scale=s/N**0.5)
print("przedzial ufnosci [%.3g,%.3g]"%(lo,hi))

""" Odp.: [math]p=0,\!003\lt \alpha=0,\!01[/math], zatem odrzucamy H0.

Odp. (alt.): [math]\mu = 25,\!2 \not\in[25,\!3; 26,\!9][/math], zatem odrzucamy H0.

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.

Rozwiązanie:

import numpy as np
import scipy.stats as st
import numpy.random as rnd

N=50
odp_na_tak=30
odp_na_nie=20

#przypadek graniczny, który daje najwyższe prawdopodobieństwo uzyskania fałszywego wyniku na 'tak',
#to mieszkańcy niezdecydowani, czyli propocja 1:1
#(propocja przechylona w kierunki 'nie', da nam niższe oszacowanie
#a propocja przechylona w kierunku 'tak' nie oznaczałaby fałszywego wyniku sondy)
#zatem hipoteza zerowa brzmi, że nasz wynik pochodzi z populacji niezdecydowanej:
P=0.5 #prawdopodobieństwo, że respondent udzieli odpowiedzi na 'tak'

#rozwiązanie parametryczne wymaga użycia rozkładu dwumianowego:
p = 1 - st.binom.cdf(29,50,P) #od 1 odejmujemy sumaryczne prawdopodobieństwo wszystkich przypadków udzielenia mniej niż 30 odpowiedzi na tak (czyli do 29 włącznie)

print("poziom p (parametrycznie)",p)


#rozwiązanie przez symulację:
Nrand=100000
R=0 #rezultat, czyli liczba przypadków, w których przynajmniej 30 osób odpowiedziało tak
for i in range(Nrand):
	odpowiedzi = (rnd.random(N)<P) #losujemy odpowiedzi od N osób
	if odpowiedzi.sum() >=odp_na_tak: #jeśli zdażyło się co najmniej 30 odpowiedzi na tak, to dodajemy jeden do liczby takich przypadków
		R+=1
	
print("poziom p (z symulacji)",R/Nrand) #na koniec dzielimy przez liczbę wszystkich sond w symulacji

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.

Rozwiązanie:

import numpy as np
import scipy.stats as st
import numpy.random as rnd

# rozwiązanie parametryczne
A = 840
B = 660
N = A+B

#przypadek graniczny, który daje najwyższe prawdopodobieństwo uzyskania fałszywego wyniku na 'tak',
#to mieszkańcy niezdecydowani, czyli propocja 1:1
#propocja przechylona w kierunki 'nie', da nam niższe oszacowanie
#a propocja przechylona w kierunku 'tak' nie oznaczałaby fałszywego wyniku sondy
#zatem hipoteza zerowa brzmi, że nasz wynik pochodzi z populacji niezdecydowanej:
p_sukces = 0.5

#rozwiązanie parametryczne
p = 1 - st.binom.cdf(A-1,N,p_sukces)

print("poziom p (parametrycznie)",p)


#rozwiązanie przez symulację:

#wiemy już, że musimy oszacować bardzo niskie prawdopodobieństwo, zatem symulacja musi być długa
#(gdybyśmy nie wiedzieli, to łatwo byśmy zauważyli, że krótka symulacja nie daje w ogóle zdarzeń sprzyjających dając wynik p=0)
#niestety trzeba pokusić się o optymalizację (mimo to kod będzie się wykonywał kilkadziesiąt sekund do kilku minut):
#chemy mieć przynajmniej 10^7 sond, ale każda sonda to 1500 odpowiedzi, zatem łącznie musimy wylosować 1.5*10^10 liczb
#losowanie w pętli 10^7 razy kolejnych sond po 1500 odpowiedzi zajmie bardzo bardzo dużo czasu
#tablica 1.5*10^10 odpowiedzi wygenerowana jednym poleceniem nie zmieści się w pamięci
#musimy skorzystać z rozwiązania pośredniego
#będziemy losować 100 razy po 1.5*10^8 odpowiedzi i dodamy wyniki

Nrand_1=100    # 10^2
Nrand_2=100000 # 10^5
N = A+B        # 1.5*10^3

Nrand_total = Nrand_1*Nrand_2

R=0 #liczba sond, w których kandydat A uzyskał poparcie 840 lub większe
#po każdym przebiegu pętli będziemy dodawać liczbę takich przypadków (jeśli wystąpiły w danej iteracji)

for i in range(Nrand_1):
	#w każdej iteracji losujemy odpowiedzi od N osób w Nrand_2 sondach
	#iteracji będzie w sumie Nrand_1, co da nam łącznie N osób w Nrand_total sondach
	#jako, że zakładamy równe prawdopodobieństwo, to może możemy wylosować tablicę zawierającą tylko 0 i 1
	#gdzie 1 oznacza sukces (głos na kandydata A), a 0 porażkę
	#korzystając z funksji randint, która zwraca (pseudo)losowe liczby całkowite z zakresu [a,b)
	a,b = 0,2 
	odpowiedzi = rnd.randint(a,b,size=(N,Nrand_2)) 
	S = np.sum(odpowiedzi,axis=0)  #wyniki kolejnych sond (sumujemy sukcesy w wierszach)

	R+=np.sum(S>=A) #dodajemy liczbę sond, w których kandydat A uzyskał poparcie 840 lub większe (w danej iteracji pętli)

print("poziom p (z symulacji)",R/Nrand_total) #na koniec dzielimy przez liczbę wszystkich sond w symulacji

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?

Rozwiązanie:

import scipy.stats as st
import numpy as np

# rozwiązanie parametryczne
za1 = int(0.52*1500)
za2 = int(0.46*1500)
N = 1500

# H0: stosunek się nie zmienił zatem p_za było w obu anlietach takie samo:
p_za = (0.52+0.46)/2

# Jakie jest prawdopodobieństwo zaobserwowoania wyników ankiet o takiej bądź większej 
# różnicy głosów za, jeśli obie są zgodne z H0 ?
D = za1 - za2

#musimy zsumować wszystkie możliwości, w których wystąpiła taka różnica lub większa
#zrobimy to w pętli
p_binom = 0.0
for k in range(N-D+1):
	p1 = st.binom.pmf(k,N,p_za) # prawdopodobieństwo uzyskania dokładnie k głosów za
	
	# sumujemy prawdopodobieństwa uzyskania liczby głosów różnej przynajmniej o D
	przypadki_o_roznicy_przynajmniej_D = np.arange(k+D,N+1)
	p2 = np.sum(st.binom.pmf(przypadki_o_roznicy_przynajmniej_D,N,p_za))
	#co jest równoważne:
	#p2 = 1-st.binom.cdf(k+D-1,N,p_za)
	
	#sumujemy dla kolejnych k
	p_binom += p1*p2 # mnożymy prawdopodobieństwa ponieważ są wyniki obu sond są niezależne
	
p=2*p_binom #mnożymy przez 2, ponieważ robimy test dwustronny (różnica mogła wystąpić również w drugą stronę)
print("p =",p)

# rozwiązanie symulacyjne (tym razem skorzystamy z randsample)
def randsample(x, N):
	'''zwraca wektor 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

N_rep = 100000
# świat zgodny z H0:
H0 = np.concatenate((np.zeros(N*(1-p_za)), np.ones(N*p_za)))

a = np.zeros(N_rep)
for i in range(N_rep):
	w1 = np.sum(randsample(H0,N))
	w2 = np.sum(randsample(H0,N))
	a[i] = np.abs(w1-w2)
p = np.sum(a>=D)/N_rep
print(p)

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?

Rozwiązanie:

import scipy.stats as st
import numpy as np

# rozwiązanie parametryczne
wysoki = 135
niski = 470
N = 605
wysoki_zawaly = 10
niski_zawaly = 21
N_zaw = wysoki_zawaly + niski_zawaly 

# H0: proporcja zawałowców w obu grupach taka sama:
p_zawalu = N_zaw/N
# H1: proporcja jest różna, test jednostronny (wysoki cholesterol ZWIĘKSZA ryzyko)

# Jakie jest prawdopodobieństwo zaobserwowania wyników badania o takiej bądź większej 
# różnicy proporcji zawałów, jeśli obie są zgodne z H0 ?
roznica = wysoki_zawaly/wysoki - niski_zawaly/niski

p_binom = 0.0
for k_wys in range(wysoki+1):
	p_wys = st.binom.pmf(k_wys,wysoki,p_zawalu) # prawdopodobieństwo uzyskania k_wys zawałów grupie o liczebności wysoki
	for k_nis in range(niski+1): 
		if k_wys/wysoki - k_nis/niski>=roznica: # czy liczebności k_wys i k_nis dają większą bądź równą różnice proporcji?
			p_nis = st.binom.pmf(k_nis,niski,p_zawalu) # prawdopodobieństwo uzyskania k_nis zawałów w grupie o niskim cholesterolu
			p_binom += p_wys*p_nis # zdarzenia są niezależne

print(p_binom)

# rozwiązanie symulacyjne
def randsample(x, N):
	'''zwraca wektor 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

N_rep = 100000
# świat zgodny z H0:
H0 = np.concatenate((np.zeros(N-N_zaw), np.ones(N_zaw)))

wyn = np.zeros(N_rep)
for i in range(N_rep):
	wysoki_zawaly_r = np.sum(randsample(H0,wysoki))
	niski_zawaly_r = np.sum(randsample(H0,niski))
	wyn[i] = wysoki_zawaly_r/wysoki - niski_zawaly_r/niski
p = np.sum(wyn>=roznica)/N_rep
print(p)

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)

Rozwiązanie:

import numpy as np
import itertools as itt
PLANETY = np.array([0.68,0.94,0.71,0.24,0.12])
N=len(PLANETY)

r=[] #używamy listy, bo nie chcemy obliczać ile będzie przypadków (aczkolwiek można to zrobić)

for planety in itt.permutations(PLANETY):
	r.append(np.mean(planety[:2])-np.mean(planety[2:]))

R=np.mean(PLANETY[:2])-np.mean(PLANETY[2:]) #prawdziwa różnica

r=np.array(r) #porównanie z liczbą działa tylko dla tablic numpy, a nie dla list

print("p =",np.sum(r>=R)/len(r))

Elektrownia jądrowa

Przed wybudowaniem elektrowni jądrowej przeprowadzono pomiary intensywności promieniowania jonizującego w pobliżu planowanego budynku reaktora. Powtórzono te pomiary po uruchomieniu reaktora. Czy zebrane dane pozwalają stwierdzić, że poziom promieniowania istotnie wzrósł? Dane pomiarowe wczytaj z pliku.

Rozwiązanie:

import numpy as np
import scipy.stats as st

alfa = 0.01 # poziom istotności

PRZED, PO = np.loadtxt('Pomiary_skazen.txt') #wczytanie danych (dwa zbiory przed uruchomieniem i po)

t, p = st.ttest_ind(PRZED, PO) # wykorzystanie gotowej funkcji wykonującej DWUSTRONNY test t na różnicę średnich
                            # dla danych NIESPAROWANYCH (niezależnych)
                            # ttest_ind  - ind od independent 
                            # ttest_rel  - rel od related 

p/=2 # dzielimy przez dwa ponieważ chcemy zrobić test jednostronny
print("Poziom p",p)

Odp. [math]p=0,\!39\gt \alpha=0,\!01[/math], zatem nie mamy podstaw do odrzucenia hipotezy [math]H_0[/math], że poziom promieniowania istotnie wzrósł.