python - power - En Scipy, ¿cómo y por qué curve_fit calcula la covarianza de las estimaciones de los parámetros?

optimal parameters not found: number of calls to function has reached maxfev=1000.

Question

He estado usando scipy.optimize.leastsq para ajustar algunos datos. Me gustaría obtener algunos intervalos de confianza en estas estimaciones, así que veo la salida de cov_x , pero la documentación no es muy clara en cuanto a qué es esto y cómo obtener la matriz de covarianza para mis parámetros.

En primer lugar, dice que es un jacobiano, pero en las notes también dice que " cov_x es una aproximación jacobiana al hessiano", de modo que no es realmente un jacobiano sino un hessiano que utiliza una aproximación del jacobiano. ¿Cuál de estas afirmaciones es correcta?

En segundo lugar esta frase para mí es confusa:

Esta matriz se debe multiplicar por la varianza residual para obtener la covarianza de las estimaciones del parámetro, consulte curve_fit .

De hecho, voy a mirar el código fuente de curve_fit donde lo hacen:

s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0))
pcov = pcov * s_sq

que corresponde a multiplicar cov_x por s_sq pero no puedo encontrar esta ecuación en ninguna referencia. ¿Alguien puede explicar por qué esta ecuación es correcta? Mi intuición me dice que debería ser al revés, ya que se supone que cov_x es un derivado (jacobiano o hessiano), así que estaba pensando: cov_x * covariance(parameters) = sum of errors(residuals) donde sigma(parameters) es el cosa que quiero

¿Cómo conecto la cosa que curve_fit está haciendo con lo que veo en, por ejemplo? wikipedia: http://en.wikipedia.org/wiki/Propagation_of_uncertainty#Non-linear_combinations

Answer 1

Encontré esta solución durante mi búsqueda de una pregunta similar, y solo tengo una pequeña mejora en la respuesta de HansHarhoff. La salida completa de leastsq proporciona un valor de retorno infodict, que contiene infodict [''fvec''] = f (x) -y. Por lo tanto, para calcular el chi cuadrado reducido = (en la notación anterior)

s_sq = (infodict[''fvec'']**2).sum()/ (Nn)

Por cierto Gracias HansHarhoff por hacer la mayor parte del trabajo pesado para resolver esto.

Answer 2

OK, creo que encontré la respuesta. Primero, la solución: cov_x * s_sq es simplemente la covarianza de los parámetros que es lo que desea. Tomar el cuadrado de los elementos diagonales le dará una desviación estándar (¡pero tenga cuidado con las covarianzas!).

Varianza residual = chi cuadrado reducido = s_sq = suma [(f (x) -y) ^ 2] / (Nn), donde N es el número de puntos de datos yn es el número de parámetros de ajuste. Plaza de chi reducida .

El motivo de mi confusión es que cov_x, tal como se indica con leastsq, no es en realidad lo que se llama cov (x) en otros lugares, sino que es cov reducido (x) o cov fraccional (x). La razón por la que no aparece en ninguna de las otras referencias es que se trata de un cambio de escala simple que es útil en los cálculos numéricos, pero no es relevante para un libro de texto.

Acerca de Hessian versus Jacobian, la documentación está mal redactada. Es la arpillera la que se calcula en ambos casos, como es obvio, ya que el jacobiano es cero como mínimo. Lo que quieren decir es que están utilizando una aproximación al jacobiano para encontrar al Hessiano.

Una nota más. Parece que el resultado de curve_fit no tiene en cuenta el tamaño absoluto de los errores, sino que solo tiene en cuenta el tamaño relativo de los sigmas proporcionados. Esto significa que el pcov devuelto no cambia incluso si las barras de error cambian por un factor de un millón. Esto, por supuesto, no es correcto, pero parece ser una práctica estándar, es decir. Matlab hace lo mismo cuando usa su caja de herramientas de ajuste de curva. El procedimiento correcto se describe aquí: https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Parameter_errors_and_correlation

Parece bastante sencillo hacer esto una vez que se ha encontrado el óptimo, al menos para los cuadrados mínimos lineales.