Guia 9, Ejercicio 14¶
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
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)
errdict = dict(ecolor='r', elinewidth=3, capsize=5, capthick=3, alpha=0.75)
formato = dict(color='r', alpha=0.5, error_kw=errdict)
Los siguientes son los datos del ejercicio 7 de la guia 5, los cuales se ajusta por cuadrados minimos a una recta. Se utiliza la función polyfit sin considerar el error en la variable y ya que es igual para todos y por lo tanto no cambia el resultado.
X = [2.00, 2.10, 2.20, 2.30, 2.40, 2.50, 2.60, 2.70, 2.80, 2.90, 3.00] #Datos en x
Y = [2.78, 3.29, 3.29, 3.33, 3.23, 3.69, 3.46, 3.87, 3.62, 3.40, 3.99] #Datos en y
Sy = 0.3 # Error en y
def recta(x, p):
x = np.asarray(x)
return p[1] + p[0] * x
p = np.polyfit(X, Y, 1) # Ajusto por una función lineal
print('Pendiente: {:.3f}, Ordenada: {:.3f}'.format(p[0], p[1]))
Se considera esta recta obtenida como la teórica y se generan a partir de esta un nuevo conjunto $\{\hat{y}_i\}$ de $11$ datos, donde cada $\hat{y}_i$ tiene distribución normal $N(y_i = f(x_i), \, \sigma)$.
A partir este conjunto se realizan 2 calculos de $\chi^2$:
El primero es de los nuevos puntos respecto de la recta original, es decir:
$$ \chi^2 = \sum_{i=1}^{11} \left(\frac{\hat{y}_i - y_i}{\sigma}\right)^2 = \sum_{i=1}^{11} \chi^2_{1} = \chi^2_{11} $$Este estadístico tiene una distribución $\chi^2$ de $11$ grados de libertad.
Por otro lado, a los nuevos $\hat{y}_i$ se les ajusta una nueva recta de la forma $y = ax + b$ y se obtienen unos nuevos estimadores $\hat{a}$ y $\hat{b}$. A partir de estos se calcula entonces:
$$ \chi^2 = \sum_{i=1}^{11} \left(\frac{\hat{y}_i - (\hat{a}x_i + \hat{b})}{\sigma}\right)^2 = \chi^2_{9} $$Se puede demostrar que al estar estimando $p = 2$ parametros con $n = 11$ datos, este ultimo estadístico es una variable $\chi^2$ de $\nu = n-p = 9$ grados de libertad. Las hipotesis para que esto suceda es que los errores en las variables $y_i$ sean gaussianos, y la función a estimar sea lineal en los parametros.
Se repite este proceso $1000$ veces y se generan 2 histogramas de los estadisticos $\chi^2$. Para el primero se superpone con su distribución correspondiente, mientras que el segundo histograma se muestra con las distribuciones $\chi^2$ de $9$ y $11$ grados de libertad. Es decir, la correcta y la cual no se tuvo en cuenta que se ajustaron parametros, respectivamente.
A simple vista, para el primer caso, la $\chi^2_{11}$ es representativa de los datos, mientras que en el segundo caso no lo es. Se realiza un test de Kolmogorov cuya hipotesis nula es que los datos provengan de una distribución $\chi^2_\nu$, y se observa que los p-values para los $\nu$ incorrectos son insignificantes.
def chi2(f, p, x, y, sy):
return np.sum(((y-f(x, p))/sy)**2)
chi2_teo = []
chi2_fit = []
for i in range(1000):
Yr = st.norm(loc=recta(X, p), scale=Sy).rvs(len(X)) # Genero datos con distribución normal
pr = np.polyfit(X, Yr, 1) # Ajusto la recta
chi2_teo.append(chi2(recta, p, X, Yr, Sy))
chi2_fit.append(chi2(recta, pr, X, Yr, Sy))
x = np.linspace(0, 35, 100)
hist(chi2_teo, normed=True, bins=20, **formato)
plt.plot(x, st.chi2.pdf(x, len(X)), label='$\chi^2_{(11)}$')
plt.xlabel('$\chi^2$')
plt.legend()
print('Kolmogorov-Smirnov test\n')
print('{} | {}'.format('ν'.rjust(2), 'p-value'))
print('------------')
for i in range(9, 14):
p = st.kstest(chi2_teo, 'chi2', [i])[1]
print('{} | {:.3g}'.format(str(i).rjust(2), p))
hist(chi2_fit, normed=True, bins=20, **formato)
plt.plot(x, st.chi2.pdf(x, len(X)), label='$\chi^2_{(11)}$')
plt.plot(x, st.chi2.pdf(x, len(X)-2), label='$\chi^2_{(9)}$')
plt.xlabel('$\chi^2$')
plt.legend()
print('Kolmogorov-Smirnov test\n')
print('{} | {}'.format('ν'.rjust(2), 'p-value'))
print('------------')
for i in range(7, 12):
p = st.kstest(chi2_fit, 'chi2', [i])[1]
print('{} | {:.3g}'.format(str(i).rjust(2), p))
