Scipy curve_fit error: se divide por cero

He estado tratando de ajustar una función a algunos datos durante un tiempo usando scipy.optimize.curve_fit:

from __future__ import (print_function, division, unicode_literals, absolute_import) import numpy as np from scipy.optimize import curve_fit import matplotlib.pyplot as mpl x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]) y = np.array([20.8, 20.9, 22.9, 25.2, 26.9, 28.3, 29.5, 30.7, 31.8, 32.9, 34.0, 35.3, 36.4, 37.5, 38.6, 39.6, 40.6, 41.6, 42.5, 43.2, 44.2, 45.0, 45.8, 46.5, 47.3, 48.0, 48.6, 49.2, 49.8, 50.4]) def f(x, a, b, c): return a/(1+b*x**c) popt, pcov = curve_fit(f, x, y) print(popt, np.sqrt(np.diag(pcov)), sep='\n') 

Pero siempre aparece un error:

 RuntimeWarning: divide by zero encountered in power return a/(1+b*x**c) 

Tal vez alguien me pueda ayudar a evitarlo? Cualquier ayuda sería muy apreciada. ¡Aclamaciones!

Muy bien, dos trucos útiles.

Primero, reemplace 0 en su x con un número realmente pequeño, como 1e-8 (no se ría, hay un paquete central en R que hace esto, escrito por his name shall not be spoken y la gente lo usará todo el tiempo) .) En realidad no conseguí su RuntimeWarning en absoluto. Estoy ejecutando scipy 0.12.0 y numpy 1.7.1 . Tal vez esto depende de la versión.

Pero tendremos un ajuste muy malo:

 In [41]: popt, pcov Out[41]: (array([ 3.90107143e+01, -3.08698757e+07, -1.52971609e+02]), inf) 

Entonces, el truco 2, en lugar de optimizar la función f , definimos una función g :

 In [38]: def g(x, a, b, c): ....: return b/a*x**c+1/a ....: In [39]: curve_fit(g, x, 1/y) #Better fit Out[39]: (array([ 19.76748582, -0.14499508, 0.44206688]), array([[ 0.29043958, 0.00899521, 0.01650935], [ 0.00899521, 0.00036082, 0.00070345], [ 0.01650935, 0.00070345, 0.00140253]])) 

Ahora podemos usar el vector de parámetros resultante como vector de inicio para optimizar f() . Como curve_fit es un método de mínimos cuadrados no lineales, el parámetro optimiza g() no es necesario, el parámetro optimiza f() , pero esperamos que esté cerca. Las matrices de covarianza son muy diferentes, por supuesto.

 In [78]: curve_fit(f, x, y, p0=curve_fit(g, x, 1/y)[0]) #Alternative Fit Out[78]: (array([ 18.0480446 , -0.22881647, 0.31200106]), array([[ 1.14928169, 0.03741604, 0.03897652], [ 0.03741604, 0.00128511, 0.00136315], [ 0.03897652, 0.00136315, 0.00145614]])) 

La comparación de los resultados:

introduzca la descripción de la imagen aquí

Ahora el resultado es bastante bueno.

Sus valores de x comienzan en 0. Si, por alguna razón, el parámetro c es negativo durante el cálculo, entonces evaluará 0 elevado a un exponente negativo, que es una división por cero: para p>0 tenemos

 0**(-p) = 1/(0**p) = 1/0