python nonlineare minimi quadrati che si adattano

Question

Sono un po 'fuori dalla mia profondità in termini di matematica coinvolta nel mio problema, quindi mi scuso per qualsiasi nomenclatura errata.

Stavo cercando di usare la funzione scipy leastsq, ma non sono sicuro che sia la funzione corretta. Ho la seguente equazione:

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

Ho dati (8 set) per tutti i termini tranne kd (PLP, p0, l0). Devo trovare il valore di kd mediante la regressione non lineare dell'equazione di cui sopra. Dagli esempi che ho letto, leastsq sembra non consentire l'inserimento dei dati, per ottenere l'output di cui ho bisogno.

Grazie per il vostro aiuto

Answer 1

Questo è un esempio ridotto all'osso di come utilizzare 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 somma dei quadrati della residuals è la funzione di kd stiamo cercando per minimizzare:

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

Qui elabro alcuni dati casuali. Vorresti caricare qui i tuoi dati reali.

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) 

cede qualcosa come

3.49914274899 

Questo è il miglior valore in forma per kd trovata da optimize.leastsq.

Qui generare il valore della PLP utilizzando il valore per kd abbiamo appena trovato:

PLP_fit=func(kd,p0,l0) 

Qui di seguito è un grafico della PLP contro p0. La linea blu proviene dai dati, la linea rossa è la curva più adatta.

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

Answer 2

Un'altra opzione è quella di utilizzare lmfit.

Forniscono un ottimo example per iniziare :.

#!/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>