Cómo incluir correctamente las incertidumbres en el ajuste con python

Question

Estoy tratando de encajar algunos puntos de datos con y incertidumbres en python. Los datos están etiquetados en python como x, y y yerr. Necesito hacer un ajuste lineal en esa información en la escala loglog. Como referencia, si los resultados aptos son adecuadamente, comparo los resultados de pitón con los de SciDAVis

que probé curve_fit con

def func(x, a, b): 
    return np.exp(a* np.log(x)+np.log(b)) 

popt, pcov = curve_fit(func, x, y,sigma=yerr) 

, así como kmpfit con

def funcL(p, x): 
    a,b = p 
    return (np.exp(a*np.log(x)+np.log(b))) 

def residualsL(p, data): 
    a,b=p 
    x, y, errorfit = data 
    return (y-funcL(p,x))/errorfit 

a0=1 
b0=0.1 
p0 = [a0,b0] 
fitterL = kmpfit.Fitter(residuals=residualsL, data=(x,y,yerr)) 
fitterL.parinfo = [{}, {}] 
fitterL.fit(params0=p0) 

y cuando Estoy tratando de ajustar los datos con uno de esos sin incertidumbres (es decir, establecer yerr = 1), todo funciona bien y los resultados son idénticos a los de scidavis. Pero si configuro el yerr a las incertidumbres del archivo de datos, obtengo algunos resultados inquietantes. En python obtengo, por ejemplo, a = 0,86 y en scidavis a = 0,14. Leí algo acerca de que los errores están incluidos como ponderaciones. ¿Debo cambiar algo para calcular el ajuste correctamente? ¿O qué estoy haciendo mal?

edición: aquí es un ejemplo de un archivo de datos (x, y, yerr)

3.942387e-02 1.987800e+00 5.513165e-01 
6.623142e-02 7.126161e+00 1.425232e+00 
9.348280e-02 1.238530e+01 1.536208e+00 
1.353088e-01 1.090471e+01 7.829126e-01 
2.028446e-01 1.023087e+01 3.839575e-01 
3.058446e-01 8.403626e+00 1.756866e-01 
4.584524e-01 7.345275e+00 8.442288e-02 
6.879677e-01 6.128521e+00 3.847194e-02 
1.032592e+00 5.359025e+00 1.837428e-02 
1.549152e+00 5.380514e+00 1.007010e-02 
2.323985e+00 6.404229e+00 6.534108e-03 
3.355974e+00 9.489101e+00 6.342546e-03 
4.384128e+00 1.497998e+01 2.273233e-02 

y el resultado:

in python: 
    without uncertainties: a=0.06216 +/- 0.00650 ; b=8.53594 +/- 1.13985 
    with uncertainties: a=0.86051 +/- 0.01640 ; b=3.38081 +/- 0.22667 
in scidavis: 
    without uncertainties: a = 0.06216 +/- 0.08060; b = 8.53594 +/- 1.06763 
    with uncertainties: a = 0.14154 +/- 0.005731; b = 7.38213 +/- 2.13653 

Answer 1

Debo estar mal entendido algo. Sus datos publicado no se parece en nada

f(x,a,b) = np.exp(a*np.log(x)+np.log(b)) 

La línea roja es el resultado de scipy.optimize.curve_fit, la línea verde es el resultado de SciDAVis.

Supongo que ninguno de los algoritmos está convergiendo hacia un buen ajuste, por lo que no es sorprendente que los resultados no coincidan.

---

No puedo explicar cómo SciDAVis encuentra sus parámetros, pero de acuerdo con las definiciones como yo a entender, scipy es encontrar parámetros con menores residuos mínimos cuadrados de scidavis:

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

def func(x, a, b): 
    return np.exp(a* np.log(x)+np.log(b)) 

def sum_square(residuals): 
    return (residuals**2).sum() 

def residuals(p, x, y, sigma): 
    return 1.0/sigma*(y - func(x, *p)) 

data = np.loadtxt('test.dat').reshape((-1,3)) 
x, y, yerr = np.rollaxis(data, axis = 1) 
sigma = yerr 

popt, pcov = optimize.curve_fit(func, x, y, sigma = sigma, maxfev = 10000) 
print('popt: {p}'.format(p = popt)) 
scidavis = (0.14154, 7.38213) 
print('scidavis: {p}'.format(p = scidavis)) 

print('''\ 
sum of squares for scipy: {sp} 
sum of squares for scidavis: {d} 
'''.format(
      sp = sum_square(residuals(popt, x = x, y = y, sigma = sigma)), 
      d = sum_square(residuals(scidavis, x = x, y = y, sigma = sigma)) 
    )) 

plt.plot(x, y, 'bo', x, func(x,*popt), 'r-', x, func(x, *scidavis), 'g-') 
plt.errorbar(x, y, yerr) 
plt.show() 

rendimientos

popt: [ 0.86051258 3.38081125] 
scidavis: (0.14154, 7.38213) 
sum of squares for scipy: 53249.9915654 
sum of squares for scidavis: 239654.84276