rvs - Implementando una prueba de Kolmogorov Smirnov en python scipy

rvs python

Question

Tengo un conjunto de datos en N números que quiero probar para la normalidad. Sé que scipy.stats tiene una función kstest pero no hay ejemplos sobre cómo usarla y cómo interpretar los resultados. ¿Hay alguien aquí que esté familiarizado con eso que pueda darme algún consejo?

Según la documentación, el uso de kstest devuelve dos números, el estadístico de prueba KS D y el valor p. Si el valor de p es mayor que el nivel de significación (por ejemplo, 5%), no podemos rechazar la hipótesis de que los datos provienen de la distribución dada.

Cuando hago una ejecución de prueba extrayendo 10000 muestras de una distribución normal y probando la gaussianidad:

import numpy as np
from scipy.stats import kstest

mu,sigma = 0.07, 0.89
kstest(np.random.normal(mu,sigma,10000),''norm'')

Me sale el siguiente resultado:

(0.04957880905196102, 8.9249710700788814e-22)

El valor p es inferior al 5%, lo que significa que podemos rechazar la hipótesis de que los datos se distribuyen normalmente. ¡Pero las muestras fueron extraídas de una distribución normal!

¿Alguien puede entenderme y explicarme la discrepancia aquí?

(¿Las pruebas de normalidad suponen que mu = 0 y sigma = 1? Si es así, ¿cómo puedo probar que mis datos están distribuidos de manera gaussiana pero con una mu y sigma diferentes?)

Answer 1

Como complemento a la respuesta de @unutbu, también podría proporcionar los parámetros de distribución para la distribución de prueba en kstest. Supongamos que tenemos algunas muestras de una variable (y las llamamos datax), y quisiéramos comprobar si esas muestras no podrían provenir de un lognormal, un uniforme o una normal. Tenga en cuenta que para las estadísticas de scipy, la forma en que se toman los parámetros de entrada para cada distribución varía un poco. Ahora, gracias a "args" (tupla o secuencia) en kstest, es posible proporcionar los argumentos para la distribución scipy.stats contra la que desea probar.

:) También agregué la opción de usar una prueba de dos muestras, en caso de que quisieras hacerlo de cualquier manera:

import numpy as np
from math import sqrt
from scipy.stats import kstest, ks_2samp, lognorm
import scipy.stats

def KSSeveralDists(data,dists_and_args,samplesFromDists=100,twosampleKS=True):
    returnable={}
    for dist in dists_and_args:
        try:
            if twosampleKS:
                try:
                    loc=dists_and_args[dist][0]
                    scale=dists_and_args[dist][1]
                    expression=''scipy.stats.''+dist+''.rvs(loc=loc,scale=scale,size=samplesFromDists)''
                    sampledDist=eval(expression)
                except:
                    sc=dists_and_args[dist][0]
                    loc=dists_and_args[dist][1]
                    scale=dists_and_args[dist][2]
                    expression=''scipy.stats.''+dist+''.rvs(sc,loc=loc,scale=scale,size=samplesFromDists)''
                    sampledDist=eval(expression)
                D,p=ks_2samp(data,sampledDist)
            else:
                D,p=kstest(data,dist,N=samplesFromDists,args=dists_and_args[dist])
        except:
            continue
        returnable[dist]={''KS'':D,''p-value'':p}
    return returnable

a=lambda m,std: m-std*sqrt(12.)/2.
b=lambda m,std: m+std*sqrt(12.)/2.
sz=2000

sc=0.5 #shape 
datax=lognorm.rvs(sc,loc=0.,scale=1.,size=sz)
normalargs=(datax.mean(),datax.std())

#suppose these are the parameters you wanted to pass for each distribution
dists_and_args={''norm'':normalargs,
               ''uniform'':(a(*normalargs),b(*normalargs)),
               ''lognorm'':[0.5,0.,1.]
              }
print "two sample KS:"
print KSSeveralDists(datax,dists_and_args,samplesFromDists=sz,twosampleKS=True)
print "one sample KS:"
print KSSeveralDists(datax,dists_and_args,samplesFromDists=sz,twosampleKS=False)

Lo que da como salida algo como:

dos muestras KS: {''lognorm'': {''KS'': 0.02349999999999999965, ''p-value'': 0.63384188886455217}, ''norm'': {''KS'': 0.1060000000000000004, ''p-value'': 2.918766666723155e-10}, ''uniforme'' '': {'' KS '': 0.15300000000000002,'' p-value '': 6.443660021191129e-21}}

una muestra KS: {''lognorm'': {''KS'': 0.01763415915126032, ''p-value'': 0.56275820961065193}, ''norm'': {''KS'': 0.10792612430093562, ''p-value'': 0.0}, ''uniform'': { ''KS'': 0.14910036159697559, ''p-value'': 0.0}}

Nota: Para la distribución uniforme de scipy.stats, a y b se toman como a = loc y b = loc + escala (consulte la documentation ).

Answer 2

Sus datos fueron generados con mu = 0.07 y sigma = 0.89. Está probando estos datos contra una distribución normal con media 0 y desviación estándar de 1.

La hipótesis nula ( H0 ) es que la distribución de la cual sus datos son una muestra es igual a la distribución normal estándar con media 0, desviación estándar 1.

El pequeño valor p indica que se esperaría un estadístico de prueba tan grande como D con el valor p de probabilidad.

En otras palabras, (con un valor de p ~ 8.9e-22) es altamente improbable que H0 sea ​​cierto.

Eso es razonable, ya que los medios y las desviaciones estándar no coinciden.

Compara tu resultado con:

In [22]: import numpy as np
In [23]: import scipy.stats as stats
In [24]: stats.kstest(np.random.normal(0,1,10000),''norm'')
Out[24]: (0.007038739782416259, 0.70477679457831155)

Para probar que sus datos son gaussianos, puede cambiarlos y volver a escalarlos para que sea normal con media 0 y desviación estándar 1:

data=np.random.normal(mu,sigma,10000)
normed_data=(data-mu)/sigma
print(stats.kstest(normed_data,''norm''))
# (0.0085805670733036798, 0.45316245879609179)

Advertencia: ( muchas gracias a user333700 (también conocido como desarrollador scipy Josef Perktold )) Si no conoce mu y sigma , estimar los parámetros hace que el valor p sea inválido:

import numpy as np
import scipy.stats as stats

mu = 0.3
sigma = 5

num_tests = 10**5
num_rejects = 0
alpha = 0.05
for i in xrange(num_tests):
    data = np.random.normal(mu, sigma, 10000)
    # normed_data = (data - mu) / sigma    # this is okay
    # 4915/100000 = 0.05 rejects at rejection level 0.05 (as expected)
    normed_data = (data - data.mean()) / data.std()    # this is NOT okay
    # 20/100000 = 0.00 rejects at rejection level 0.05 (not expected)
    D, pval = stats.kstest(normed_data, ''norm'')
    if pval < alpha:
        num_rejects += 1
ratio = float(num_rejects) / num_tests
print(''{}/{} = {:.2f} rejects at rejection level {}''.format(
    num_rejects, num_tests, ratio, alpha))     

huellas dactilares

20/100000 = 0.00 rejects at rejection level 0.05 (not expected)

que muestra que stats.kstest no puede rechazar el número esperado de hipótesis nulas si la muestra se normaliza utilizando la media y la desviación estándar de la muestra

normed_data = (data - data.mean()) / data.std()    # this is NOT okay

Answer 3

También puede considerar el uso de la prueba de Shapiro-Wilk, que "prueba la hipótesis nula de que los datos se extrajeron de una distribución normal". También se implementa en scipy :

http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.shapiro.html

Tendrá que pasar sus datos directamente a la función.

import scipy

W, p = scipy.stats.shapiro(dataset)
print("Shapiro-Wilk test statistic, W:", W, "/n", "p-value:", p)

Lo que devuelve algo como:

 Shapiro-Wilk test statistic, W: 0.7761164903640747 
 p-value: 6.317247641091492e-37

Con p << 0.01 (o 0.05, si lo prefiere, no importa), tenemos buenas razones para rechazar la hipótesis nula de que estos datos se obtuvieron de una distribución normal.

Answer 4

Una actualización sobre la respuesta de unutbu:

Para las distribuciones que solo dependen de la ubicación y la escala pero que no tienen un parámetro de forma, las distribuciones de varias estadísticas de prueba de bondad de ajuste son independientes de la ubicación y los valores de escala. La distribución no es estándar, sin embargo, se puede tabular y usar con cualquier ubicación y escala de la distribución subyacente.

La prueba de Kolmogorov-Smirnov para la distribución normal con ubicación y escala estimadas también se llama prueba de Lilliefors .

Ahora está disponible en statsmodels, con valores p aproximados para el rango de decisión relevante.

>>> import numpy as np
>>> mu,sigma = 0.07, 0.89
>>> x = np.random.normal(mu, sigma, 10000)
>>> import statsmodels.api as sm
>>> sm.stats.lilliefors(x)
(0.0055267411213540951, 0.66190841161592895)

La mayoría de los estudios de Monte Carlo muestran que la prueba de Anderson-Darling es más poderosa que la prueba de Kolmogorov-Smirnov. Está disponible en scipy.stats con valores críticos y en modelos de estadísticas con valores p aproximados:

>>> sm.stats.normal_ad(x)
(0.23016468240712129, 0.80657628536145665)

Ninguna de las pruebas rechaza la hipótesis nula de que la muestra tiene una distribución normal. Mientras que el kstest en la pregunta rechaza la hipótesis nula de que la muestra tiene una distribución normal estándar .