poprawny sposób uzyskać przedział ufności z scipy

Question

mam 1-wymiarową tablicę danych:

a = np.array([1,2,3,4,4,4,5,5,5,5,4,4,4,6,7,8]) 

dla którego chcę, aby uzyskać 68% przedział ufności (IE: 1 sigma).

Pierwszy komentarz w this answer stwierdza, że ​​cel ten można osiągnąć stosując scipy.stats.norm.interval z funkcji scipy.stats.norm poprzez:

from scipy import stats 
import numpy as np 
mean, sigma = np.mean(a), np.std(a) 

conf_int = stats.norm.interval(0.68, loc=mean, scale=sigma) 

Ale komentarz w this post stwierdza, że ​​rzeczywisty poprawny sposób uzyskiwania przedział ufności jest :

conf_int = stats.norm.interval(0.68, loc=mean, scale=sigma/np.sqrt(len(a))) 

czyli współczynnik 1/np.sqrt(len(a)) stosowany jest w sigmie.

Pytanie brzmi: która wersja jest poprawna?

Answer 1

68% przedział ufności dla jeden zwrócić z rozkładu normalnego o średnią i odchyleniem standardowym sigma

stats.norm.interval(0.68, loc=mu, scale=sigma) 

68% przedział ufności dla średnią N zwraca z rozkład normalny ze średnią i odchyleniem std Sigma

stats.norm.interval(0.68, loc=mu, scale=sigma/sqrt(N)) 

Intuitiv Eee, te formuły mają sens, ponieważ jeśli trzymasz słoik z galaretowatymi ziarnami i poprosisz dużą liczbę ludzi, aby odgadł liczbę galaretek, każdy z nich może być odstąpiony o wiele - to samo odchylenie standardowe sigma - ale średnia z domysłów wykona niezwykle dobrą robotę szacowania faktycznej liczby, co znajduje odzwierciedlenie w odchyleniu standardowym średniego obkurczania o współczynnik 1/sqrt(N).

---

Jeśli pojedyncze losowanie ma wariancję sigma**2, następnie przez Bienaymé formula suma Nnieskorelowane rysuje ma wariancję N*sigma**2.

Średnia jest równa sumie podzielonej przez N. Przy pomnożeniu zmiennej losowej (takiej jak suma) przez stałą, wariancja jest mnożona przez stałą kwadratową. To jest

Var(cX) = c**2 * Var(X) 

Więc wariancja średniej równa

(variance of the sum)/N**2 = N * sigma**2/N**2 = sigma**2/N 

a więc odchylenie standardowe od średniej (który to pierwiastek kwadratowy z wariancji) równa

sigma/sqrt(N). 

ten jest źródłem sqrt(N) w mianowniku.

---

Oto przykład kodu, na podstawie Toma kodu, który demonstruje roszczeń powyżej:

import numpy as np 
from scipy import stats 

N = 10000 
a = np.random.normal(0, 1, N) 
mean, sigma = a.mean(), a.std(ddof=1) 
conf_int_a = stats.norm.interval(0.68, loc=mean, scale=sigma) 

print('{:0.2%} of the single draws are in conf_int_a' 
     .format(((a >= conf_int_a[0]) & (a < conf_int_a[1])).sum()/float(N))) 

M = 1000 
b = np.random.normal(0, 1, (N, M)).mean(axis=1) 
conf_int_b = stats.norm.interval(0.68, loc=0, scale=1/np.sqrt(M)) 
print('{:0.2%} of the means are in conf_int_b' 
     .format(((b >= conf_int_b[0]) & (b < conf_int_b[1])).sum()/float(N))) 

wydruków

68.03% of the single draws are in conf_int_a 
67.78% of the means are in conf_int_b 

Pamiętaj, że jeśli zdefiniujesz conf_int_b z szacowaną mean i sigma na podstawie próbki a, średnia nie może spaść w conf_int_b z wybraną częstotliwością .

---

Jeśli wziąć próbkę z dystrybucją i obliczyć próbkę średnią i odchylenie std,

mean, sigma = a.mean(), a.std() 

uważać, aby pamiętać, że nie ma gwarancji, że będą one równa populacja średnia i odchylenie standardowe oraz że jesteśmy przy założeniu, że populacja jest normalnie dystrybuowana - to nie są automatyczne dodatki!

Jeśli wziąć próbkę i chcą szacunków populacja średnia i odchylenie standardowe , należy użyć

mean, sigma = a.mean(), a.std(ddof=1) 

ponieważ ta wartość sigma jest unbiased estimator dla odchylenia standardowego populacji.

Answer 2

Przetestowałem Twoje metody za pomocą tablicy o znanym przedziale ufności. numpy.random.normal (mu, std, size) zwraca tablicę wyśrodkowaną na mu ze standardowym odchyleniem std (w the docs, jest to zdefiniowane jako Standard deviation (spread or “width”) of the distribution.).

from scipy import stats 
import numpy as np 
from numpy import random 
a = random.normal(0,1,10000) 
mean, sigma = np.mean(a), np.std(a) 
conf_int_a = stats.norm.interval(0.68, loc=mean, scale=sigma) 
conf_int_b = stats.norm.interval(0.68, loc=mean, scale=sigma/np.sqrt(len(a))) 


conf_int_a 
(-1.0011149125527312, 1.0059797764202412) 
conf_int_b 
(-0.0076030415111100983, 0.012467905378619625) 

Ponieważ wartość Sigma należy -1 do 1, przy czym sposób / np.sqrt(len(a)) jest nieprawidłowy.

Edit

Ponieważ nie mam reputacji wypowiedzenia powyżej będę wyjaśnić jak krawaty Ta odpowiedź język dokładnej odpowiedzi unutbu użytkownika. Jeśli wypełnisz losową tablicę z rozkładem normalnym, 68% całości będzie mieściło się w przedziale 1- σ średniej. W powyższym przypadku, jeśli sprawdzeniu, że widać

b = a[np.where((a>-1)&(a <1))] 
len(a) 
> 6781 

lub 68% populacji mieści się w 1 σ. Cóż, około 68%. Gdy użyjesz większej i większej macierzy, zbliżysz się do 68% (w próbie 10, 9 było pomiędzy -1 a 1). Dzieje się tak, ponieważ 1- σ jest nieodłączną dystrybucją danych i im więcej danych masz, tym lepiej możesz go rozwiązać.

Zasadniczo moja interpretacja pytania brzmiała: Jeśli mam próbkę danych, które chcę wykorzystać do opisania rozkładu, z którego zostały narysowane, jaka jest metoda znalezienia odchylenia standardowego tych danych? podczas gdy interpretacja unutbu wydaje się być bardziej Jaki jest przedział, w którym mogę umieścić średnią z 68% pewnością?. Co oznaczałoby, dla żelki, odpowiedziałem Jak się domyślam i unutbu odpowiedział Co ich domysły mówią nam o żelki.

Answer 3

Właśnie sprawdziłem, w jaki sposób R i GraphPad obliczają przedziały ufności i zwiększają interwał w przypadku małej liczności próbki (n). Np. Ponad 6-krotnie dla n = 2 w porównaniu do dużego n. Ten kod (na podstawie shasan na answer) pasuje do ich przedziałów ufności:

import numpy as np, scipy.stats as st 

# returns confidence interval of mean 
def confIntMean(a, conf=0.95): 
    mean, sem, m = np.mean(a), st.sem(a), st.t.ppf((1+conf)/2., len(a)-1) 
    return mean - m*sem, mean + m*sem 

Dla R, sprawdziłem przed t.test (a). Strona GraphPad confidence interval of a mean zawiera informacje o "poziomie użytkownika" dotyczące zależności wielkości próby.

Tutaj wyjście na przykład Gabriela:

In [2]: a = np.array([1,2,3,4,4,4,5,5,5,5,4,4,4,6,7,8]) 

In [3]: confIntMean(a, 0.68) 
Out[3]: (3.9974214366806184, 4.877578563319382) 

In [4]: st.norm.interval(0.68, loc=np.mean(a), scale=st.sem(a)) 
Out[4]: (4.0120010966037407, 4.8629989033962593) 

Należy zauważyć, że różnica między okresami confIntMean()st.norm.interval() i jest stosunkowo niewielka tutaj; len (a) == 16 nie jest za mały.