conexión de mínimos cuadrados no lineales python

Question

Estoy un poco fuera de mi profundidad en términos de las matemáticas involucradas en mi problema, así que me disculpo por cualquier nomenclatura incorrecta.

Estaba viendo el uso de la función scraps, leastsq, pero no estoy seguro si es la función correcta. tengo la siguiente ecuación:

eq = lambda PLP,p0,l0,kd : 0.5*(-1-((p0+l0)/kd) + np.sqrt(4*(l0/kd)+(((l0-p0)/kd)-1)**2)) 

tengo de datos (8 juegos) para todos los términos excepto para kd (PLP, p0, L0). Necesito encontrar el valor de kd por regresión no lineal de la ecuación anterior. De los ejemplos que he leído, leastsq parece no permitir la entrada de datos, para obtener la salida que necesito.

Gracias por su ayuda

Answer 1

Este es un ejemplo muy básicos de cómo utilizar scipy.optimize.leastsq:

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

def func(kd,p0,l0): 
    return 0.5*(-1-((p0+l0)/kd) + np.sqrt(4*(l0/kd)+(((l0-p0)/kd)-1)**2)) 

La suma de los cuadrados de la residuals es la función de kd que estamos tratando para minimizar:

def residuals(kd,p0,l0,PLP): 
    return PLP - func(kd,p0,l0) 

Aquí genero algunos datos aleatorios. En su lugar, querrá cargar sus datos reales.

N=1000 
kd_guess=3.5 # <-- You have to supply a guess for kd 
p0 = np.linspace(0,10,N) 
l0 = np.linspace(0,10,N) 
PLP = func(kd_guess,p0,l0)+(np.random.random(N)-0.5)*0.1 

kd,cov,infodict,mesg,ier = optimize.leastsq(
    residuals,kd_guess,args=(p0,l0,PLP),full_output=True,warning=True) 

print(kd) 

produce algo así como

3.49914274899 

Este es el mejor valor de ajuste para kd encontrado por optimize.leastsq.

Aquí generar el valor de PLP utilizando el valor de kd que acabamos de encontrar:

PLP_fit=func(kd,p0,l0) 

A continuación se muestra un gráfico de PLP frente p0. La línea azul proviene de datos, la línea roja es la mejor curva de ajuste.

plt.plot(p0,PLP,'-b',p0,PLP_fit,'-r') 
plt.show() 

Answer 2

Otra opción es utilizar lmfit.

Proporcionan un gran example para empezar :.

#!/usr/bin/env python 
#<examples/doc_basic.py> 
from lmfit import minimize, Minimizer, Parameters, Parameter, report_fit 
import numpy as np 

# create data to be fitted 
x = np.linspace(0, 15, 301) 
data = (5. * np.sin(2 * x - 0.1) * np.exp(-x*x*0.025) + 
     np.random.normal(size=len(x), scale=0.2)) 

# define objective function: returns the array to be minimized 
def fcn2min(params, x, data): 
    """ model decaying sine wave, subtract data""" 
    amp = params['amp'] 
    shift = params['shift'] 
    omega = params['omega'] 
    decay = params['decay'] 
    model = amp * np.sin(x * omega + shift) * np.exp(-x*x*decay) 
    return model - data 

# create a set of Parameters 
params = Parameters() 
params.add('amp', value= 10, min=0) 
params.add('decay', value= 0.1) 
params.add('shift', value= 0.0, min=-np.pi/2., max=np.pi/2) 
params.add('omega', value= 3.0) 


# do fit, here with leastsq model 
minner = Minimizer(fcn2min, params, fcn_args=(x, data)) 
kws = {'options': {'maxiter':10}} 
result = minner.minimize() 


# calculate final result 
final = data + result.residual 

# write error report 
report_fit(result) 

# try to plot results 
try: 
    import pylab 
    pylab.plot(x, data, 'k+') 
    pylab.plot(x, final, 'r') 
    pylab.show() 
except: 
    pass 

#<end of examples/doc_basic.py>