¿Cuál es el error de numpy.polyfit?

Quiero usar numpy.polyfit para cálculos físicos, por lo tanto, necesito la magnitud del error.

Si especifica full=True en su llamada a polyfit , incluirá información adicional:

 >>> x = np.arange(100) >>> y = x**2 + 3*x + 5 + np.random.rand(100) >>> np.polyfit(x, y, 2) array([ 0.99995888, 3.00221219, 5.56776641]) >>> np.polyfit(x, y, 2, full=True) (array([ 0.99995888, 3.00221219, 5.56776641]), # coefficients array([ 7.19260721]), # residuals 3, # rank array([ 11.87708199, 3.5299267 , 0.52876389]), # singular values 2.2204460492503131e-14) # conditioning threshold 

El valor residual devuelto es la sum de los cuadrados de los errores de ajuste, sin saber si esto es lo que está buscando:

 >>> np.sum((np.polyval(np.polyfit(x, y, 2), x) - y)**2) 7.1926072073491056 

En la versión 1.7 también hay una palabra clave cov que devolverá la matriz de covarianza para sus coeficientes, que podría utilizar para calcular la incertidumbre de los coeficientes de ajuste.

Como se puede ver en la documentación :

 Returns ------- p : ndarray, shape (M,) or (M, K) Polynomial coefficients, highest power first. If `y` was 2-D, the coefficients for `k`-th data set are in ``p[:,k]``. residuals, rank, singular_values, rcond : present only if `full` = True Residuals of the least-squares fit, the effective rank of the scaled Vandermonde coefficient matrix, its singular values, and the specified value of `rcond`. For more details, see `linalg.lstsq`. 

Lo que significa que si puedes hacer un ajuste y obtener los residuos como:

  import numpy as np x = np.arange(10) y = x**2 -3*x + np.random.random(10) p, res, _, _, _ = numpy.polyfit(x, y, deg, full=True) 

Luego, la p son sus parámetros de ajuste, y la res serán los residuales, como se describió anteriormente. Los _ son porque no necesita guardar los últimos tres parámetros, por lo que puede guardarlos en la variable _ que no usará. Esto es una convención y no es obligatorio.


La respuesta de @ Jaime explica lo que significa el residual. Otra cosa que puedes hacer es mirar esas desviaciones al cuadrado como una función (la sum de las cuales es res ). Esto es particularmente útil para ver una tendencia que no encaja lo suficiente. res puede ser grande debido al ruido estadístico, o posiblemente a la mala adaptación sistemática, por ejemplo:

 x = np.arange(100) y = 1000*np.sqrt(x) + x**2 - 10*x + 500*np.random.random(100) - 250 p = np.polyfit(x,y,2) # insufficient degree to include sqrt yfit = np.polyval(p,x) figure() plot(x,y, label='data') plot(x,yfit, label='fit') plot(x,yfit-y, label='var') 

Así que en la figura, note el mal ajuste cerca de x = 0 :
Polyfit