Episodio III

Willy Pregliasco y Mariano Gómez Berisso, Ago 2016

Comparación de diferentes modelos

Vamos a considerar que tenemos una medición y varios modelos $M_i$ para contrastar con los datos $D$.

Lo que nos interesa es ver qué nos dicen los datos acerca de cuál modelo es el que está mejor fundamentado por las mediciones.

Las ideas que voy a aplicar ideas que están descriptas en el libro de Phil Gregory (cap 3).

Finalmente, vamos a realizar la cuenta para hacer un ajuste polinómico y evaluando el grado del polinomio que hace falta usar para describir los datos.

1. Comparación de modelos y factor de Occam

Si tengo varios modelos para elegir mutuamente excluyentes $M_i$ y mediciones de dator $D$ con información contextual $I$, el Teorema de Bayes nos dice que:

$$ P(M_i\ | DI) = \frac{P(D\ |M_i I)}{P(D\ |I)} \ P(M_i\ |I) $$

El denominador siempre molesta, pero es el mismo para todos los modelos, y como los modelos son exhaustivos:

$$ \sum_i P(M_i\ | DI) = 1 $$

de manera que podemos calcular obviando el factor de normalización y finalmente normalizando:

$$ P(M_i\ | DI) \propto P(D\ |M_i I) \ P(M_i\ |I) $$

El segundo factor representa el valor prior de cada modelo. Su influencia es evidente y no es interesante considerarla en esta instancia. Vamos a suponer, sin hacer mucho lío, que todos los modelos son igualmente probables antes de considerar los datos.

1.1 Marginalización

El problema que tiene esta expresión, es que cada modelo tiene parámetros $\theta_i$ para estimar.
Cuando estudiamos los ajustes, vimos que los datos determinan una distribución de los parámetros, y eso nos permite marginalizar la expresión:

$$ P(D\ |M_i I) = \int P(D\ |M_i \theta_i I)\ P(\theta_i\ |I) \ \mathrm{d}\theta_i = \int {\cal{L}}_i(\theta_i) \ \mathrm{d}\theta_i $$

Lo bonito de esta expresión es que esta integral puedo calcularla porque mi modelo me dice cuál es la probabilidad de los datos, dado un conjunto de parámetros.

1.2 Cálculo

Para calcular la integral, es preciso establecer una distribución prior de los parámetros. En este caso haremos una distribución ancha y uniforme, pero se puede calcular poniendo una función con la forma que sea necesaria (y que será ancha y generosa).

Si a los datos les asignamos una distribución normal alrededor de la curva de ajuste:

$$ {\cal{L}}_i(\theta_i) = P(D\ |M_i\theta_iI) P(\theta_i\ |I) = \left( \prod_k \mathrm{Gauss}\left( y_k, f(x_k, \theta_i \right), \sigma) \right) \ P(\theta_i\ |I) $$

Y con esto estaría en condiciones de calcular lo que hace falta

1.3 Factor de Occam

Sin hacer la integral, podemos ver una propiedad cualitativa de este resultado.

Como dijimos antes, $P(\theta_i\ |I)$ es la distribución prior de los parámetros, que es una función ancha y aproximadamente constante sobre la integral. El ancho representa la incerteza que tengo sobre los parámetros antes de hacer el ajuste. Vamos a llamarla $\Delta \theta_i$. Como la prior es una densidad de probabilidad en $\theta_i$, tiene que estar normalizada, y se puede aproximar $$ P(\theta_i\ |I) \approx 1/\Delta\theta_i $$ y sacarlo fuera de la integral.

El integrando que queda ya no está normalizado. Como función de $\theta_i$ es una likelihood que tiene un máximo y un ancho $\delta\theta_i$ que es el ancho de los parámetros que obtengo de considerar los datos. De estas consideraciones resulta que:

$$ P(D\ |M_i I) = \cal{L}_{i - max}\ \left(\frac{\delta\theta_i}{\Delta\theta_i}\right) $$

El primer factor habla de la bondad del ajuste y será mayor cuanto mejor ajuste la curva a los datos experimentales. Pero el segundo factor es un número muy chiquito, que penaliza la existencia misma del parámetro. Ese cociente se llama factor de Occam porque requiere que el likelihood mejore mucho para que se justifique la inclusión de nuevos parámetros. Vale aclarar que si $\theta_i$ encapsula varios parámetros, ambos delta se refieren al volumen en el espacio de los parámetros. Y el cociente de Occam da más chico cuanta mayor sea la dimensión en ese espacio.

_entia non sunt multiplicanda praeter necessitatem
(no hay que multiplicar los entes sin necesidad)_
Guillermo de Occam

2. Ejemplo: ajuste polinómico

Vamos a hacer un ajuste polinómico de datos del calentamiento global. Se mide la diferencia de temperaturas respecto de un valor promedio de referencia en varias estaciones meteorológicas y se hace una extrapolación a escala mundial. Ese valor se denomina la anomalía de temperatura. Hay tablas con datos anuales. Elegimos aquí un subconjunto de datos cada 5 años desde 1990.

Los datos tienen una variación que es menor que 0.1 $C$, de manera que los datos deberían desviarse de la curva del orden de $\sigma\approx .07$. Les vamos a asignar un error con distribución normal.

In [1]:
%pylab inline

### Datos del calentamiento global http://climate.nasa.gov/vital-signs/global-temperature/
anio = [1990,1995,2000,2005,2010,2015]
anomalia = [0.44,0.46,0.42,0.69,0.72,0.87]
sigma = .07

### Necesito menos datos
X, Y = array(anio), array(anomalia)

### Gráfico
figure(figsize=(10,4))
plot(X,Y,'ko', markersize=8)
errorbar(X, Y, sigma, linestyle="None")
xlabel('año',fontsize=20)
ylabel('anomalia [C]',fontsize=20)
title('Global Warming',fontsize=20)
grid()
xlim(1985,2020)
ylim(.3,1);

Populating the interactive namespace from numpy and matplotlib

2.1 Priors

Cómo estimo las distribuciones prior de los parámetros del polinomio? Este es un tema que requiere cierta delicadeza, pero que no lleva nada muy conceptual. Puede esquivárselo sin misericordia ni remordimiento.

Vamos a elegir priors uniformes y simétricas respecto del origen.

Los datos varían entre $a_0 = 0.4\ C$ y $a_1=0.9\ C$ entre $t_0=1990$ y $t_1=2020$ y uno espera que los términos del polinomio representen esas variaciones. El problema es que puedo tener grandes variaciones en un término que estén compensadas por grandes variaciones en la dirección opuesta hechas por otro término. Este es un comportamiento que tengo que permitir, aunque poniéndole alguna limitación.

Otro inconveniente es que es difícil imaginarse el valor de los coeficientes del polinomio, a menos que centre el origen de los valores de $x$ en el centro de los datos $(x_0, y_0)$. En ese caso, la variación de $x$ alrededor de $x_0$ será menor que $\pm\Delta x$ y lo mismo pasará en la otra dirección.

La variación de un término no debería superar unas $\beta=10$ veces la variación total de la función:

$$ \left| \theta'_n \right| \ (\Delta x)^n < \beta \ \Delta y $$

y el polinomio tiene la forma:

$$ y-y_0 = \sum_i \theta'_i \ (x-x_0)^i $$

En nuestro caso, los límites resultan:

In [2]:
#Límites
beta = 10
x0, y0 = mean(X)  , mean(Y)
dx, dy = max(X-x0), max(Y-y0)

l0 = beta * dy
l1 = l0 / dx
l2 = l1 / dx
l3 = l2 / dx
l4 = l3 / dx
l5 = l4 / dx

l = array([l0,l1,l2,l3,l4,l5])
print('n=0', l0)
print('n=1', l1)
print('n=2', l2)
print('n=3', l3)
print('n=4', l4)
print('n=5', l5)

Xc, Yc = X-x0, Y-y0

n=0 2.7
n=1 0.216
n=2 0.01728
n=3 0.0013824
n=4 0.000110592
n=5 8.84736e-06

Esto suena muy razonable.

La prior para cada parámetro será una función uniforme de la forma $$ P(\theta_i) = \left \{\begin{array}{ll} \frac{1}{2l_i} & | \theta_i | < l_i \\ 0 & | \theta_i | > l_i \end{array} \right. $$

In [3]:
### Defino las funciones necesarias para integrar

from scipy.optimize import minimize
from scipy.stats    import norm

### Función de ajuste: Polinomio con parámetros tita

def fn(x,tita):
    out=0
    for i,t in enumerate(tita):
        out += t * x**i
    return out

### Logaritmo de la Prior

def lnPrior(tita):
    out=0
    for i,t in enumerate(tita):
        if abs(t)<l[i]:
            out -= log(2*l[i])
        else:
            out -= inf
    return out

### Logaritmo de Likehood

def lnLike(tita,X,Y):  
    ## Quiero poder aceptar un array de tita's
    try:
        iterTest=iter(tita[0])  #test de si es un array
        
    except TypeError:
        # not iterable
        out = 0
        for x,y in zip(X,Y):
            #print (x,y,sigma)
            out += norm.logpdf(y,loc=fn(x,tita),scale=sigma)
        out += lnPrior(tita)
        return out
    
    else:
        # itero sobre los elementos
        out=[]
        for t in tita:
            out.append(lnLike(t,X,Y))
        return array(out)

### Valores Óptimos

def optimo(grado):
    coef  = polyfit(Xc, Yc, grado)  ## Valores iniciales espectacularmente buenos
    t0 = coef[::-1]
    t0n = t0/l[:len(t0)]
    opFn = lambda titan, *arg: -lnLike(titan*l[:len(t0)],*arg)
    
    result = minimize(opFn, x0=t0n, args=(Xc,Yc), method='BFGS', tol =1e-3)
    
    if result['success']:
        return result['x']*l[:len(t0)]
    else:
        print (result)
        return t0*nan

### Chequeo con polinomios de distinto grado
print('Check Norris:')
#set_printoptions(precision=2)
set_printoptions(formatter={'float': '{:+9.2g}'.format})
for g in range(6):
    print('\tgrado',g, optimo(g))
    
set_printoptions(suppress=True)

Check Norris:
	grado 0 [ +1.1e-17]
	grado 1 [ -2.8e-18    +0.018]
	grado 2 [   -0.048    +0.018  +0.00066]
	grado 3 [   -0.048    +0.025  +0.00066  -5.6e-05]
	grado 4 [   -0.049    +0.025  +0.00072  -5.6e-05  -3.3e-07]
	grado 5 [   -0.049     +0.06  +0.00072    -0.001  -3.3e-07  +4.9e-06]

2.2 Integral

La función $\log {\cal{L}(\theta)}$ podemos aproximarla por un paraboloide, de la forma:

$$ \log {\cal{L}(\theta)} \approx \log {\cal{L}_{max}} - \frac{1}{2} \ (\theta -\hat{\theta})^T \ H \ (\theta -\hat{\theta}) $$

donde $H$ es el Hessiano (la matriz de derivadas segundas) de la función $\log {\cal{L}(\theta)}$:

$$ H_{ij} = \left( \frac{\partial^2 }{\partial \theta_i\ \partial \theta_j} \right)\ \log \cal{L}(\theta) $$

Con esta aproximación, la likehood queda proporcional a la densidad de una distribución gaussiana multivariada, y no es difícil ver que

$$ \int {\cal{L}(\theta)} d\theta \approx {\cal{L}_{max}} \ \frac{\sqrt{(2\pi)^k\ }}{\left| H \right|} $$
In [4]:
import numdifftools as nd

def BayesFactor(grado):
    
    ## Valores Óptimos
    tita0 = optimo(grado)
    scale = l[:len(tita0)]
    
    ## Lmax
    Lmax = exp(lnLike(tita0,Xc,Yc))
    
    ## Hessiano (en coordenadas escaleadas)
    func = lambda tin: lnLike(tin*scale,Xc,Yc)
    #func = lambda tin: lnLike(tin,Xc,Yc)
    H = nd.Hessian(func,step=0.0001)(tita0/scale)
    
    ## Cambio de escaleo
    for i,hh in enumerate(H):
        for j,h in enumerate(hh):
            
            H[i,j] /= scale[i]*scale[j]

    ## Salida
    out = Lmax * sqrt((2*pi)**len(tita0)) / linalg.det(H)
    
    return abs(out)

bf = zeros(6)
for g in arange(6):
    bf[g] = BayesFactor(g)
    
## Si los modelos tienen igual prior, los factores tienen que normalizar a 1

bf = array(bf)
bf = bf / sum(bf)
In [5]:
### Grafico resultado

figure(figsize(15,5))
subplot(121)
title(r'$P\ (D\ |M_i I)$', fontsize=20)
xlabel(r'Grado del polinomio: $i$', fontsize=20)
bar(arange(6), bf, align='center',color='r',alpha=.8)
ylim(0,1.15)
grid()

subplot(122)
title(r'$\log P\ (D\ |M_i I)$', fontsize=20)
xlabel(r'Grado del polinomio: $i$', fontsize=20)
bar(arange(6), bf, align='center',color='r',alpha=.8)
#plot(arange(6), bf,'ro')
#plot(arange(6), bf,'r-', alpha=.2)
grid()
yscale('log')
ylim(1e-21,9);
_qué bonitoooo!_
In [6]:
def ploty(grado):
    plot(X,Y,'ko')
    errorbar(X, Y, sigma, linestyle="None")
    xx=linspace(1985,2025,100)
    plot(xx,fn(xx-x0,optimo(grado))+y0,'r-')
    ylim(0,2)
    grid()

figure(figsize=(15,20))

for i in range(6):
    subplot(6,2,i+1)
    title('grado {:d}'.format(i),fontsize=10)
    ploty(i)

3. Aprendizajes

  • Es posible cuantificar la probablidad de cada modelo según los datos medidos.
  • La probabilidad del modelo resulta del compromiso entre un buen ajuste y usar pocos parámetros, lo que queda señalado con los factores de Occam
  • La integral se puede calcular sin problemas, viendo la curvatura del logaritmo en el máximo.
  • El polinomio que mejor ajusta a nuestros datos es una recta con pendiente. Según este criterio, hay evidencia de que la cosa se va calentando en el planeta.
  • Si aumentamos el valor de $\sigma$ de los datos, comienza a ser más probable que el modelo de temperatura constante sea el más apropiado. Eso se debe a que si los errores de la medición son grandes, no es necesario hilar tan fino las variaciones.
  • La pretensión de que un ajuste pase por todos los puntos es
    • innecesaria
    • costosa
    • peligrosas
    • psicótica
    • genera predicciones erróneas
_basta de ajuste ! . . ._