Python – calculando líneas de tendencia con errores

Así que tengo algunos datos almacenados en dos listas, y los graficé usando

plot(datasetx, datasety) 

Entonces establezco una línea de tendencia

 trend = polyfit(datasetx, datasety) trendx = [] trendy = [] for a in range(datasetx[0], (datasetx[-1]+1)): trendx.append(a) trendy.append(trend[0]*a**2 + trend[1]*a + trend[2]) plot(trendx, trendy) 

Pero tengo una tercera lista de datos, que es el error en la seguridad de datos original. Estoy bien con el trazado de las barras de error, pero lo que no sé es usar esto, cómo encontrar el error en los coeficientes de la línea de tendencia polinomial.

Entonces, digamos que mi línea de tendencia resultó ser 5x ^ 2 + 3x + 4 = y, tiene que haber algún tipo de error en los valores 5, 3 y 4.

¿Hay alguna herramienta que use NumPy que calcule esto para mí?

Creo que puedes usar la función curve_fit de scipy.optimize ( documentación ). Un ejemplo básico del uso:

 import numpy as np from scipy.optimize import curve_fit def func(x, a, b, c): return a*x**2 + b*x + c x = np.linspace(0,4,50) y = func(x, 5, 3, 4) yn = y + 0.2*np.random.normal(size=len(x)) popt, pcov = curve_fit(func, x, yn) 

Siguiendo la documentación, pcov da:

La covarianza estimada del popt. Las diagonales proporcionan la varianza de la estimación del parámetro.

Entonces, de esta manera usted puede calcular una estimación de error en los coeficientes. Para tener la desviación estándar puede tomar la raíz cuadrada de la varianza.

Ahora tiene un error en los coeficientes, pero solo se basa en la desviación entre los datos y el ajuste. En caso de que también quiera tener en cuenta un error en los propios datos y, la función curve_fit proporciona el argumento sigma :

sigma: Ninguna o secuencia de longitud N

Si no es Ninguno, representa la desviación estándar de ydata. Este vector, si se da, se usará como ponderaciones en el problema de los mínimos cuadrados.

Un ejemplo completo:

 import numpy as np from scipy.optimize import curve_fit def func(x, a, b, c): return a*x**2 + b*x + c x = np.linspace(0,4,20) y = func(x, 5, 3, 4) # generate noisy ydata yn = y + 0.2 * y * np.random.normal(size=len(x)) # generate error on ydata y_sigma = 0.2 * y * np.random.normal(size=len(x)) popt, pcov = curve_fit(func, x, yn, sigma = y_sigma) # plot import matplotlib.pyplot as plt fig = plt.figure() ax = fig.add_subplot(111) ax.errorbar(x, yn, yerr = y_sigma, fmt = 'o') ax.plot(x, np.polyval(popt, x), '-') ax.text(0.5, 100, r"a = {0:.3f} +/- {1:.3f}".format(popt[0], pcov[0,0]**0.5)) ax.text(0.5, 90, r"b = {0:.3f} +/- {1:.3f}".format(popt[1], pcov[1,1]**0.5)) ax.text(0.5, 80, r"c = {0:.3f} +/- {1:.3f}".format(popt[2], pcov[2,2]**0.5)) ax.grid() plt.show() 

resultado


Luego algo más , sobre el uso de matrices numpy. Una de las principales ventajas de usar numpy es que puede evitar los bucles porque las operaciones en las matrices se aplican elementwise. Por lo tanto, el for-loop en su ejemplo también se puede hacer como sigue

 trendx = arange(datasetx[0], (datasetx[-1]+1)) trendy = trend[0]*trendx**2 + trend[1]*trendx + trend[2] 

Donde uso arange lugar de rango, ya que devuelve una matriz numpy en lugar de una lista. En este caso también puedes usar la función polyval :

 trendy = polyval(trend, trendx) 

No he podido encontrar ninguna forma de obtener los errores en los coeficientes incorporados para numpy o python. Tengo una herramienta simple que escribí basada en la Sección 8.5 y 8.6 de An Introduction to Error Analysis de John Taylor. Quizás esto sea suficiente para su tarea (tenga en cuenta que el rendimiento predeterminado es la varianza, no la desviación estándar). Puede obtener grandes errores (como en el ejemplo proporcionado) debido a una covarianza significativa.

 def leastSquares(xMat, yMat): ''' Purpose ------- Perform least squares using the procedure outlined in 8.5 and 8.6 of Taylor, solving matrix equation X a = Y Examples -------- >>> from scipy import matrix >>> xMat = matrix([[ 1, 5, 25], [ 1, 7, 49], [ 1, 9, 81], [ 1, 11, 121]]) >>> # matrix has rows of format [constant, x, x^2] >>> yMat = matrix([[142], [168], [211], [251]]) >>> a, varCoef, yRes = leastSquares(xMat, yMat) >>> # a is a column matrix, holding the three coefficients a, b, c, corresponding to >>> # the equation a + b*x + c*x^2 Returns ------- a: matrix best fit coefficients varCoef: matrix variance of derived coefficents yRes: matrix y-residuals of fit ''' xMatSize = xMat.shape numMeas = xMatSize[0] numVars = xMatSize[1] xxMat = xMat.T * xMat xyMat = xMat.T * yMat xxMatI = xxMat.I aMat = xxMatI * xyMat yAvgMat = xMat * aMat yRes = yMat - yAvgMat var = (yRes.T * yRes) / (numMeas - numVars) varCoef = xxMatI.diagonal() * var[0, 0] return aMat, varCoef, yRes