python - power - scipy gaussian fit

SciPy Curve Fit falla Ley de Poder

Question

Entonces, intento ajustar un conjunto de datos con una ley de poder del siguiente tipo:

def f(x,N,a): # Power law fit
    if a >0:
        return N*x**(-a)
    else:
        return 10.**300

par,cov = scipy.optimize.curve_fit(f,data,time,array([10**(-7),1.2]))

donde la condición else es solo para forzar a ser positivo. Usar scipy.optimize.curve_fit produce un ajuste horrible (línea verde) , devolviendo valores de 1.2e + 04 y 1.9e0-7 para N y a, respectivamente, sin ninguna intersección con los datos. De los ajustes que he puesto manualmente, los valores deberían ubicarse alrededor de 1e-07 y 1.2 para N y a, respectivamente, aunque ponerlos en curve_fit como parámetros iniciales no cambia el resultado. Eliminar la condición para que un resultado sea positivo da como resultado un ajuste peor, ya que elige un negativo, lo que conduce a un ajuste con la pendiente del signo equivocada.

No puedo encontrar la forma de obtener un ajuste creíble, ni mucho menos confiable, de esta rutina, pero no puedo encontrar ninguna otra buena rutina de ajuste de curva Python. ¿Tengo que escribir mi propio algoritmo de mínimos cuadrados o hay algo que estoy haciendo mal aquí?

Answer 1

ACTUALIZAR

En la publicación original, mostré una solución que usa lmfit que permite asignar límites a tus parámetros. Comenzando con la versión 0.17, scipy también permite asignar límites a sus parámetros directamente (ver documentación ). Encuentre esta solución a continuación después de EDITAR, que con suerte puede servir como un ejemplo mínimo sobre cómo usar scipy''s curve_fit con límites de parámetros.

Publicación original

Como lo sugirió @Warren Weckesser, podría usar lmfit para realizar esta tarea, lo que le permite asignar límites a sus parámetros y evitar esta cláusula ''fea''.

Como no proporciona ningún dato, creé algunos que se muestran aquí:

Siguen la ley f(x) = 10.5 * x ** (-0.08)

Los ajusto, como sugiere @ roadrunner66, al transformar la ley de potencia en una función lineal:

y = N * x ** a
ln(y) = ln(N * x ** a)
ln(y) = a * ln(x) + ln(N)

Así que primero uso np.log en los datos originales y luego hago el ajuste. Cuando ahora uso lmfit, obtengo el siguiente resultado:

[[Variables]]
    lN:   2.35450302 +/- 0.019531 (0.83%) (init= 1.704748)
    a:   -0.08035342 +/- 0.005158 (6.42%) (init=-0.5)

Entonces a está bastante cerca del valor original y np.exp(2.35450302) da 10.53 que también está muy cerca del valor original.

La trama se ve de la siguiente manera; como puede ver, el ajuste describe muy bien los datos:

Aquí está el código completo con un par de comentarios en línea:

import numpy as np
import matplotlib.pyplot as plt
from lmfit import minimize, Parameters, Parameter, report_fit

# generate some data with noise
xData = np.linspace(0.01, 100., 50.)
aOrg = 0.08
Norg = 10.5
yData = Norg * xData ** (-aOrg) + np.random.normal(0, 0.5, len(xData))
plt.plot(xData, yData, ''bo'')
plt.show()

# transform data so that we can use a linear fit
lx = np.log(xData)
ly = np.log(yData)
plt.plot(lx, ly, ''bo'')
plt.show()

def decay(params, x, data):

    lN = params[''lN''].value
    a = params[''a''].value

    # our linear model
    model = a * x + lN
    return model - data # that''s what you want to minimize

# create a set of Parameters
params = Parameters()
params.add(''lN'', value=np.log(5.5), min=0.01, max=100)  # value is the initial value
params.add(''a'', value=-0.5, min=-1, max=-0.001)  # min, max define parameter bounds

# do fit, here with leastsq model
result = minimize(decay, params, args=(lx, ly))

# write error report
report_fit(params)

# plot data
xnew = np.linspace(0., 100., 5000.)
# plot the data
plt.plot(xData, yData, ''bo'')
plt.plot(xnew, np.exp(result.values[''lN'']) * xnew ** (result.values[''a'']), ''r'')
plt.show()

EDITAR

Asumiendo que tiene scipy 0.17 instalado, también puede hacer lo siguiente usando curve_fit . Lo muestro para su definición original de la ley de potencia (línea roja en el siguiente diagrama), así como para los datos logarítmicos (línea negra en el gráfico a continuación). Los datos se generan de la misma manera que arriba. La trama se ve de la siguiente manera:

Como puede ver, los datos se describen muy bien. Si imprime popt y popt_log , obtiene una array([ 10.47463426, 0.07914812]) y una array([ 2.35158653, -0.08045776]) , respectivamente (nota: para la letra uno, tendrá que tomar el exponencial del primer argumento - np.exp(popt_log[0]) = 10.502 que está cerca de los datos originales.

Aquí está el código completo:

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

# generate some data with noise
xData = np.linspace(0.01, 100., 50)
aOrg = 0.08
Norg = 10.5
yData = Norg * xData ** (-aOrg) + np.random.normal(0, 0.5, len(xData))

# get logarithmic data
lx = np.log(xData)
ly = np.log(yData)

def f(x, N, a):
    return N * x ** (-a)

def f_log(x, lN, a):
    return a * x + lN

# optimize using the appropriate bounds
popt, pcov = curve_fit(f, xData, yData, bounds=(0, [30., 20.]))
popt_log, pcov_log = curve_fit(f_log, lx, ly, bounds=([0, -10], [30., 20.]))

xnew = np.linspace(0.01, 100., 5000)

# plot the data
plt.plot(xData, yData, ''bo'')
plt.plot(xnew, f(xnew, *popt), ''r'')
plt.plot(xnew, f(xnew, np.exp(popt_log[0]), -popt_log[1]), ''k'')
plt.show()