import numpy as np
from scipy import stats as st
from matplotlib import pyplot as plt
% matplotlib inline
Defino abajo un par de funciones que voy a usar más adelante.
Noten que uso la funcion kstest con la alternativa 'two-sided', que no significa que es un test a dos colas (o sí, según como interpreten el estadistico). Lo que hace es devolver el máximo del valor absoluto de la diferencia (la interpretacion a dos colas sería sí además lo devolviera con el signo).
Pueden ver el codigo fuente (está hecho en Python también, es cortito y entendible) yendo a la ayuda y tocando un link arriba a la derecha que los lleva a GitHub: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kstest.html
Para ver cual es la diferencia entre las distintas alternativas que tiene el test, pueden leer acá: http://blog.thegrandlocus.com/2014/09/mis-using-the-ks-test-for-p-hacking
Nota: fijense que calculo el p-value (es decir, la integral de la PDF) a partir de la distribucion acumulada (empirica). Se tiene un poco más de resolución y no hay arbitrariedades de como eligen los bins del histograma.
def ks(n, it, rvs, cdf):
'''Calcula el estadistico de Kolmogorov-Smirnov.'''
ks = np.empty(it)
for k in range(it):
ks[k] = st.kstest(rvs(size=n), cdf, alternative='two-sided').statistic
return ks
def ks_mod(n, it, rvs, cdf):
'''Calcula el estadistico de Kolmogorov-Smirnov modificado.'''
ks = np.empty(it)
for k in range(it):
x = rvs(size=n)
x = (x - x.mean()) / x.std(ddof=1)
ks[k] = st.kstest(x, cdf, alternative='two-sided').statistic
return ks
def ecdf(data):
'''Distribución cumulativa empírica.
Agrego 2 puntos al inicio y al final (se que va a ir entre 0 y 1) para poder graficar más fácil.
'''
x, y = np.empty(data.size+2), np.arange(data.size+2) / data.size
x[0] = 0
x[1:-1] = np.sort(data)
x[-1] = 1
y[-1] = 1
return x, y
def pvalue(t, ecdf):
'''Calcula el p-value (cola a derecha) interpolando en la distribución cumulativa empirica.'''
return 1 - np.interp(t, *ecdf)
def t_critico(pvalue, ecdf):
'''Calcula el valor crítico interpolando en la distribución cumulativa empirica.'''
return np.interp(1 - pvalue, *reversed(ecdf))
ks = ks(n=10, it=10000, rvs=st.norm.rvs, cdf=st.norm.cdf)
ksm = ks_mod(n=10, it=10000, rvs=st.norm.rvs, cdf=st.norm.cdf)
f, ax = plt.subplots(ncols=2, figsize=(12,4))
for k, label in ((ks, 'KS'), (ksm, 'KS modificado')):
tc = t_critico(0.05, ecdf(k))
print('El valor crítico para {label} es {tc:.2f}'.format(label=label, tc=tc))
hist, bins = np.histogram(k, bins=100, density=True)
bins = bins[1:]
ax[0].step(bins, hist, where='post')
cond = bins >= tc
ax[0].fill_between(bins[cond], hist[cond], step='post')
ax[1].plot(*ecdf(k), label=label, zorder=10)
ax[1].hlines(1-0.05, 0, tc, linestyle='--')
ax[1].vlines(tc, 0, 1-0.05, linestyle='--')
ax[1].scatter(tc, 1-0.05)
ax[0].set_title('PDF')
ax[1].set_title('CDF')
ax[1].legend()
Para el punto 2 habia 4 combinaciones posibles, donde algunos hicieron unas y otros otras: calcular el estadistico comun o el modificado y calcular el pvalue segun la distribucion del comun o el modificado.
Cuando calculan el pvalue usando la distribucion correcta, obtienen una distribución uniforme (por definicion, a mi me da perfectamente uniforme porque lo calculé con los mismos datos que use para generar la distribucion empirica del estadistico).
En cambio, si calculan el pvalue con la distribucion incorrecta, se obtienen otras cosas. Por ejemplo, el caso interesante es el estadistico modificado pero calculando el pvalue con la distribucion de KS normal (abajo a la izquierda). La distribucion se corre hacia el 1, y por lo tanto, esperarian aceptar más de lo normal la hipotesis nula. Al calcular el estadistico modificado, están comparando contra una distribucion que se "ajusta" mejor a los datos, porque le ponen el $\mu$ y $\sigma$ calculados a partir de estos.
f, ax = plt.subplots(ncols=2, nrows=2, figsize=(12,6), sharex=True)
hist_opt = dict(bins=100, normed=True, histtype='bar')
ax[0,0].hist(pvalue(ks, ecdf(ks)), **hist_opt)
ax[0,1].hist(pvalue(ks, ecdf(ksm)), **hist_opt)
ax[1,0].hist(pvalue(ksm, ecdf(ks)), **hist_opt)
ax[1,1].hist(pvalue(ksm, ecdf(ksm)), **hist_opt)
ax[0,0].set_title('KS statistic & KS test')
ax[0,1].set_title('KS statistic & Modified KS test')
ax[1,0].set_title('Modified KS statistic & KS test')
ax[1,1].set_title('Modified KS statistic & Modified KS test')
for a in ax[1]: a.set_xlabel('p-value')
Los datos de Julian¶
Con el test comun probablemente aceptarian la hipotesis de que son gaussianos, mientras que con el test modificado probablemente la rechazarian.
x = np.array([float(k) for k in '6.55 4.39 3.80 3.53 5.84 3.51 3.76 3.61 3.91 4.23'.split(' ')])
xm, xs = x.mean(), x.std(ddof=1)
tj_ks = st.kstest(x, st.norm(loc=xm, scale=xs).cdf).statistic
p_ks = pvalue(tj_ks, ecdf(ks))
p_ksm = pvalue(tj_ks, ecdf(ksm))
print('El valor medio es {mean:.3f} y la desviación estandar es {std:.3f}'.format(mean=xm, std=xs))
print('El estadistico es {:.3f}'.format(tj_ks))
print('El p-value para KS es {p_ks:.3f} y para KS modificado es {p_ksm:.3f}'.format(p_ks=p_ks, p_ksm=p_ksm))
El poder del test modificado¶
Uno espera que a medida que aumente la cantidad de datos, pueda estimar mejor las distribuciones y pueda diferenciar entre una gaussiana y una t-student de 15 grados de libertad. Para esto, calculamos el poder. Primero generamos las distribuciones (que dependen de la cantidad de datos) del estadistico modificado asumiendo que la distribucion es normal o una t de Student.
from tqdm import tqdm
N = [10, 100, 1000, 10000]
it = 10000
ks_tstud = np.empty((len(N), it))
ks_norm = np.empty((len(N), it))
for ix, n in enumerate(tqdm(N)):
ks_norm[ix] = ks_mod(n=n, it=it, rvs=st.norm.rvs, cdf=st.norm.cdf)
ks_tstud[ix] = ks_mod(n=n, it=it, rvs=st.t(df=15).rvs, cdf=st.norm.cdf)
Si graficamos la distribucion del KS modificado en funcion de la cantidad de datos, vemos que se desplaza hacia el 0. Al tener mas datos, estamos estimando mejor la distribucion (empirica) y difiere cada vez menos de la teorica.
hist_opt = dict(bins=100, histtype='step')
for x, n in zip(ks_norm, N):
plt.hist(x, label=str(n), **hist_opt)
plt.legend()
plt.title('KS modificado en función del numero de datos')
Si comparamos las distribuciones de los estadisticos KS modificados para la normal y para la t de Student, vemos que a medida que crece la cantidad de datos, hay cada vez menos overlap entre estos. Si no hubiese overlap, serian completamente distinguibles: si cae en una dada region la hipotesis correcta es una y si cae en otra region seria la otra. Es decir, el poder sería 1.
f, axes = plt.subplots(ncols=4, figsize=(16,4))
hist_opt = dict(bins=100, normed=True, histtype='step')
for k in range(4):
axes[k].hist(ks_norm[k], label='Normal', **hist_opt)
axes[k].hist(ks_tstud[k], label='t-Student(15)', **hist_opt)
axes[k].set_title('N={}'.format(N[k]))
axes[k].legend()
El poder depende del valor critico que elijamos, o equivalentemente, del la significancia. Para una significancia de 0.05, calculamos el valor critico correspondiente (para H0: provienen de una gaussiana) y con ese valor critico calculamos la integral (sobre ls distribucion correspondiente a H1: provienen de una t de Student) y obtenermos el poder.
for k in range(4):
tc = t_critico(0.05, ecdf(ks_norm[k])) # Calculo el t_critico
poder = pvalue(tc, ecdf(ks_tstud[k])) # Calculo el poder (la integral es en la misma zona que el p-value y reuso esa funcion)
print('Para tiras de {n:>5} datos, el valor crítico del test modificado es {tc:.4f} y el poder es {power:.4f}'.format(n=N[k], tc=tc, power=poder))
Podemos hacer esto mismo para muchos valores de significancia, y ver como varia el poder en funcion de la cantidad de datos.
p = np.linspace(0, 1, 1000)
for k, n in enumerate(N):
tc = t_critico(p, ecdf(ks_norm[k]))
poder = pvalue(tc, ecdf(ks_tstud[k]))
plt.plot(p, poder, label=n)
plt.legend()
plt.xlabel('Significancia $\\alpha$')
plt.ylabel('Poder $1-\\beta$')
