WnioskowanieStatystyczne/ Przedziały ufności

Wstęp

Średnia (lub inny parametr) estymowany na podstawie próby pobranej z populacji nie jest równa rzeczywistej wartości tego parametru w populacji.

• Wylosuj 10 razy 100 liczb z rozkładu normalnego [math]N(0,1)[/math] i oblicz średnią. Czy udało ci się zaobserwować wartość 0?

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

for i in range(10):
	x=st.norm.rvs(size=100)
	print(np.mean(x))

Rozbieżność między uzyskanym wynikiem a rzeczywistą średnią populacji zależy od wielkości badanej grupy oraz zmienności badanej cechy w jej obrębie. Jeśli badana grupa jest niewielka i ma dużą zmienność analizowanej cechy wówczas rozbieżność między średnią uzyskaną, a rzeczywistą może być znaczna. Natomiast, jeśli badana grupa jest duża z niewielką zmiennością danych, wówczas uzyskana średnia będzie prawdopodobnie bardzo bliska średniej populacyjnej. Przedział ufności ((ang. confidence interval — CI)) odzwierciedla zarówno wielkość badanej grupy jak i zmienność analizowanej cechy wewnątrz tej grupy. Przedział ufności jest określany z różnym procentem „zaufania”, np. 90 czy też 95%. Najczęściej używa się 95% przedziału ufności. Oznacza to, że jeśli grupa badana była zgromadzona w sposób losowy to rzeczywisty parametr populacji z 95% pewnością znajduje się w tym przedziale. Innymi słowy, możemy sobie wyobrazić, że mamy dostępne wyniki dowolnie wielu badań analogicznych do analizowanego. Oczekujemy, że w 95% przypadków badana cecha znajduje się wewnątrz wyznaczonego przez nas przedziału ufności.

Pewnych problemów koncepcyjnych nastręcza konstrukcja przedziałów ufności. Dla danej znanej populacji możemy obliczyć średnią populacji oraz prawdopodobieństwo uzyskania konkretnej wartości średniej przy losowaniu próby o zadanej liczebności z tejże populacji, możemy więc określić prawdopodobieństwo [math]P[/math], że odległość średniej z próby i średniej z populacji jest [math]D[/math]. Mając do dyspozycji tylko próbę możemy albo założyć, że pochodzi ona z jakiegoś znanego rozkładu i wyliczyć z niego wartości krytyczne albo zakładając, że jest ona reprezentatywna możemy metodą bootstrapu „wytworzyć” wiele innych prób z badanej populacji i oszacować jakie są granice, w które wpada żądana frakcja średnich (np.:90%, 95%).

Bootstrap jest związany z pobieraniem próby. Najkorzystniejszą sytuacją jest ta, w której dla oszacowania różnych parametrów statystycznych populacji mamy możliwość pobierania z tej populacji wielu prób. Jeśli jest to niemożliwe możemy posłużyć się pobieraniem wielokrotnie prób z tej próby, którą posiadamy. Postępowanie takie jest sensowne pod warunkiem, że próba, która służy nam do generowania innych możliwych pobrań próby jest reprezentatywna. W bootstrapie losujemy ze zwracaniem (dlaczego?).

Przedział ufności dla średniej

Przypadek 1: znana wariancja rozkładu

Rozważamy zmienną losową podlegającą rozkładowi [math]N(\mu,\sigma^2)[/math] czyli rozkładowi normalnemu o średniej [math] \mu[/math] i znanym odchyleniu standardowym [math] \sigma[/math]. Pobieramy próby tych zmiennych losowych o rozmiarze [math]n[/math]. Przypomnijmy, że średnia z takiej próby też może być traktowana jako zmienna losowa. Podlega ona rozkładowi normalnemu o średniej [math] \mu[/math] i standardowym odchyleniu [math]\frac{\sigma}{\sqrt{n}}[/math] czyli [math]N\left( \mu,\frac{\sigma^2}{n} \right)[/math].

Oznacza to, że zmienna [math]\frac{\bar x - \mu}{\sigma/\sqrt{n}} [/math] podlega rozkładowi [math] N(0,1)[/math] czyli [math] Z[/math]. Konstrukcja przedziału ufności dla średniej polega zatem na znalezieniu takich wartości [math]z_{\alpha/2}[/math], że:

[math] P \left( \left|\frac{\bar x - \mu}{\sigma/\sqrt{n}} \right| \le z_{\alpha/2}\right) = 1-\alpha [/math]

Zatem przedział ufności [math] (1-\alpha)*100[/math]% :

[math] \left[ \bar x - z_{\alpha/2}\frac{\sigma}{\sqrt n} ,\;\;\; \bar x+ z_{\alpha/2}\frac{\sigma}{\sqrt n} \right] [/math].

Przykład

Wyciągamy losową próbę [math] (n = 25)[/math] z populacji o rozkładzie normalnym. Dostajemy średnią z próby [math] \bar x = 122[/math]. Załóżmy, że znamy standardowe odchylenie populacji [math] \sigma = 20[/math].

• Oblicz przedział ufności 95% dla średniej populacji [math] \mu[/math].
• Co zrobić aby zmniejszyć obliczony przedział 10-krotnie?

Korzystając z tablic znajdujemy [math] z_{\alpha/2} = z_{2,5%} = -1,96[/math].

Wartość tą możemy uzyskać wołając w pythonie:

st.norm.ppf(0.025)

[math] \bar x \pm 1,96 \frac{\sigma}{\sqrt{n}} = 122 \pm 1,96\cdot\frac{20}{\sqrt{25}} = 122 \pm 7,84 = [\mathrm{114,16} \;\;\; 129,84][/math]

Możemy być pewni na 95%, że nieznana średnia populacji [math] \mu[/math] znajduje się pomiędzy 114,16 a 129,84. Jeśli chcemy zmniejszyć przedział ufności 10-krotnie, musimy pobrać 100 razy większą próbę tj. n = 2500.

Pracując z modułem scipy.stats mamy do dyspozycji rozkłady normalne [math]N(\mu,\sigma)[/math] o dowolnych parametrach [math]\mu[/math] i [math]\sigma[/math].

Możemy więc uzyskać interesujący nas przedział ufności przy pomocy pojedynczego wywołania:

st.norm.ppf([0.025, 0.975],loc = 122, scale = 20.0/25**0.5)

gdzie loc = [math]\mu[/math] oraz scale = [math]\sigma[/math].

Zadanie

Importer win musi zbadać średnią zawartość alkoholu w nowej partii win francuskich. Z doświadczenia z poprzednimi gatunkami wina, przyjmuje on, że standardowe odchylenie w populacji wynosi 1,2%. Importer wybrał losową próbę 60 butelek nowego wina i otrzymał średnią z próby 9,3%. Znaleźć przedział ufności 90% dla średniej zawartości alkoholu w nowej partii win.

Odp.: [9,0452 9,5548].

import scipy.stats as st

sigma=1.2
N=60
x_sr=9.3
sigma_sr=sigma/N**0.5
alfa=0.1

#korzystając ze wzoru
print(st.norm.ppf([alfa/2,1-alfa/2]) * sigma_sr + x_sr)
#lub korzystając z rozkładu o zadanych parametrach
print(st.norm.ppf([alfa/2,1-alfa/2], loc=x_sr, scale=sigma_sr) )

Przypadek 2: nieznana wariancja rozkładu

Sytuacja zmienia się jeśli zakładamy, że próba pochodzi z rozkładu normalnego o śreniej [math]\mu[/math] i nieznanej wariancji. Musimy wówczas zamiast wariancji posłużyć się jej estymatorem [math]s^2 = \frac{1}{n-1} \sum_{i=1}^n {\left(x_i-\mu\right)^2}[/math]. Zmienna losowa [math]\frac{\bar x - \mu}{s/\sqrt{n}} [/math] podlega rozkładowi [math]t_{n-1}[/math] (czytaj: rozkładowi t o n−1 stopniach swobody). Konstrukcja przedziału ufności dla średniej polega zatem na znalezieniu takich wartości [math]t_{\alpha/2}[/math] w rozkładzie t o n−1 stopniach swobody, że:

[math] P \left( \left|\frac{\bar x - \mu}{s/\sqrt{n}} \right| \le t_{\alpha/2}\right) = 1-\alpha [/math]

Zatem przedział ufności [math] (1-\alpha)\cdot 100[/math]% dla średniej [math]\mu[/math], gdy nie znamy odchylenia standardowego [math]\sigma[/math], a jedynie jego estymatę:

[math] \left[ \bar x- t_{\frac{\alpha}{2}} \frac{s}{\sqrt{n}}, \,\,\,\,\, \bar x+ t_{\frac{\alpha}{2}} \frac{s}{\sqrt{n}}\right] [/math]

gdzie [math] t_{\alpha/2}[/math] jest wartością, która odcina obszar α/2 z rozkładu [math] t[/math] z n−1 stopniami swobody.

Uwaga

w Pythonie możemy skorzystać z funkcji modułu numpy do obliczania wartości estymatorów:

• s numpy.std
• [math]\bar x[/math] numpy.mean

Przykład

Lekarz chce zbadać średni czas trwania kuracji tj. od podania leku do ustąpienia objawów w pewnej chorobie. Losowa próba 15 pacjentów dała średni czas [math]\bar x = 10,37[/math] dnia i odchylenie standardowe [math]s = 3,5[/math] dnia. Zakładając normalny rozkład w populacji czasów trwania kuracji znaleźć 95% przedział ufności dla średniego czasu trwania kuracji.

Znajdujemy wartość z rozkładu [math] t[/math] o n−1 (=14) stopniach swobody, która odcina obszar [math] \alpha/2 = 0,025[/math]. [math] t_{0,025} = 2,145[/math].

w pythonie:

st.t.ppf(0.025,14)

Dostajemy więc [math] \bar x \pm t_{\alpha/2}\frac{s}{\sqrt{n}} = 10,37 \pm 2,145\cdot\frac{3,5}{\sqrt{15}} = [8,\!43\ 12,\!31][/math].

Podobnie jak w poprzednim przykładzie możemy też skorzystać z tego, że w scipy.stats definiowane są rozkłady [math]t[/math] o dowolnych parametrach i wykonać całe obliczenia jednym wywołaniem funkcji:

st.t.ppf([0.025, 0.975], 14, loc = 10.37, scale = 3.5/15**0.5)

Lekarz może być pewny, że w 95% przypadków od podania leku do ustąpienia objawów upłynie czas pomiędzy 8,43 a 12,31 dnia.

Zadanie

Producent opon rowerowych chce oszacować średni dystans jaki można przejechać na oponie pewnego rodzaju zanim opona się zużyje. Pobrano losową próbę 32 opon, opona jest używana aż do przetarcia i odległość przejechana na każdej oponie jest rejestrowana. Dane (w tysiącach kilometrów) są następujące:

32, 33, 28, 37, 29, 30, 25, 27, 39, 40, 26, 26, 27, 30, 25, 30, 31, 29, 24, 36, 25, 37, 37, 20, 22, 35, 23, 28, 30, 36, 40, 41

Znaleźć 99% przedział ufności dla średniego przebiegu opon tego rodzaju. Zadanie rozwiązać analogicznie jak powyższy przykład.

Odp: [ 27,76 33,36].

import scipy.stats as st
import numpy as np

X=np.array([32, 33, 28, 37, 29, 30, 25, 27, 39, 40, 26, 26, 27, 30, 25, 30, 31, 29, 24, 36, 25, 37, 37, 20, 22, 35, 23, 28, 30, 36, 40, 41])
alfa=0.01

N=len(X)
x_sr = X.mean()
s=X.std(ddof=1)
s_sr=s/N**0.5

print(st.t.ppf([alfa/2,1-alfa/2],df=N-1,loc=x_sr,scale=s_sr))

-->

Zadanie to można rozwiązać także w oparciu o bootstrap, wykorzystując funkcję zwracającą wartość zadanego kwantyla dla danego rozkładu. Rozkład jest tu zadany empirycznie w postaci zbioru.

 st.scoreatpercentile(zbiór, kwantyl)

Proszę zwrócić uwagę na rozbieżności w wynikach. Z czego one wynikają?

Rozmiar próby

Gdy pobieramy próbę, często chcielibyśmy znać minimalny rozmiar próby, który zapewni nam żądaną precyzję wyniku.

Musimy odpowiedzieć sobie na trzy pytania:

• Jak nasze oszacowanie nieznanego parametru ma być bliskie prawdziwej wartości? Odpowiedź oznaczmy [math] D[/math] (dystans).
• Jaki chcemy mieć poziom ufności, że nasze oszacowanie i prawdziwa wartość parametru są od siebie oddalone o nie więcej niż D?
• Jakie jest oszacowanie wariancji w populacji?

Jeśli nie znamy odpowiedzi na pkt. 3 przeprowadzamy tzw. pilot study i szacujemy [math] \sigma[/math] na podstawie odchylenia std. próby.

Średnia [math]\bar x[/math] podlega rozkładowi normalnemu [math]N(\mu,\sigma^2/\sqrt{n})[/math]. Wymagana odległość pomiędzy [math]\bar x -\mu =D[/math]. Korzystając z transformacji do rozkładu standardowego Z możemy zapisać [math]\bar x = \mu +z \sigma/\sqrt{n}[/math]. Podstawiając do poprzedniego wyrażenia otrzymujemy minimalny rozmiar próby potrzebny do oszacowania średniej populacji [math] \mu[/math], który wynosi:

[math] n=\frac{z_{\alpha/2}^2\sigma^2}{D^2} [/math]

W przypadku nieznajomości [math]\sigma[/math] korzystamy z jej nieobciążonego estymatora [math]s[/math].

Zadanie

Biuro podróży chce oszacować średnią ilość pieniędzy wydaną na wakacje przez osoby korzystające z jego usług. Ludzie przeprowadzający analizę chcieliby móc oszacować średni koszt wakacji z dokładnością do 200 zł na poziomie ufności 95%. Z poprzednich doświadczeń tego biura podróży wynika, że odchylenie standardowe w populacji wynosi [math] \sigma = 400[/math] zł. Jaka będzie minimalna wielkość próby?

Odp: [math] n = 15,366[/math] więc wielkość próby wynosi 16 (zaokrąglamy w górę).

import numpy as np
import scipy.stats as st
import pylab as py
 
def randsample(x,ile):
	ind = st.randint.rvs(0,len(x),size = ile)
	y = x[ind]
	return y
 
alfa = 0.05
x=np.ones(100)
x[56:]=0
Nboot=100000
A=np.zeros(Nboot)
for i in range(Nboot):
    ankieta=randsample(x,1500)
    A[i]=np.sum(ankieta)/1500.0
lo = st.scoreatpercentile(A, per = alfa/2*100)
hi = st.scoreatpercentile(A, per = (1-alfa/2)*100)
print('przedzial ufnosci: %(lo).3f - %(hi).3f\n'%{'lo':lo,'hi':hi})
szer_binu = (hi-lo)/10
biny = np.arange(lo-10*szer_binu, hi+10*szer_binu, szer_binu)
(n,x,patch) = py.hist(A,bins = biny )
py.plot([lo, lo] , [0, np.max(n)] ,'r' )
py.plot([hi,hi],[0, np.max(n)],'r')
py.show()

Przykład z bootstrapem

Rozważmy sondę przedwyborczą, mamy dwóch kandydatów na prezydenta. Ankietowano 1500 osób. 840 osób deklarowało poparcie dla kandydata A zaś 660 dla kandydata B. Na ile pewny może być kandydat A swojego zwycięstwa?

• Jak dokładnie brzmi pytanie? W terminologii przedziałów ufności możemy je sformułować następująco: Jaki jest 95% przedział ufności dla poparcia kandydata A w całej populacji? Czy też innymi słowami: W jakim przedziale na 95% znajduje się proporcja glosujących popierających kandydata A.
• Nasze najlepsze mniemanie o własnościach „świata” z którego pochodzą dane otrzymujemy ze zwykłej proporcji. Wynika z niej, że kandydat A ma poparcie 56% zaś kandydat B poparcie 44% wyborców.
• Przypiszmy do kandydata A „1” zaś do B „0” (w ten sposób tworzymy zmienną losową: oddany głos).
• Pobranie ankiety modelujemy przez pobranie losowo 1500 próbek z modelu naszego „świata” czyli wektora złożonego z 56 jedynek i 44 zer. Wynikiem jednej ankiety jest proporcja popierających kandydata A (lub B)
• Zbieramy rozkład proporcji - musimy w tym celu „przeprowadzić” wielokrotnie ankietę. Narysujmy histogram.
• Chcemy znaleźć 95% przedział ufności musimy znaleźć kwantyl 2,5 oraz 97,5.

Liczby te stanowią poszukiwany przedział ufności.

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

def randsample(x,ile):
	ind = st.randint.rvs(0,len(x),size = ile)
	y = x[ind]
	return y

alfa = 0.05
x=np.ones(100)
x[56:]=0
Nboot=10000
A=np.zeros(Nboot)
for i in range(Nboot):
    ankieta=randsample(x,1500)
    A[i]=np.sum(ankieta)/1500.0
lo = st.scoreatpercentile(A, per = alfa/2*100)
hi = st.scoreatpercentile(A, per = (1-alfa/2)*100)
print('przedzial ufnosci: %(lo).3f - %(hi).3f\n'%{'lo':lo,'hi':hi})
szer_binu = (hi-lo)/10
biny = np.arange(lo-10*szer_binu, hi+11*szer_binu, szer_binu)
(n,x,patch) = py.hist(A,bins = biny )
py.plot([lo, lo] , [0, np.max(n)] ,'r' )
py.plot([hi,hi],[0, np.max(n)],'r')
py.show()

Zadania

Przyrost masy w nowej diecie

Producent karmy dla zwierząt chciał przetestować nowy rodzaj karmy. Próbki podawał 12 zwierzakom przez 4 tygodnie. Po tym czasie zanotował następujące przyrosty masy:

15.43, 16.92, 14.43, 12.94, 15.92, 17.42, 18.91, 16.92, 14.93, 14.49, 15.92, 15.43 kg

średni przyrost wynosi 15.80 kg. Producent widzi jednak, że w próbie jest dość znaczny rozrzut pomiędzy poszczególnymi zwierzętami 12,94-18,91 i nie jest pewien czy można reklamować nowy produkt podając średni przyrost 15,8 kg. Podejrzewa, że inna grupa zwierząt może mieć zupełnie inną średnią.

• Używając powyższych danych znajdź 95% przedział ufności na średni przyrost masy.
• Wynik zilustruj przy pomocy histogramu.
• Jaki byłby wynik przy założeniu, że masy zwierząt pochodzą z rozkładu normalnego?

Odp:

zakładając, że x pochodzi z rozkładu normalnego:

[14,80 16,81];

bootstrap: 14,94-16,67.

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

def randsample(x,ile):
	ind = st.randint.rvs(0,len(x),size = ile)
	y = x[ind]
	return y

X=np.array([15.43, 16.92, 14.43, 12.94, 15.92, 17.42, 18.91, 16.92, 14.93, 14.49, 15.92, 15.43])

alfa = 0.05
Nboot=int(1e5)

Xsr=X.mean()
s=X.std(ddof=1)
N=len(X)
print("średnia",Xsr)

M=np.zeros(Nboot)
for i in range(Nboot):
	xb=randsample(X,N)
	M[i]=xb.mean()

lo,hi = st.t.ppf([alfa/2,1-alfa/2],df=N-1,loc=Xsr,scale=s/N**0.5)
print('przedzial ufnosci z założenia o normalności: %(lo).3f - %(hi).3f\n'%{'lo':lo,'hi':hi})

lo = st.scoreatpercentile(M, per = alfa/2*100)
hi = st.scoreatpercentile(M, per = (1-alfa/2)*100)
print('przedzial ufnosci z bootstrapu: %(lo).3f - %(hi).3f\n'%{'lo':lo,'hi':hi})

szer_binu = (hi-lo)/10
biny = np.arange(lo-10*szer_binu, hi+10*szer_binu, szer_binu)
(n,x,patch)=py.hist(M,bins=biny)
py.plot([lo, lo] , [0, np.max(n)] ,'r' )
py.plot([hi,hi],[0, np.max(n)],'r')

py.show()

Zawartość aluminium w Tebańskich naczyniach.

Zawartość procentowa aluminium w 18 antycznych naczyniach z Teb była następująca:
11.4, 13.4, 13.5, 13.8, 13.9, 14.4, 14.5, 15, 15.1, 15.8, 16, 16.3, 16.5, 16.9, 17, 17.2, 17.5, 19.0
Jaka jest mediana procentowej zawartości aluminium i jaki jest 95% przedział ufności?

Odp:

mediana 15,45 i przedział ufności: 14,20 16,70.

Średnica drzew

Ogrodnik eksperymentuje z nowym gatunkiem drzew. Posadził 20 sztuk i po dwóch latach zmierzył następujące średnice pni (w cm): 8.5, 7.6, 9.3, 5.5, 11.4, 6.9, 6.5, 12.9, 8.7, 4.8, 4.2, 8.1, 6.5, 5.8, 6.7, 2.4, 11.1, 7.1, 8.8, 7.2

• Proszę znaleźć średnią średnicę i 90% przedział ufności dla średniej.
• Proszę znaleźć medianę i 90% przedział ufności dla mediany.
• Obydwa wyniki zilustrować przy pomocy histogramów.

Odp:

średnia: 7,50 i przedzial ufnosci dla średniej: [6,61 8,41];

mediana 7,15 i przedzial ufnosci dla mediany: [6,50 8,50].

Przedział ufności dla proporcji

W badaniach nad cholesterolem u ludzi stwierdzono, że w grupie 135 badanych z wysokim poziomem cholesterolu 10 osób przeszło zawał serca.

Pytanie: Jaki jest 95% przedział ufności dla proporcji 10/135?

• Proszę wykorzystać metodę bootstrapu.
• Proszę wykorzystać rozkład dwumianowy.

Odp: proporcja 10/135 = 0,07 i jej 95% przedział ufności [0,03 0,12].

Bezrobotni

W próbce 200 osób 7 procent jest bezrobotnych. Określić 95% przedział ufności dla odsetka bezrobotnych w populacji.

Odp: Średnia 7% i jej 95% przedział ufności [3,50 10,50].

Żywotność baterii

W próbce 20 testowanych baterii stwierdzono średni czas życia 28,85 miesiąca. Określić 95% przedział ufności dla średniej. Wartości dla badanej próbki były następujące:
30, 32, 31, 28, 31, 29, 29, 24, 30, 31, 28, 28, 32, 31, 24, 23, 31, 27, 27, 31 miesięcy
Obejrzeć rozkład przy pomocy histfit i zbadać jaki wpływ na przedział ufności ma przyjęcie założenia o normalności rozkładu czasów życia.

Odp: Średnia 28,85 i jej 95% przedział ufności [27,65 29,95].

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)
	# Tu 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 array z centrami binów korzystając z granic słupków
	# zwróconych przez py.hist w macierzy bins
 
	y = st.norm.pdf( bincenters, loc = np.mean(x), scale = np.std(x))
	# obliczamy wartości w normalnym rozkładzie gęstości prawdopodobieństwa
	# o średniej mi i wariancji sigma**2 dla wartości bincenters
 
	l = py.plot(bincenters, y, 'r--', linewidth=1)
	# do histogramu dorysowujemy linię

Pomiary

Mamy 10 pomiarów pewnej wielkości:
0.02, 0.026, 0.023, 0.017, 0.022, 0.019, 0.018, 0.018, 0.017, 0.022
Proszę znaleźć średnią i 95% przedział ufności.

Odp:Średnia 0,020 i jej 95% przedział ufności [0,0185 0,0220].

Czy pomiarów jest wystarczająco dużo aby sensownie wyznaczyć średnią i przedział ufności?

Wskazówka: Obliczyć średnie dla 1 000 000 powtórzeń i obejrzeć histogramy dla 10, 20, 30, 100, i 200 przedziałów.