¿Cómo interpreto la matriz de covarianza de un ajuste de curva?

Question

No soy muy bueno en estadística, así que pido disculpas si esta es una pregunta simplista. Estoy ajustando una curva a algunos datos, y a veces mis datos se ajustan mejor a una exponencial negativa en la forma $a * e^{(-b * x)} + c$, y a veces el ajuste se acerca más a $a * e^{(-b * x^2)} + c$. Sin embargo, a veces ambos fallan, y me gustaría recurrir a un ajuste lineal. Mi pregunta es, ¿cómo puedo determinar qué modelo se ajusta mejor a un conjunto de datos en particular a partir de la matriz de varianza-covarianza resultante que se devuelve de la función scipy.optimize.curve_fit()? Creo que la varianza está en una de las diagonales de esta matriz, pero no estoy seguro de cómo interpretar eso.

ACTUALIZACIÓN: Basándome en una pregunta similar, espero que la matriz de varianza-covarianza pueda decirme cuál de los tres modelos que estoy intentando se ajusta mejor a los datos (estoy intentando ajustar muchos conjuntos de datos a uno de estos tres modelos).

Las matrices resultantes se ven así para el ejemplo dado:

pcov_lin 
[[  2.02186921e-05  -2.02186920e-04]
 [ -2.02186920e-04   2.76322124e-03]]
pcov_exp
[[  9.05390292e+00  -7.76201283e-02  -9.20475334e+00]
 [ -7.76201283e-02   6.69727245e-04   7.90218415e-02]
 [ -9.20475334e+00   7.90218415e-02   9.36160310e+00]]
pcov_exp_2 
[[  1.38338049e-03  -7.39204594e-07  -7.81208814e-04]
 [ -7.39204594e-07   8.99295434e-09   1.92970700e-06]
 [ -7.81208814e-04   1.92970700e-06   9.14746758e-04]]

Aquí hay un ejemplo de lo que estoy haciendo:

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
import scipy.optimize

def exp_func(x, a, b, c):
    return a * np.exp(-b * x) + c

def exp_squared_func(x, a, b, c):
    return a * np.exp(-b * x*x*x) + c

def linear_func(x, a, b):
    return a*x + b

def main():
    x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], np.float)
    y = np.array([1, 1, 1, 1, 0.805621, 0.798992, 0.84231, 0.728796, 0.819471, 0.570414, 0.355124, 0.276447, 0.159058, 0.0762189, 0.0167807, 0.0118647, 0.000319948, 0.00118267, 0, 0, 0], np.float)

    p0 = [0.7746042467213462, 0.10347274384077858, -0.016253458007293588]
    popt_lin, pcov_lin      = scipy.optimize.curve_fit(linear_func, x, y)
    popt_exp, pcov_exp      = scipy.optimize.curve_fit(exp_func, x, y)
    popt_exp_2, pcov_exp_2  = scipy.optimize.curve_fit(exp_squared_func, x, y)

    plt.figure()
    plt.plot(x, y, 'ko', label="Datos originales")
    plt.plot(x, linear_func(x, *popt_lin), 'r-', label='lineal')
    plt.plot(x, exp_func(x, *popt_exp), 'b-', label='exponencial')
    plt.plot(x, exp_squared_func(x, *popt_exp_2), 'g-', label='exponencial cuadrada')
    plt.legend()
    plt.show()

if __name__ == '__main__':
    main()

Answer 1

Como aclaración, la variable pcov de scipy.optimize.curve_fit es la covarianza estimada de la estimación del parámetro, es decir, en pocas palabras, dada la data y un modelo, cuánta información hay en la data para determinar el valor de un parámetro en el modelo dado. Por lo tanto, realmente no te dice si el modelo elegido es bueno o no. Ver también esto.

El problema de qué es un buen modelo es, de hecho, un problema difícil. Como argumentan los estadísticos

Todos los modelos son incorrectos, pero algunos son útiles

Por lo tanto, los criterios a utilizar en la comparación de diferentes modelos dependen de lo que se quiera lograr.

Por ejemplo, si se desea una curva que esté lo más "cerca posible" de la data, se podría seleccionar un modelo que dé el menor residual. En su caso, sería el modelo func y los parámetros estimados popt que tengan el valor más bajo al calcular

numpy.linalg.norm(y-func(x, *popt))

Sin embargo, si se selecciona un modelo con más parámetros, el residual disminuirá automáticamente, a costa de una mayor complejidad del modelo. Entonces, todo vuelve a depender de cuál sea el objetivo del modelo.