python - lorentzian - scipy curve fit chi square

¿Por qué scipy.optimize.curve_fit no se ajusta a los datos?

Question

He estado tratando de ajustar un exponencial a algunos datos durante un tiempo usando scipy.optimize.curve_fit pero tengo una dificultad real. Realmente no puedo ver ninguna razón por la que esto no funcionaría, pero solo produce una línea recta, ¡ni idea de por qué!

Cualquier ayuda sería muy apreciada.

from __future__ import division
import numpy
from scipy.optimize import curve_fit
import matplotlib.pyplot as pyplot
def func(x,a,b,c):
return a*numpy.exp(-b*x)-c
yData = numpy.load(''yData.npy'')
xData = numpy.load(''xData.npy'')
trialX = numpy.linspace(xData[0],xData[-1],1000)
# Fit a polynomial
fitted = numpy.polyfit(xData, yData, 10)[::-1]
y = numpy.zeros(len(trailX))
for i in range(len(fitted)):
y += fitted[i]*trialX**i
# Fit an exponential
popt, pcov = curve_fit(func, xData, yData)
yEXP = func(trialX, *popt)
pyplot.figure()
pyplot.plot(xData, yData, label=''Data'', marker=''o'')
pyplot.plot(trialX, yEXP, ''r-'',ls=''--'', label="Exp Fit")
pyplot.plot(trialX, y, label = ''10 Deg Poly'')
pyplot.legend()
pyplot.show()

xData = [1e-06, 2e-06, 3e-06, 4e-06,
5e-06, 6e-06, 7e-06, 8e-06,
9e-06, 1e-05, 2e-05, 3e-05,
4e-05, 5e-05, 6e-05, 7e-05,
8e-05, 9e-05, 0.0001, 0.0002,
0.0003, 0.0004, 0.0005, 0.0006,
0.0007, 0.0008, 0.0009, 0.001,
0.002, 0.003, 0.004, 0.005,
0.006, 0.007, 0.008, 0.009, 0.01]
yData =
[6.37420666067e-09, 1.13082012115e-08,
1.52835756975e-08, 2.19214493931e-08, 2.71258852882e-08, 3.38556130078e-08, 3.55765277358e-08,
4.13818145846e-08, 4.72543475372e-08, 4.85834751151e-08, 9.53876562077e-08, 1.45110636413e-07,
1.83066627931e-07, 2.10138415308e-07, 2.43503982686e-07, 2.72107045549e-07, 3.02911771395e-07,
3.26499455951e-07, 3.48319349445e-07, 5.13187669283e-07, 5.98480176303e-07, 6.57028222701e-07,
6.98347073045e-07, 7.28699930335e-07, 7.50686502279e-07, 7.7015576866e-07, 7.87147246927e-07,
7.99607141001e-07, 8.61398763228e-07, 8.84272900407e-07, 8.96463883243e-07, 9.04105135329e-07,
9.08443443149e-07, 9.12391264185e-07, 9.150842683e-07, 9.16878548643e-07, 9.18389990067e-07]

Answer 1

Los algoritmos numéricos tienden a funcionar mejor cuando no se alimentan con números extremadamente pequeños (o grandes).

En este caso, el gráfico muestra que sus datos tienen valores x e y extremadamente pequeños. Si los escala, el ajuste es notable mejor:

xData = np.load(''xData.npy'')*10**5
yData = np.load(''yData.npy'')*10**5

from __future__ import division
import os
os.chdir(os.path.expanduser(''~/tmp''))
import numpy as np
import scipy.optimize as optimize
import matplotlib.pyplot as plt
def func(x,a,b,c):
return a*np.exp(-b*x)-c
xData = np.load(''xData.npy'')*10**5
yData = np.load(''yData.npy'')*10**5
print(xData.min(), xData.max())
print(yData.min(), yData.max())
trialX = np.linspace(xData[0], xData[-1], 1000)
# Fit a polynomial
fitted = np.polyfit(xData, yData, 10)[::-1]
y = np.zeros(len(trialX))
for i in range(len(fitted)):
y += fitted[i]*trialX**i
# Fit an exponential
popt, pcov = optimize.curve_fit(func, xData, yData)
print(popt)
yEXP = func(trialX, *popt)
plt.figure()
plt.plot(xData, yData, label=''Data'', marker=''o'')
plt.plot(trialX, yEXP, ''r-'',ls=''--'', label="Exp Fit")
plt.plot(trialX, y, label = ''10 Deg Poly'')
plt.legend()
plt.show()

Tenga en cuenta que después de xData yData xData y yData , los parámetros devueltos por curve_fit también deben ser escalados. En este caso, a , b deben dividirse por 10 ** 5 para obtener los parámetros ajustados para los datos originales.

Una objeción que podría tener a lo anterior es que la escala debe elegirse en lugar de "cuidadosamente". (Leer: ¡No todas las elecciones razonables de escala funcionan!)

Puede mejorar la robustez de curve_fit proporcionando una estimación inicial razonable para los parámetros. Por lo general, usted tiene algún conocimiento a priori sobre los datos que pueden motivar las conjeturas del tipo de campo / fondo de juego para valores de parámetros razonables.

Por ejemplo, llamando a curve_fit con

guess = (-1, 0.1, 0)
popt, pcov = optimize.curve_fit(func, xData, yData, guess)

ayuda a mejorar el rango de escalas en las que curve_fit tiene éxito en este caso.

Answer 2

Una mejora (leve) de esta solución, sin tener en cuenta el conocimiento a priori de los datos, podría ser la siguiente: tome la media inversa del conjunto de datos y utilícelo como el "factor de escala" para pasarlo a los lesssq subyacentes () llamado por curve_fit (). Esto permite que el instalador funcione y devuelve los parámetros en la escala original de los datos.

La línea relevante es:

popt, pcov = curve_fit(func, xData, yData)

que se convierte en:

popt, pcov = curve_fit(func, xData, yData,
diag=(1./xData.mean(),1./yData.mean()) )

Aquí está el ejemplo completo que produce esta imagen:

from __future__ import division
import numpy
from scipy.optimize import curve_fit
import matplotlib.pyplot as pyplot
def func(x,a,b,c):
return a*numpy.exp(-b*x)-c
xData = numpy.array([1e-06, 2e-06, 3e-06, 4e-06, 5e-06, 6e-06,
7e-06, 8e-06, 9e-06, 1e-05, 2e-05, 3e-05, 4e-05, 5e-05, 6e-05,
7e-05, 8e-05, 9e-05, 0.0001, 0.0002, 0.0003, 0.0004, 0.0005,
0.0006, 0.0007, 0.0008, 0.0009, 0.001, 0.002, 0.003, 0.004, 0.005
, 0.006, 0.007, 0.008, 0.009, 0.01])
yData = numpy.array([6.37420666067e-09, 1.13082012115e-08,
1.52835756975e-08, 2.19214493931e-08, 2.71258852882e-08,
3.38556130078e-08, 3.55765277358e-08, 4.13818145846e-08,
4.72543475372e-08, 4.85834751151e-08, 9.53876562077e-08,
1.45110636413e-07, 1.83066627931e-07, 2.10138415308e-07,
2.43503982686e-07, 2.72107045549e-07, 3.02911771395e-07,
3.26499455951e-07, 3.48319349445e-07, 5.13187669283e-07,
5.98480176303e-07, 6.57028222701e-07, 6.98347073045e-07,
7.28699930335e-07, 7.50686502279e-07, 7.7015576866e-07,
7.87147246927e-07, 7.99607141001e-07, 8.61398763228e-07,
8.84272900407e-07, 8.96463883243e-07, 9.04105135329e-07,
9.08443443149e-07, 9.12391264185e-07, 9.150842683e-07,
9.16878548643e-07, 9.18389990067e-07])
trialX = numpy.linspace(xData[0],xData[-1],1000)
# Fit a polynomial
fitted = numpy.polyfit(xData, yData, 10)[::-1]
y = numpy.zeros(len(trialX))
for i in range(len(fitted)):
y += fitted[i]*trialX**i
# Fit an exponential
popt, pcov = curve_fit(func, xData, yData,
diag=(1./xData.mean(),1./yData.mean()) )
yEXP = func(trialX, *popt)
pyplot.figure()
pyplot.plot(xData, yData, label=''Data'', marker=''o'')
pyplot.plot(trialX, yEXP, ''r-'',ls=''--'', label="Exp Fit")
pyplot.plot(trialX, y, label = ''10 Deg Poly'')
pyplot.legend()
pyplot.show()

Answer 3

el modelo a*exp(-b*x)+c encaja bien con los datos, pero sugiero una pequeña modificación:
usa esto en su lugar

a*x*exp(-b*x)+c

buena suerte