python - fit - Manera correcta de obtener intervalo de confianza con scipy

scipy stats fit

Question

Tengo una matriz de datos unidimensional:

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

para lo cual quiero obtener el intervalo de confianza del 68% (es decir, el sigma 1 ).

El primer comentario en esta respuesta indica que esto se puede lograr usando scipy.stats.norm.interval desde la función scipy.stats.norm , a través de:

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)

Pero un comentario en esta publicación indica que la forma correcta real de obtener el intervalo de confianza es:

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

es decir, sigma se divide por la raíz cuadrada del tamaño de la muestra: np.sqrt(len(a)) .

La pregunta es: ¿qué versión es la correcta?

Answer 1

Acabo de comprobar cómo R y GraphPad calculan los intervalos de confianza y aumentan el intervalo en caso de un tamaño de muestra pequeño (n). Por ejemplo, más de 6 veces para n = 2 en comparación con una n grande. Este código (basado en la answer de shasan) coincide con sus intervalos de confianza:

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

Para R, he comprobado contra t.test (a). El intervalo de confianza de GraphPad de una página media tiene información de "nivel de usuario" sobre la dependencia del tamaño de la muestra.

Aquí la salida para el ejemplo de Gabriel:

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)

Tenga en cuenta que la diferencia entre los confIntMean() y st.norm.interval() es relativamente pequeña aquí; len (a) == 16 no es demasiado pequeño.

Answer 2

El intervalo de confianza del 68% para un solo sorteo de una distribución normal con mu media y desviación estándar sigma es

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

El intervalo de confianza del 68% para la media de N se extrae de una distribución normal con la media mu y la desviación estándar sigma es

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

De manera intuitiva, estas fórmulas tienen sentido, ya que si levanta un frasco de caramelos y le pide a un gran número de personas que adivinen el número de caramelos, cada individuo puede perder mucho - la misma desviación estándar sigma - pero el promedio de las conjeturas hará un trabajo notablemente bueno al estimar el número real y esto se refleja en la desviación estándar de la reducción de la media por un factor de 1/sqrt(N) .

Si un solo sorteo tiene varianza sigma**2 , entonces, según la fórmula de Bienaymé , la suma de N sorteos no correlacionados tiene varianza N*sigma**2 .

La media es igual a la suma dividida por N. Cuando multiplicas una variable aleatoria (como la suma) por una constante, la varianza se multiplica por la constante al cuadrado. Es decir

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

Así que la varianza de la media es igual a

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

y así la desviación estándar de la media (que es la raíz cuadrada de la varianza) es igual a

sigma/sqrt(N).

Este es el origen del sqrt(N) en el denominador.

Aquí hay un código de ejemplo, basado en el código de Tom, que demuestra las afirmaciones hechas anteriormente:

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)))

huellas dactilares

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

Tenga en cuenta que si define conf_int_b con las estimaciones de mean y sigma basadas en la muestra a , es posible que la media no caiga en conf_int_b con la frecuencia deseada.

Si toma una muestra de una distribución y calcula la media de la muestra y la desviación estándar,

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

tenga cuidado de notar que no hay garantía de que estos sean iguales a la media de la población y la desviación estándar y que estamos asumiendo que la población está normalmente distribuida, ¡eso no se da automáticamente!

Si toma una muestra y desea estimar la media de la población y la desviación estándar, debe usar

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

ya que este valor para sigma es el estimador imparcial para la desviación estándar de la población.

Answer 3

Probé sus métodos utilizando una matriz con un intervalo de confianza conocido. numpy.random.normal (mu, std, size) devuelve una matriz centrada en mu con una desviación estándar de std (en los documentos , esto se define como la 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)

Como el valor sigma debe ser -1 a 1, el / np.sqrt(len(a)) parece ser incorrecto.

Editar

Como no tengo la reputación de comentar más arriba, aclararé cómo esta respuesta se relaciona con la respuesta completa de unutbu. Si llena una matriz aleatoria con una distribución normal, el 68% del total estará dentro de 1-σ de la media. En el caso anterior, si verifica que ve

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

o el 68% de la población cae dentro de 1σ. Bueno, alrededor del 68%. A medida que usa una matriz más y más grande, se acercará al 68% (en una prueba de 10, 9 estaban entre -1 y 1). Esto se debe a que 1-σ es la distribución inherente de los datos, y cuantos más datos tenga, mejor podrá resolverlos.

Básicamente, mi interpretación de su pregunta fue: Si tengo una muestra de datos que quiero usar para describir la distribución de la que se obtuvieron, ¿cuál es el método para encontrar la desviación estándar de esos datos? mientras que la interpretación de Unutbu parece ser más ¿Cuál es el intervalo en el que puedo colocar la media con un 68% de confianza? . Lo que significaría, para jelly beans, respondí Cómo están adivinando y unutbu contestó ¿Qué nos dicen sus suposiciones acerca de jelly beans?