Curva exponencial ajustada en ciencia ficción

Tengo dos matrices NumPy x y y . Cuando bash ajustar mis datos utilizando la función exponencial y curve_fit (SciPy) con este código simple

 #!/usr/bin/env python from pylab import * from scipy.optimize import curve_fit x = np.array([399.75, 989.25, 1578.75, 2168.25, 2757.75, 3347.25, 3936.75, 4526.25, 5115.75, 5705.25]) y = np.array([109,62,39,13,10,4,2,0,1,2]) def func(x, a, b, c, d): return a*np.exp(bc*x)+d popt, pcov = curve_fit(func, x, y) 

Me sale mal coeficientes popt

 [a,b,c,d] = [1., 1., 1., 24.19999988] 

¿Cuál es el problema?

Primer comentario: dado que a*exp(b - c*x) = (a*exp(b))*exp(-c*x) = A*exp(-c*x) , a o b es redundante. Dejaré caer y usaré:

 def func(x, a, c, d): return a*np.exp(-c*x)+d 

Ese no es el problema principal. El problema es simplemente que curve_fit no puede converger en una solución a este problema cuando se usa la estimación inicial predeterminada (que es todo 1s). Compruebe pcov ; Verás que es inf . Esto no es sorprendente, porque si c es 1, la mayoría de los valores de exp(-c*x) desbordan a 0:

 In [32]: np.exp(-x) Out[32]: array([ 2.45912644e-174, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, 0.00000000e+000]) 

Esto sugiere que c debería ser pequeño. Una mejor estimación inicial es, digamos, p0 = (1, 1e-6, 1) . Entonces me sale:

 In [36]: popt, pcov = curve_fit(func, x, y, p0=(1, 1e-6, 1)) In [37]: popt Out[37]: array([ 1.63561656e+02, 9.71142196e-04, -1.16854450e+00]) 

Esto parece razonable:

 In [42]: xx = np.linspace(300, 6000, 1000) In [43]: yy = func(xx, *popt) In [44]: plot(x, y, 'ko') Out[44]: [] In [45]: plot(xx, yy) Out[45]: [] 

parcela de ajuste

En primer lugar, recomendaría modificar su ecuación a a*np.exp(-c*(xb))+d , de lo contrario, la exponencial siempre estará centrada en x=0 que no siempre es el caso. También debe especificar condiciones iniciales razonables (el cuarto argumento de curve_fit especifica condiciones iniciales para [a,b,c,d] ).

Este código encaja muy bien:

 from pylab import * from scipy.optimize import curve_fit x = np.array([399.75, 989.25, 1578.75, 2168.25, 2757.75, 3347.25, 3936.75, 4526.25, 5115.75, 5705.25]) y = np.array([109,62,39,13,10,4,2,0,1,2]) def func(x, a, b, c, d): return a*np.exp(-c*(xb))+d popt, pcov = curve_fit(func, x, y, [100,400,0.001,0]) print popt plot(x,y) x=linspace(400,6000,10000) plot(x,func(x,*popt)) show()