¿Cómo determinar la incertidumbre de los parámetros de ajuste con Python?

Tengo los siguientes datos para xey:

xy 1.71 0.0 1.76 5.0 1.81 10.0 1.86 15.0 1.93 20.0 2.01 25.0 2.09 30.0 2.20 35.0 2.32 40.0 2.47 45.0 2.65 50.0 2.87 55.0 3.16 60.0 3.53 65.0 4.02 70.0 4.69 75.0 5.64 80.0 7.07 85.0 9.35 90.0 13.34 95.0 21.43 100.0 

Para los datos anteriores, estoy tratando de ajustar los datos en el formulario:

fórmula

Sin embargo, hay ciertas incertidumbres asociadas con x e y, donde x tiene una incertidumbre del 50% de x e y tiene una incertidumbre fija. Estoy tratando de determinar la incertidumbre en los parámetros de ajuste con este paquete de incertidumbres . Pero, estoy teniendo problemas con el ajuste de curvas con la función de ajuste de curvas de optimización de scipy. Obtuve el siguiente error:

minpack.error: El resultado de la llamada de función no es una matriz adecuada de flotadores.

¿Cómo puedo corregir el siguiente error y determinar la incertidumbre de los parámetros de ajuste (a, b y n)?

MWE

 from __future__ import division import numpy as np import re from scipy import optimize, interpolate, spatial from scipy.interpolate import UnivariateSpline from uncertainties import unumpy def linear_fit(x, a, b): return a * x + b uncertainty = 0.5 y_error = 1.2 x = np.array([1.71, 1.76, 1.81, 1.86, 1.93, 2.01, 2.09, 2.20, 2.32, 2.47, 2.65, 2.87, 3.16, 3.53, 4.02, 4.69, 5.64, 7.07, 9.35, 13.34, 21.43]) x_uncertainty = x * uncertainty x = unumpy.uarray(x, x_uncertainty) y = np.array([0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0, 95.0, 100.0]) y = unumpy.uarray(y, y_error) n = np.arange(0, 5, 0.005) coefficient_determination_on = np.empty(shape = (len(n),)) for j in range(len(n)): n_correlation = n[j] x_fit = 1 / ((x) ** n_correlation) y_fit = y fit_a_raw, fit_b_raw = optimize.curve_fit(linear_fit, x_fit, y_fit)[0] x_prediction = (fit_a_raw / ((x) ** n_correlation)) + fit_b_raw y_residual_squares = np.sum((x_prediction - y) ** 2) y_total_squares = np.sum((y - np.mean(y)) ** 2) coefficient_determination_on[j] = 1 - (y_residual_squares / y_total_squares) 

Permítanme primero prefaciar esto, ya que este problema es imposible de resolver “a la perfección”, dado que usted quiere resolver para a , b n . Esto se debe a que para una n fija, su problema admite una solución de forma cerrada, mientras que si deja n libre, no lo hace, y de hecho, el problema puede tener varias soluciones. Por lo tanto, el análisis de error clásico (como el que utilizan las no uncertanities ) se rompe y hay que recurrir a otros métodos.

El caso n arreglado

Si n es fijo, su problema es que las bibliotecas a las que llama no son compatibles con uarray , por lo que debe hacer una solución. Afortunadamente, el ajuste lineal (bajo la distancia de l2) es simplemente mínimos cuadrados lineales que admite una solución de forma cerrada, que solo requiere rellenar los valores con unos y luego resolver las ecuaciones normales .

introduzca la descripción de la imagen aquí

Dónde:

introduzca la descripción de la imagen aquí

Podrías hacerlo así:

 import numpy as np from uncertainties import unumpy uncertainty = 0.5 y_error = 1.2 n = 1.0 # Define x and y x = np.array([1.71, 1.76, 1.81, 1.86, 1.93, 2.01, 2.09, 2.20, 2.32, 2.47, 2.65, 2.87, 3.16, 3.53, 4.02, 4.69, 5.64, 7.07, 9.35, 13.34, 21.43]) # Take power of x values according to n x_pow = x ** n x_uncertainty = x_pow * uncertainty x_fit = unumpy.uarray(np.c_[x_pow, np.ones_like(x)], np.c_[x_uncertainty, np.zeros_like(x_uncertainty)]) y = np.array([0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0, 95.0, 100.0]) y_fit = unumpy.uarray(y, y_error) # Use normal equations to find coefficients inv_mat = unumpy.ulinalg.pinv(x_fit.T.dot(x_fit)) fit_a, fit_b = inv_mat.dot(x_fit.T.dot(y_fit)) print('fit_a={}, fit_b={}'.format(fit_a, fit_b)) 

Resultado:

 fit_a=4.8+/-2.6, fit_b=28+/-10 

El caso n desconocido

Con n desconocido, realmente tiene algún problema ya que el problema no es convexo. Aquí, el análisis de error lineal (como lo hacen las uncertainties ) no funcionará.

Una solución es realizar una inferencia bayesiana , usando algún paquete como pymc . Si estás interesado en esto, podría intentar hacer un informe, pero no sería tan limpio como el anterior.

Siguiendo un poco el caso de una función lineal , creo que podría hacerse de manera similar. La solución para el lagrangiano, sin embargo, parece ser muy tediosa, aunque es posible, por supuesto. Un inventado hizo una medida diferente que parece plausible y debería dar un resultado muy similar. Tomando la elipse de error, la vuelvo a escalar de modo que la gráfica se convierta en tangente. Tomo la distancia hasta ese punto de contacto ( X_k,Y_k ) como medida para chi-cuadrado, que se calcula a partir de (x_k-X_k/sx_k)**2+(y_k-Y_k/sy_k)**2 . Es plausible, ya que, en el caso de los errores puros, este es el ajuste por mínimos cuadrados estándar. Para puros errores de x simplemente cambia. Para iguales errores de x,y daría la regla perpendicular, es decir, la distancia más corta. Con la función chi-cuadrado correspondiente, scipy.optimize.leastsq ya proporciona la matriz de covarianza aproximada de Hesse. En detalle uno tiene que escalarlo , sin embargo. También tenga en cuenta que los parámetros están fuertemente correlacionados.

Mi procedimiento es el siguiente:

 import matplotlib matplotlib.use('Qt5Agg') import matplotlib.pyplot as plt import numpy as np import myModules.multipoleMoments as mm from random import random from scipy.optimize import minimize,leastsq ###for gaussion distributed errors def boxmuller(x0,sigma): u1=random() u2=random() ll=np.sqrt(-2*np.log(u1)) z0=ll*np.cos(2*np.pi*u2) z1=ll*np.cos(2*np.pi*u2) return sigma*z0+x0, sigma*z1+x0 ###function to fit def f(t,c,d,n): return c+d*np.abs(t)**n ###to create some test data def test_data(c,d,n, xList,const_sx,rel_sx,const_sy,rel_sy): yList=[f(t,c,d,n) for t in xList] xErrList=[ boxmuller(x,const_sx+x*rel_sx)[0] for x in xList] yErrList=[ boxmuller(y,const_sy+y*rel_sy)[0] for y in yList] return xErrList,yErrList ###how to rescale the ellipse to make fitfunction a tangent def elliptic_rescale(x,c,d,n,x0,y0,sa,sb): y=f(x,c,d,n) r=np.sqrt((x-x0)**2+(y-y0)**2) kappa=float(sa)/float(sb) tau=np.arctan2(y-y0,x-x0) new_a=r*np.sqrt(np.cos(tau)**2+(kappa*np.sin(tau))**2) return new_a ###for plotting ellipses def ell_data(a,b,x0=0,y0=0): tList=np.linspace(0,2*np.pi,150) k=float(a)/float(b) rList=[a/np.sqrt((np.cos(t))**2+(k*np.sin(t))**2) for t in tList] xyList=np.array([[x0+r*np.cos(t),y0+r*np.sin(t)] for t,r in zip(tList,rList)]) return xyList ###residual function to calculate chi-square def residuals(parameters,dataPoint):#data point is (x,y,sx,sy) c,d,n= parameters theData=np.array(dataPoint) best_t_List=[] for i in range(len(dataPoint)): x,y,sx,sy=dataPoint[i][0],dataPoint[i][1],dataPoint[i][2],dataPoint[i][3] ###getthe point on the graph where it is tangent to an error-ellipse ed_fit=minimize(elliptic_rescale,0,args=(c,d,n,x,y,sx,sy) ) best_t=ed_fit['x'][0] best_t_List+=[best_t] best_y_List=[f(t,c,d,n) for t in best_t_List] ##weighted distance not squared yet, as this is done by scipy.optimize.leastsq wighted_dx_List=[(x_b-x_f)/sx for x_b,x_f,sx in zip( best_t_List,theData[:,0], theData[:,2] ) ] wighted_dy_List=[(x_b-x_f)/sx for x_b,x_f,sx in zip( best_y_List,theData[:,1], theData[:,3] ) ] return wighted_dx_List+wighted_dy_List ###some start values cc,dd,nn=2.177,.735,1.75 ssaa,ssbb=1,3 xx0,yy0=2,3 csx,rsx=.1,.05 csy,rsy=.4,.00 ####some data xThData=np.linspace(0,3,15) yThData=[ f(t, cc,dd,nn) for t in xThData] ###some noisy data xNoiseData,yNoiseData=test_data(cc,dd,nn, xThData, csx,rsx, csy,rsy) xGuessdError=[csx+rsx*x for x in xNoiseData] yGuessdError=[csy+rsy*y for y in yNoiseData] ###plot that fig1 = plt.figure(1) ax=fig1.add_subplot(1,1,1) ax.plot(xThData,yThData) ax.errorbar(xNoiseData,yNoiseData, xerr=xGuessdError, yerr=yGuessdError, fmt='none',ecolor='r') #### and plot... for i in range(len(xNoiseData)): ###...an elliple on the error values el0=ell_data(xGuessdError[i],yGuessdError[i],x0=xNoiseData[i],y0=yNoiseData[i]) ax.plot(*zip(*el0),linewidth=1, color="#808080",linestyle='-') ####...as well as a scaled version that touches the original graph. This gives the error shortest distance to that graph ed_fit=minimize(elliptic_rescale,0,args=(cc,dd,nn,xNoiseData[i],yNoiseData[i],xGuessdError[i],yGuessdError[i]) ) best_t=ed_fit['x'][0] best_a=elliptic_rescale(best_t,cc,dd,nn,xNoiseData[i],yNoiseData[i],xGuessdError[i],yGuessdError[i]) el0=ell_data(best_a,best_a/xGuessdError[i]*yGuessdError[i],x0=xNoiseData[i],y0=yNoiseData[i]) ax.plot(*zip(*el0),linewidth=1, color="#a000a0",linestyle='-') ###Now fitting zipData=zip(xNoiseData,yNoiseData, xGuessdError, yGuessdError) estimate = [2.0,1,2.5] bestFitValues, cov,info,mesg, ier = leastsq(residuals, estimate,args=(zipData), full_output=1) print bestFitValues ####covariance matrix ####note eg: https://stackoverflow.com/questions/14854339/in-scipy-how-and-why-does-curve-fit-calculate-the-covariance-of-the-parameter-es s_sq = (np.array(residuals(bestFitValues, zipData))**2).sum()/(len(zipData)-len(estimate)) pcov = cov * s_sq print pcov #### and plot the result... ax.plot(xThData,[f(x,*bestFitValues) for x in xThData]) for i in range(len(xNoiseData)): ####...as well as a scaled ellipses that touches the fitted graph. ed_fit=minimize(elliptic_rescale,0,args=(bestFitValues[0],bestFitValues[1],bestFitValues[2],xNoiseData[i],yNoiseData[i],xGuessdError[i],yGuessdError[i]) ) best_t=ed_fit['x'][0] best_a=elliptic_rescale(best_t,bestFitValues[0],bestFitValues[1],bestFitValues[2],xNoiseData[i],yNoiseData[i],xGuessdError[i],yGuessdError[i]) el0=ell_data(best_a,best_a/xGuessdError[i]*yGuessdError[i],x0=xNoiseData[i],y0=yNoiseData[i]) ax.plot(*zip(*el0),linewidth=1, color="#f0a000",linestyle='-') #~ ax.grid(None) plt.show() 

introduzca la descripción de la imagen aquí El gráfico azul es la función original. Los puntos de datos con barras de error rojas se calculan a partir de esto. Una elipse gris muestra la “línea de densidad de probabilidad constante”. Las elipses púrpuras tienen el gráfico original como tangente, las elipses naranjas tienen el ajuste como tangente

Aquí los mejores valores de ajuste son (no sus datos):

 [ 2.16146783 0.80204967 1.69951865] 

La matriz de la covarianza tiene la forma:

 [[ 0.0644794 -0.05418743 0.05454876] [-0.05418743 0.07228771 -0.08172885] [ 0.05454876 -0.08172885 0.10173394]] 

Editar Pensando en la “distancia elíptica”, creo que esto es exactamente lo que hace el enfoque lagrangiano en el documento vinculado. Solo que para una línea recta puede anotar una solución exacta, mientras que en este caso no puede.

Actualización tuve algunos problemas con los datos del OP. Aunque funciona con reescalado. Como la pendiente y el exponente están altamente correlacionados, primero tuve que descubrir cómo cambia la matriz de covarianza para los datos reescalados. Los detalles se pueden encontrar en J. Phys. Chem. 105 (2001) 3917

Usando las definiciones de funciones de arriba, el tratamiento de datos se vería así:

 ###some start values cc,dd,nn=.2,1,7.5 csx,rsx=2.0,.0 csy,rsy=0,.5 ###some noisy data yNoiseData=np.array([1.71,1.76, 1.81, 1.86, 1.93, 2.01, 2.09, 2.20, 2.32, 2.47, 2.65, 2.87, 3.16, 3.53, 4.02, 4.69, 5.64, 7.07, 9.35,13.34,21.43]) xNoiseData=np.array([0.0,5.0,10.0,15.0,20.0,25.0,30.0,35.0,40.0,45.0,50.0,55.0,60.0,65.0,70.0,75.0,80.0,85.0,90.0,95.0,100.0]) xGuessdError=[csx+rsx*x for x in xNoiseData] yGuessdError=[csy+rsy*y for y in yNoiseData] ###plot that fig1 = plt.figure(1) ax=fig1.add_subplot(1,2,2) bx=fig1.add_subplot(1,2,1) ax.errorbar(xNoiseData,yNoiseData, xerr=xGuessdError, yerr=yGuessdError, fmt='none',ecolor='r') ####rescaling print "\n++++++++++++++++++++++++ scaled ++++++++++++++++++++++++\n" alpha=.05 beta=.01 xNoiseDataS = [ beta*x for x in xNoiseData ] yNoiseDataS = [ alpha*x for x in yNoiseData ] xGuessdErrorS = [ beta*x for x in xGuessdError ] yGuessdErrorS = [ alpha*x for x in yGuessdError ] xtmp=np.linspace(0,1.1,25) bx.errorbar(xNoiseDataS,yNoiseDataS, xerr=xGuessdErrorS, yerr=yGuessdErrorS, fmt='none',ecolor='r') ###Now fitting zipData=zip(xNoiseDataS,yNoiseDataS, xGuessdErrorS, yGuessdErrorS) estimate = [.1,1,7.5] bestFitValues, cov,info,mesg, ier = leastsq(residuals, estimate,args=(zipData), full_output=1) print bestFitValues plt.plot(xtmp,[ f(x,*bestFitValues)for x in xtmp]) ####covariance matrix s_sq = (np.array(residuals(bestFitValues, zipData))**2).sum()/(len(zipData)-len(estimate)) pcov = cov * s_sq print pcov #### scale back print "\n++++++++++++++++++++++++ scaled back ++++++++++++++++++++++++\n" realBestFitValues= [bestFitValues[0]/alpha, bestFitValues[1]/alpha*(beta)**bestFitValues[2],bestFitValues[2] ] print realBestFitValues uMX = np.array( [[1/alpha,0,0],[0,beta**bestFitValues[2]/alpha,bestFitValues[1]/alpha*beta**bestFitValues[2]*np.log(beta)],[0,0,1]] ) uMX_T = uMX.transpose() realCov = np.dot(uMX, np.dot(pcov,uMX_T)) print realCov for i,para in enumerate(["b","a","n"]): print para+" = "+"{:.2e}".format(realBestFitValues[i])+" +/- "+"{:.2e}".format(np.sqrt(realCov[i,i])) ax.plot(xNoiseData,[f(x,*realBestFitValues) for x in xNoiseData]) plt.show() 

Datos de OP Así que los datos están razonablemente ajustados. Creo que, sin embargo, también hay un término lineal.

La salida proporciona:

 ++++++++++++++++++++++++ scaled ++++++++++++++++++++++++ [ 0.09788886 0.69614911 5.2221032 ] [[ 1.25914194e-05 2.86541696e-05 6.03957467e-04] [ 2.86541696e-05 3.88675025e-03 2.00199108e-02] [ 6.03957467e-04 2.00199108e-02 1.75756532e-01]] ++++++++++++++++++++++++ scaled back ++++++++++++++++++++++++ [1.9577772055183396, 5.0064036934715239e-10, 5.2221031993990517] [[ 5.03656777e-03 -2.74367539e-11 1.20791493e-02] [ -2.74367539e-11 8.69854174e-19 -3.90815222e-10] [ 1.20791493e-02 -3.90815222e-10 1.75756532e-01]] b = 1.96e+00 +/- 7.10e-02 a = 5.01e-10 +/- 9.33e-10 n = 5.22e+00 +/- 4.19e-01 

Uno puede ver una fuerte correlación de pendiente y exponente en la matriz de covarianza. También tenga en cuenta que el error de la pendiente es enorme.

BTW utilizando como modelo b+a*x**n + e*x obtengo con término lineal adicional

 ++++++++++++++++++++++++ scaled ++++++++++++++++++++++++ [ 0.08050174 0.78438855 8.11845402 0.09581568] [[ 5.96210962e-06 3.07651631e-08 -3.57876577e-04 -1.75228231e-05] [ 3.07651631e-08 1.39368435e-03 9.85025139e-03 1.83780053e-05] [ -3.57876577e-04 9.85025139e-03 1.85226736e-01 2.26973118e-03] [ -1.75228231e-05 1.83780053e-05 2.26973118e-03 7.92853339e-05]] ++++++++++++++++++++++++ scaled back ++++++++++++++++++++++++ [1.6100348667765145, 9.0918698097511416e-16, 8.1184540175879985, 0.019163135651422442] [[ 2.38484385e-03 2.99690170e-17 -7.15753154e-03 -7.00912926e-05] [ 2.99690170e-17 3.15340690e-30 -7.64119623e-16 -1.89639468e-18] [ -7.15753154e-03 -7.64119623e-16 1.85226736e-01 4.53946235e-04] [ -7.00912926e-05 -1.89639468e-18 4.53946235e-04 3.17141336e-06]] b = 1.61e+00 +/- 4.88e-02 a = 9.09e-16 +/- 1.78e-15 n = 8.12e+00 +/- 4.30e-01 e = 1.92e-02 +/- 1.78e-03 

Claro, los ajustes siempre se mejoran si se agregan parámetros, pero creo que esto se ve bastante razonable aquí (incluso podría ser que es b + a*x**n+e*x**m , pero esto va demasiado lejos) .)