Alumno: Mauro Silberberg

Libreta: 408/10

Aca abajo copio el codigo para graficar histogramas que usé en el trabajo anterior.

In [1]:
import numpy as np
import scipy.optimize as opt
import scipy.stats as st
import matplotlib.pyplot as plt
%matplotlib inline

def hist(data, bins=10, normed=False, ax=None, **kwargs):
    '''Grafica histogramas con barras de error'''
    y, bins = np.histogram(data, bins=bins) #Divido los datos en bins
    binwidth = (bins[1:]-bins[:-1])
    
    N = np.sum(y) * binwidth if normed else 1
    yerr = hist_err(y, binwidth=binwidth) / N #Calculo y normalizo el error
    y = y / N #Normalizo los datos
    
    if ax == None:
        f, ax = plt.subplots(1,1)
    ax.bar(bins[:-1], y, width=binwidth, yerr=yerr, **kwargs)
    return y, bins

def hist_err(data, binwidth=1):
    '''Calcula el error (binomial) para un histograma'''
    N = np.sum(data * binwidth)
    p = data * binwidth / N
    return np.sqrt(p * (1 - p) * N)

errdict = dict(ecolor='red', elinewidth=3, capsize=5, capthick=3, alpha=0.75)
formato = dict(color='r', alpha=0.5, error_kw=errdict)

Guia 5, ejercicio 6

Si tenemos $N$ pares de mediciones $(x_i, y_i)$ y los ajustamos por cuadrados mínimos a una función lineal de la forma, $ y = a_1 x + a_0 $, la matriz de covarianza de los parametros queda:

$$ Cov(a_0,a_1) = \begin{pmatrix} Var(a_0) & Cov(a_0,a_1) \\ Cov(a_0,a_1) & Var(a_1) \end{pmatrix} = \frac{\sigma^2}{\Delta} \begin{pmatrix} \sum x_i^2 & -\sum x_i \\ -\sum x_i & N \end{pmatrix} $$

donde $\Delta = N\sum x_i^2 - (\sum x_i)^2$, $\sigma$ es el error en $y_i$ (igual para todo $i$), y hay un pequeño abuso de notacion entre la matriz $Cov(a_0,a_1)$ y el elemento $Cov(a_0,a_1)$.

Nota: Abajo se encuentran los datos $X$ e $Y$ y se calculan los parametros que mejor ajustan a la recta a traves de la funcion de curve_fit, que esta pensada para ajustes no lineales pero admite la opcion absolute_sigma, a diferencia de numpy.polyfit, que devuelve una matriz de covarianza escalada (tomando los errores como relativos), la cual tiene un factor de escala $\sim 1.5$ para estos datos. (Ver http://mail.scipy.org/pipermail/scipy-user/2013-February/034199.html)

Tambien encontré lo mismo en Origin, donde al aumentar el error $\sigma$ en las mediciones por un factor $\alpha$ (es decir, $\sigma \mapsto \alpha \sigma$), la matriz de covarianza que devuelve la funcion linear fit se mantiene igual en lugar de aumentar como $\alpha^2$ como se esperaria (ya que la matriz de covarianza es proporcional a $\sigma^2$ si los errores en la variable dependiente son todos iguales).

In [2]:
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

def recta(x, *p):
    x = np.asarray(x)
    p = np.asarray(p[0]) if len(p) == 1 else np.asarray(p)
    return p[0] + p[1] * x

p, cov = opt.curve_fit(recta, X, Y, p0=[1, 1], sigma=Sy, absolute_sigma=True) #Ajuste

Error con covarianza

Si calculamos el error $\sigma_y$ del $y$ predicho, tenemos:

$$ \sigma_y^2 = Var(y) = x^2 \, Var(a_1) + Var(a_0) + 2 x \, Cov(a_1,a_0) $$

Este se minimiza en:

$$ x_{min} = - \frac{Cov(a_1,a_0)}{Var(a_1)} = \frac{\sum x_i}{N} = \bar{x} $$

que es el promedio de los datos $x$. En este lugar, la recta obtenida pasa por el promedio de los puntos $y_i$:

$$y(\bar{x}) = \bar{y} = \frac{\sum y_i}{N}$$

En el minimo, el error tiene un valor de:

$$ \sigma_y \, (x_{min}) = \sigma_y \, (\bar{x}) = \frac{\sigma}{\sqrt{N}} $$

que es el mismo resultado que se obtiene al considerar el desvio estandar de la media. O sea, el error para $\bar{y}$ es el error de la media.

Cuando nos alejamos de la region donde se tienen las mediciones, es decir, cuando extrapolamos, el error crece.

Debajo se grafican las mediciones, la recta ajustada y bandas de error para 1, 2 y 3 $\sigma$. Se observa como crece el error conforme se aleja de $\bar{x} = 2.5$.

In [3]:
def error_recta(x, cov):
    x = np.asarray(x)
    cov = np.asarray(cov)
    return np.sqrt(cov[0,0] + cov[1, 1] * x**2 + 2 * x * cov[0, 1])

x = np.linspace(1.5, 3.5, 100)
y = recta(x, p)
ys = error_recta(x, cov)

#Recta ajustada
plt.plot(x, y, color='g', lw=3)
#Banda de 1 sigma
plt.fill_between(x, y + ys, y - ys, color='g', alpha=0.5)
#Banda de 2 sigma
plt.fill_between(x, y + 2*ys, y + 1*ys, color='y', alpha=0.5)
plt.fill_between(x, y - 2*ys, y - 1*ys, color='y', alpha=0.5)
#Banda de 3 sigma
plt.fill_between(x, y + 3*ys, y + 2*ys, color='r', alpha=0.5)
plt.fill_between(x, y - 3*ys, y - 2*ys, color='r', alpha=0.5)
#Datos
plt.errorbar(X, Y, yerr=Sy, fmt='o')
#Leyenda
p1 = plt.Rectangle((0, 0), 1, 1, fc='g', alpha=0.7)
p2 = plt.Rectangle((0, 0), 1, 1, fc='y', alpha=0.7)
p3 = plt.Rectangle((0, 0), 1, 1, fc='r', alpha=0.7)
plt.legend([p1, p2, p3], ['1$\sigma$', '2$\sigma$', '3$\sigma$'], loc=4)
Out[3]:
<matplotlib.legend.Legend at 0x7f66461b6eb8>

Distribucion no unimodal en x

Hay que tener cuidado al decir que el error es minimo en la región donde se tienen valores. Si se tiene una distribucion de las mediciones que no sea unimodal, el valor medio de los $x_i$ puede quedar en una región donde no se tienen datos, y, por lo tanto, el error ser minimo en esta región, como se puede ver en el siguiente grafico.

In [4]:
X2 = np.concatenate((np.linspace(2, 3, 5), np.linspace(5, 6, 5)))
Y2 = st.norm(loc=recta(X2, p), scale=0.3).rvs(len(X2))

p2, cov2 = opt.curve_fit(recta, X2, Y2, p0=[1, 1], sigma=0.3, absolute_sigma=True) #Ajuste

x = np.linspace(-0, 8, 100)
y = recta(x, p2)
ys = error_recta(x, cov2)

plt.plot(x, y, 'g', lw=3) #Recta ajustada
plt.fill_between(x, y + ys, y - ys, color='g', alpha=0.5) #Banda de 1 sigma
#Banda de 2 sigma
plt.fill_between(x, y + 2*ys, y + 1*ys, color='y', alpha=0.5)
plt.fill_between(x, y - 2*ys, y - 1*ys, color='y', alpha=0.5)
#Banda de 3 sigma
plt.fill_between(x, y + 3*ys, y + 2*ys, color='r', alpha=0.5)
plt.fill_between(x, y - 3*ys, y - 2*ys, color='r', alpha=0.5)
#Puntos
plt.errorbar(X2, Y2, yerr=Sy, fmt='o')
#Leyenda
p1 = plt.Rectangle((0, 0), 1, 1, fc='g', alpha=0.7)
p2 = plt.Rectangle((0, 0), 1, 1, fc='y', alpha=0.7)
p3 = plt.Rectangle((0, 0), 1, 1, fc='r', alpha=0.7)
plt.legend([p1, p2, p3], ['1$\sigma$', '2$\sigma$', '3$\sigma$'], loc=4)
Out[4]:
<matplotlib.legend.Legend at 0x7f66461c4940>

Error sin covarianza

Si nos olvidamos la covarianza de los parametros al propagar el error, es decir:

$$ \sigma_y^2 = Var(y) = x^2 \, Var(a_1) + Var(a_0) $$

el error se minimiza en:

$$ x_{min} = 0 $$

y tiene un valor de:

$$ \sigma_y \, (x_{min}) = \frac{\sigma}{N \sigma_x} $$

donde $\sigma_x$ es la desviación estandar de los datos $x$.

El error es mínimo en el origen, independiente de donde se encuentren nuestros datos, lo cual no tiene sentido.

Debajo se grafican las mediciones, la recta ajustada y las bandas de error sin considerar la covarianza. Se observa que las bandas son minimas alrededor del origen.

In [5]:
def error_recta_sin_cov(x, cov):
    x = np.asarray(x)
    cov = np.asarray(cov)
    return np.sqrt(cov[0,0] + cov[1, 1] * x**2)

x = np.linspace(-1, 4, 100)
y = recta(x, p)
ys = error_recta_sin_cov(x, cov)

#Recta ajustada
plt.plot(x, y, 'g', lw=3)
#Banda de 1 sigma
plt.fill_between(x, y + ys, y - ys, color='g', alpha=0.5)
#Banda de 2 sigma
plt.fill_between(x, y + 2*ys, y + 1*ys, color='y', alpha=0.5)
plt.fill_between(x, y - 2*ys, y - 1*ys, color='y', alpha=0.5)
#Banda de 3 sigma
plt.fill_between(x, y + 3*ys, y + 2*ys, color='r', alpha=0.5)
plt.fill_between(x, y - 3*ys, y - 2*ys, color='r', alpha=0.5)
#Datos
plt.errorbar(X, Y, yerr=Sy, fmt='o')
#Leyenda
p1 = plt.Rectangle((0, 0), 1, 1, fc='g', alpha=0.7)
p2 = plt.Rectangle((0, 0), 1, 1, fc='y', alpha=0.7)
p3 = plt.Rectangle((0, 0), 1, 1, fc='r', alpha=0.7)
plt.legend([p1, p2, p3], ['1$\sigma$', '2$\sigma$', '3$\sigma$'], loc=2)
plt.grid(b=True)

Guia 5, ejercicio 7

Para los valores de $x_i$ de las mediciones, se generan nuevos valores $y_i$ con una distribución normal $N(a_0 + a_1 x_i, \sigma)$, para los parametros hallados previamente. Para estos valores $(x_i, y_i)$, ajusta una recta, se obtienen nuevos parametros $a_0$ y $a_1$ y se predice a partir de estos el valor de $y_a(x_a = 0.5)$.

Se repite este procedimiento 1000 veces y se realiza un histograma a partir de los valores de $y_a$, el cual se lo compara con una normal centrada en el valor real de $y_a$ con ancho $\sigma_y (x_a)$, obtenidos a partir de los parametros $a_0$, $a_1$ y su matriz de covarianza originales.

In [6]:
Xa = 0.5
Ya = []

for i in range(1000):
    Yr = st.norm(loc=recta(X, p), scale=0.3).rvs(len(X))
    pr, covr = opt.curve_fit(recta, X, Yr, p0=[1, 1], sigma=0.3, absolute_sigma=True)
    
    Ya.append(recta(Xa, pr))

Ya_real = recta(Xa, p)
Ya_sigma = error_recta(Xa, cov)

x = np.linspace(-4*Ya_sigma, 4*Ya_sigma, 100) + Ya_real
y = st.norm(loc=Ya_real, scale=Ya_sigma).pdf(x)

hist(Ya, bins=20, normed=True, **formato)
plt.plot(x, y, lw=3)
plt.xlabel('$y_a$', fontsize=20)
Out[6]:
<matplotlib.text.Text at 0x7f6645fc5080>