Episodio II

Willy Pregliasco y Mariano Gómez Berisso, Ago 2016

Haciendo un ajuste por una recta

1. Repaso

En la clase anterior vimos el caso particular de un ajuste realizado para el caso de ruido gaussiano y prior uniforme: tengo maxima ignorancia sobre los posibles valores de los parametros del ajuste.
Como vimos las cosas se simplifican mucho.

Si el ruido es gaussiano, resulta que los datos responden a una distribución: $$ y_i \sim \cal{N}(\mu=f(x_i,\theta)\ ,\ \sigma=\sigma_i) $$

y eso hace que (por la definición de la gaussiana): $$ P(y_i\ |x_i \theta) \propto \exp\left( -\frac{1}{2\sigma_i^2} \left(y_i-f(x_i,\theta)\right)^2\right) $$

Entonces, podemos definir una nueva función:

$$ \chi_{(\theta)}^2 = \sum_i \ \frac{\left(y_i-f(x_i,\theta)\right)^2}{\sigma_i^2} $$

y así podemos calcular la probabilidad de los parámetros $$ P(\theta\ | D I) \propto \exp \left(-\frac{\chi^2_{(\theta)}}{2}\right) $$

Con estas definiciones, resulta inmediato ver que los valores más probables para los parámetros, dada la medición, son los que hacen que $\chi_{(\theta)}^2$ sea mínimo.
Esto le da nombre a la estrategia que se conoce como el método de cuadrados mínimos.

Cómo se calcula el valor de $\theta$ que minimiza el $\chi^2$ ? Lo mejor es dejar eso a las computadoras para el caso general.

Este problema tiene soluciones analíticas siempre que la función sea lineal en los parámetros de ajuste.

El ajuste por una recta, es algo que usamos tan a menudo, que conviene saberlo bien y verlo con más detalle.

2. Ajustando una recta

Consideremos el caso en que la funcion $$f(x,\theta) = a + b x$$ donde $\theta$ encapsula los parámetros $a,b$. Los datos son los que les entregamos en clase.

Y en su versión digitalizada:

In [1]:
%pylab inline
## Cargo los valores
x,y = loadtxt("datos-clase.png.dat", unpack=True)
s = ones_like(x) * 0.01

## Dibujito
figure(figsize=(8,5))
xlim(0.0,1.0)
ylim(0.98,1.12)
title('Datos',fontsize=15)
xlabel('X',fontsize=15)
ylabel('Y',fontsize=15)
plot(x, y,  "ok")
errorbar(x, y, s, linestyle="None")
grid()

Populating the interactive namespace from numpy and matplotlib

Para ajustar una recta, tenemos que estudiar la funcion $\chi^2$ que en este caso es:

$$ \chi_{(a,b)}^2 = \sum_i \ \frac{\left(y_i-(a + b x_i)\right)^2}{\sigma_i^2}. $$

Una cosa interesante a notar es que, dado un conjunto de datos, la funcion $\chi^2$ solo depende de $a$ y $b$. De esta forma podemos utilizar cualquier algoritmo de minimizacion. Probemos diferentes estrategias.

2.1 ¿Cómo se ve la funcion $\chi^2(a,b)$?

In [2]:
# Defino la sumatoria que genera Chi2
def Chi2(a, b):
    out = ((y - (a +  b * x))/s)**2
    return out.sum()
In [3]:
from mpl_toolkits.mplot3d import Axes3D
#from matplotlib import cm
#from matplotlib.ticker import LinearLocator, FormatStrFormatter

######################################
zoom = 0.2
a_true = 1.0
b_true = 0.1
Npun   = 91  ## conviene que sea impar 
######################################

## Grilla de datos
ag = linspace(a_true-zoom, a_true+zoom, Npun)
bg = linspace(b_true-zoom, b_true+zoom, Npun)

am, bm = meshgrid(ag, bg) #parece que esta funcion es la papa para hacer graficos 3D

Z = zeros((am.shape))
for ix in range(Z.shape[0]):
    for iy in range(Z.shape[1]):
        Z[ix, iy] = Chi2(am[ix, iy], bm[ix, iy])

## Gráfico 3D        

#### 3D
figure(figsize=(18,6))
ax=subplot(121, projection='3d')
ax.plot_surface(am, bm, Z / 10000)
title('$\chi^2 / 10000$',fontsize=15)
xlabel('a',fontsize=15)
ylabel('b',fontsize=15)       

#### Contornos
subplot(122)
levels = [ 0.05, 0.1, 0.25, 0.5, 1, 2.5, 5, 10, 25, 50, 100, 250, 500, 1000, 2500, 5000, 10000]
CS = contour(am, bm, Z/10000, levels)
title('$\chi^2 / 10000$',fontsize=15)
xlabel('a',fontsize=15)
ylabel('b',fontsize=15)
clabel(CS, inline=1, fontsize=8);

2.2 Busco el mínimo de $\chi^2(a,b)$ usando fuerza bruta

Simplemente inspecciono en qué lugar de la grilla está el valor mínimo de $\chi^2$.
Anda muy bien ! (porque la grilla tiene pasos muy chiquitos)

In [4]:
print("Chi2_min_fbruta    =",    Z[Z==Z.min()][0])
print("   a_min_fbruta    =",   am[Z==Z.min()][0])
print("   b_min_fbruta    =",   bm[Z==Z.min()][0])

Chi2_min_fbruta    = 15.417769405
   a_min_fbruta    = 0.995555555556
   b_min_fbruta    = 0.108888888889

2.3 Busco el mínimo de $\chi^2(a,b)$ usando funciones de librería

Las funciones de librería andan casi siempre fenómeno.

Una ventaja es que devuelven la inversa del Hessiano (la matriz de derivadas segundas lo que nos permite calcular los errores de los parámetros de la siguiente manera. Sabemos que $$ P(\theta\ | DI) \propto \exp \left(-\frac{\chi^2_{(\theta)}}{2}\right) $$

pero cerca del extremo se puede aproximar $$ \chi^2_{(\theta)} \approx \chi^2_{(\hat{\theta})} + \frac{1}{2} (\theta-\hat{\theta})^T\ \mathbf{H} \,(\theta-\hat{\theta}) $$

La probabilidad de los parámetros puede aproximarse como una distribución normal multivariada, cerca del extremo: $$ P(\theta\ | DI) \approx \exp \left(-\frac{1}{2} (\theta-\hat{\theta})^T\ \mathbf{\Sigma}^{-1} \,(\theta-\hat{\theta}) \right) $$

donde $\mathbf{\Sigma}$ es la matriz de correlación de los parámetros.

Estas igualdades nos permiten identificar que $$ \mathbf{\Sigma} = 2\ \mathbf{H}^{-1} = \left(\begin{array}{cc} \sigma_a^2 & \sigma_a \sigma_b \rho \\ \sigma_a \sigma_b \rho & \sigma_b^2 \end{array}\right) $$

In [5]:
# Cargo la libreria de minimizacion
from scipy.optimize import minimize

# Hago la minimizacion
result = minimize(lambda t: Chi2(*t), x0 = [0,0])

# Impresión de resultados
print("Resultado del la minimizacion")
print()
print(result)
S = 2 * result['hess_inv']
sa = sqrt(S[0,0])
sb = sqrt(S[1,1])
rho= S[0,1]/sa/sb

# Dibujito
figure(figsize=(8,5))
xlim(0.0,1.0)
ylim(0.98,1.12)
title('Datos',fontsize=15)
xlabel('X',fontsize=15)
ylabel('Y',fontsize=15)
plot(x, y,  "ok")
plot(x, result["x"][0] + result["x"][1]*x, "-r")
errorbar(x, y, s, linestyle="None")
grid()
show()

print()
print(" a_min_minimize =", result["x"][0],"+/-", sa)
print(" b_min_minimize =", result["x"][1],"+/-", sb)
print(" correlación    =", rho)

Resultado del la minimizacion

      fun: 15.317199130429994
 hess_inv: array([[  7.54250273e-06,  -1.10902313e-05],
       [ -1.10902313e-05,   2.21909193e-05]])
      jac: array([  0.00000000e+00,  -8.10623169e-06])
  message: 'Optimization terminated successfully.'
     nfev: 28
      nit: 4
     njev: 7
   status: 0
  success: True
        x: array([ 0.99560243,  0.10990709])

 a_min_minimize = 0.99560242825 +/- 0.0038839420001
 b_min_minimize = 0.109907087315 +/- 0.00666196956781
 correlación    = -0.857225773184

2.4 Minimizacion de $\chi^2(a,b)$ por el metodo clásico (algo de algebra y algo para el camino)

Busquemos la condicion de minimo de $\chi^2(a,b)$ por el metodo clásico, es decir encontremos la condicion en que las derivadas de $\chi^2(a,b)$ sean nulas:
\

\begin{align} \frac{\partial \chi^2(\hat{a},\hat{b})}{\partial a} & = & -2\ \sum \frac{y_i-(a+bx_i)}{\sigma^2_i} & = & 0, \\ \frac{\partial \chi^2(\hat{a},\hat{b})}{\partial b} & = & -2\ \sum x_i\ \frac{y_i-(a+bx_i)}{\sigma^2_i} & = & 0. \end{align}

Haciendo toda el álgebra encontramos que:

\begin{align} \hat{a} & \sum 1/\sigma^2_i & + \hat{b} & \sum x_i/\sigma^2_i & = & \sum y_i/\sigma^2_i, \\ \hat{a} & \sum x_i/\sigma^2_i & + \hat{b} & \sum x^2_i/\sigma^2_i & = & \sum y_i x_i/\sigma^2_i, \end{align}

por lo que

\begin{align} \Delta & = \begin{vmatrix} \sum 1/\sigma^2_i & \sum x_i/\sigma^2_i \\ \sum x_i/\sigma^2_i & \sum x^2_i/\sigma^2_i \end{vmatrix}, \\ \\ \hat{b} & = \frac{\begin{vmatrix} \sum 1/\sigma^2_i & \sum y_i/\sigma^2_i \\ \sum x_i\sigma^2_i & \sum x_i y_i/\sigma^2_i \end{vmatrix}}{\Delta}, \\ \\ \hat{a} & = \frac{\begin{vmatrix} \sum y_i/\sigma^2_i & \sum x_i/\sigma^2_i \\ \sum x_i y_i/\sigma^2_i & \sum x^2_i/\sigma^2_i \end{vmatrix}}{\Delta}. \\ \\ \end{align} En el caso en que $\sigma_i = \sigma$ se simplifica enormemente quedado la conocida relación: \begin{align} \\ \hat{b} & = \frac{ \sum{x_i y_i} - \frac{1}{n} \sum x_i \sum y_j}{ \sum {x_i^2} - \frac{1}{n} (\sum x_i)^2 }, \\ \\ & = \frac{ Cov[x, y] }{ Var[x] }, \\ \\ \hat{a} & = \bar{y} - \hat{b}\,\bar{x}. \end{align}

Sé que esto es horrible, pero al menos ya está hecho. Sólo aclaramos que $$ Cov[x,y] = \frac{1}{n} \sum (x_i-\overline{x})(y_i-\overline{y}) = \sigma_x \sigma_y \ \rho $$

Parece mentira las cosas que veo:

  • $\Delta$ sólo depende de los valores de $x$
  • La recta pasa por el baricentro de los datos (mirar fijo la última expresión)
  • Si todos los $\sigma_i$ son iguales, $\hat{a}$ y $\hat{b}$ no dependen de los errores de medición

Aprovechemos la situación para hacer las cuentas en nuestro ejemplo:

In [6]:
xy = x * y
xx = x * x

b_clasica = (xy.mean()-x.mean()*y.mean())/(xx.mean()-x.mean()**2)
a_clasica = y.mean() - b_clasica * x.mean()

print()
print(" a_min_clasica =", a_clasica)
print(" b_min_clasica =", b_clasica)

 a_min_clasica = 0.995602442241
 b_min_clasica = 0.109907074228

2.5 ¿Que podemos decir sobre la dispersión de $a$ y $b$?

Podemos calcular la matriz de las derivadas segunda del $\chi^2$ y resulta que:

$$ \mathbf{H}_{\chi^2} = 2 \begin{pmatrix} \sum 1/\sigma^2_i & \sum x_i/\sigma^2_i \\ \sum x_i/\sigma^2_i & \sum x^2_i/\sigma^2_i \end{pmatrix} $$

y ya utilizamos que la matriz de coviariancia es:

$$ \mathbf{\Sigma} = 2\ \mathbf{H}^{-1} = \frac{1}{\Delta} \ \begin{pmatrix} \sum x^2_i/\sigma^2_i & -\sum x_i/\sigma^2_i \\ -\sum x_i/\sigma^2_i & \sum 1/\sigma^2_i \end{pmatrix} $$

de manera que

$$ \left\{ \begin{eqnarray} \sigma_a^2 & = & + \frac{1}{\Delta} & \sum \frac{x^2_i}{\sigma^2_i} \\ \sigma_b^2 & = & + \frac{1}{\Delta} & \sum \frac{1}{\sigma^2_i} \\ \sigma_a \sigma_b \rho & = & -\frac{1}{\Delta} & \sum \frac{x_i}{\sigma^2_i} \end{eqnarray}\right. $$

La expresión de la correlación es un resultado más que interesante, ya que está solamente relacionada a los valores de $x_i$. Considerando que $\Delta$ siempre mayor que cero se ve que el valor de la covarianza se halla relacionado con el baricentro de los datos representado por $\sum x_i/\sigma_i^2$.

Esto quiere decir que datos de $x$ mayoritariamente positivos, dan una covarianza de $a$ y $b$ negativa, mientras que si los $x$ son negativos, la covarianza será positiva.

En un caso más balanceado en el cual $\sum x_i/\sigma_i^2 \approx 0$ la variaciones de $a$ y $b$ se encontraran prácticamente descorrelacionadas.

In [7]:
ls2  = 1.0/s**2
xs2  = x/s**2
x2s2 = x**2/s**2

Delta = ls2.sum()*x2s2.sum() - xs2.sum()**2
lDelta = 1.0/Delta

var_a_clasica = lDelta * x2s2.sum()
var_b_clasica = lDelta * ls2.sum()
cov_ab_clasica= -1.0 * lDelta * xs2.sum()

s_a_clasica = sqrt(var_a_clasica)
s_b_clasica = sqrt(var_b_clasica)

print()
print(" a_min_clasica =", a_clasica,"+/-", s_a_clasica)
print(" b_min_clasica =", b_clasica,"+/-", s_b_clasica)
print(" correlación   =", cov_ab_clasica/s_a_clasica/s_b_clasica)

 a_min_clasica = 0.995602442241 +/- 0.00388394202301
 b_min_clasica = 0.109907074228 +/- 0.00666196955447
 correlación   = -0.857225774717

3. A mano, como estamos

¿Como se ven los ajustes que hicimos en clase?

Si superponemos todas las rectas que dibujaron quedan muy bonitas:

Vemos que:

  • a pesar de que la estrategia de ajuste es imprecisa y subjetiva, los resultados son muy consistentes
  • las rectas se agrupan en el baricentro de los datos, tal como predice nuestro cálculo

¿y los valores de a y b?

In [8]:
# Cargo los valores
from scipy.stats import pearsonr
b_clase, a_clase = loadtxt("datos-fit.dat", unpack=True)

figure(figsize=(8,8))
levels = [ 0.01, 0.025, 0.05, 0.1, 0.25]

xlim([0.95,1.05])
ylim([0.05,0.15])
axes().set_aspect("equal")
CS = contour(am, bm, Z/10000, levels)
clabel(CS, inline=1, fontsize=8)
plot(a_clase, b_clase, "or",label='clase',alpha=.7)
plot(a_clasica, b_clasica, "ob",label='min',markersize=15)
title('$\chi^2 / 10000$',fontsize=20)
xlabel('a',fontsize=20)
ylabel('b',fontsize=20)
legend()
grid()

print()
print(" a_min_clase =", a_clase.mean(),"+/-", a_clase.std())
print(" b_min_clase =", b_clase.mean(),"+/-", b_clase.std())
print(" correlación  =", pearsonr(a_clase,b_clase)[0])

 a_min_clase = 0.993094594595 +/- 0.00345159155972
 b_min_clase = 0.112707837838 +/- 0.00855935686367
 correlación  = -0.506618260732

Vemos que los valores en la estimación de los parámetros y los errores son prácticamente idénticos. La correlación tiene el mismo signo, pero es ligeramente menor.

_nunca sabemos hasta dónde sabemos . . ._