import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = (10, 5)
plt.rcParams['lines.linewidth'] = 3
plt.rcParams['ytick.labelsize'] = 12
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['axes.labelsize'] = 20
plt.rcParams['legend.fontsize'] = 20
plt.rcParams['legend.numpoints'] = 1
Guia 9, Ejercicio 13¶
a)¶
Viendo a ojo, la figura con menos cantidad de entradas es la 1, seguida de la 3 y finalmente la 2.
b)¶
Se pide realizar un histograma de los primeros $3000$ datos, y realizarle los test de $\chi^2$, Kolmogorov y Cramer-Von Mises para ver si provienen de una distribución $N(0, 2.5)$.
El test de $\chi^2$, a diferencia de los otros, depende del histograma. En este caso se obtuvo un p-value de $0.149$. Es decir, si la hipotesis es correcta, hay un $14.9\%$ de probabilidad de haber obtenido lo que se obtuvo o un valor mas extremo.
Para los test de Kolmogorov y Cramer-Von Mises se obtuvo un p-value (p) de $0.004$ y $0.033$, respectivamente, lo que nos permite rechazar la hipotesis nula con un confidence level (CL) de $96.7\%$.
def hist(data, bins=10, normed=False, plot=True, ax=None, **kwargs):
'''Grafica histogramas con barras de error'''
y, bins = np.histogram(data, bins=bins) #Divido los datos en bins
binwidth = np.diff(bins)
yerr = hist_err(y, binwidth)
if normed:
norm = np.sum(y) * binwidth
y = y / norm
yerr = yerr / norm
if plot:
if ax is None:
f, ax = plt.subplots(1,1)
ax.bar(bins[:-1], y, width=binwidth, yerr=yerr, **kwargs)
return y, yerr, bins
def hist_err(y, binwidth):
'''Calcula el error (binomial) para un histograma'''
N = np.sum(y)
p = y * binwidth / np.sum(y * binwidth)
return np.sqrt(p * (1 - p) * N)
# Cargo los datos
data = np.loadtxt('datos-G9E13.dat')
# Grafico el histograma y la N(0, 2.5)
bins = np.linspace(-7, 7, 57)
y, yerr, bins = hist(data[:3000], bins=bins, normed=True, plot=False)
x = (bins[1:]+bins[:-1])/2
f = st.norm.pdf(x, scale=2.5)
plt.errorbar(x, y, yerr=yerr, fmt='o', label='Datos')
plt.plot(x, f, label='N(0, 2.5)')
plt.legend()
# Test de Chi cuadrado
def chi2test(y, yerr, f):
chi2 = np.sum( ( (y-f)/yerr )**2 )
p = st.chi2.sf(chi2, y.size)
return chi2, p
chi2, p = chi2test(y, yerr, f)
print('Test de χ²')
print('χ² = {:.3f}, p = {:.3f}'.format(chi2, p))
# Test de Kolmogorov-Smirnov
ks, p = st.kstest(data[:3000], 'norm', [0, 2.5])
print('Test de Kolmogorov-Smirnov')
print('K-S = {:.3f}, p = {:.3f}'.format(ks, p))
# Test de Cramer-Von Mises
# Se necesita tener instalado R y el paquete goftest.
# goftest_1.0-3
from rpy2.robjects.packages import importr
from rpy2.robjects import numpy2ri
numpy2ri.activate()
goftest = importr('goftest')
cvm = goftest.cvm_test(data[:3000], 'pnorm', mean=0, sd=2.5)
print('Test de Cramer-Von Mises')
print('ω² = {:.3f}, p = {:.3f}'.format(cvm[0][0], cvm[1][0]))
c)¶
Finalmente, repetimos los 3 tests para distintas cantidades de datos.
Para el test de $\chi^2$ (curva azul), el p-value da practicamente $1$ cuando se considera hasta los primeros $60$ datos. Luego decrece y oscila entre $0$ y $1$, para finalmente tender a $0$ a partir de los $4000$ datos.
Por su parte, los tests de Kolmogorov-Smirnov y Cramer-Von Mises también oscilan entre $0$ y $1$ para valores anteriores a los $4000$ datos, para luego tender a $0$. Pero a diferencia del test de $\chi^2$, estos empiezan con un valor cercano al $20\%$ para cantidades bajas de datos, y llegan a un maximo de entre $75\%$ y $80\%$ para cantidades de datos entre $100$ y $500$. Es decir, a diferencia del test de $\chi^2$, el p-value nunca llega a valores cercanos a $1$.
Esto se debe a que no tiene mucho sentido hacer un histograma cuya cantidad de bins sea del orden de magnitud de la cantidad de datos.
ndata = np.rint(np.logspace(1, 5, 1000)).astype(int)
ndata = np.unique(ndata)
bins = np.linspace(-7, 7, 57)
x = (bins[1:]+bins[:-1])/2
f = st.norm.pdf(x, scale=2.5)
chi2_p = []
for i in ndata:
y, yerr, bins = hist(data[:i], bins=bins, normed=True, plot=False)
yerr[yerr==0] = np.min(yerr[yerr>0])
chi2_p.append(chi2test(y, yerr, f)[1])
ks_p = [st.kstest(data[:i], 'norm', [0, 2.5])[1] for i in ndata]
cvm_p = [goftest.cvm_test(data[:i], 'pnorm', mean=0, sd=2.5)[1][0] for i in ndata]
plt.plot(ndata, chi2_p, label='$\chi^2$')
plt.plot(ndata, ks_p, label='Kolmogorov')
plt.plot(ndata, cvm_p, label='Cramer-Von Mises')
plt.xscale('log')
plt.legend()
plt.xlabel('Cantidad de datos')
plt.ylabel('p-value')
plt.ylim((-0.02, 1.02))
plt.grid()
