蟒蛇非线性最小二乘拟合

Question

我在我的问题涉及的数学方面有点超出我的深度，所以我对任何不正确的术语表道歉。

我在看使用scipy函数leastsq，但不知道它是否是正确的函数。 我有以下方程：

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

我有除了KD（PLP，P0，L0）的所有的条款数据（8套）。我需要通过上述方程的非线性回归来找到kd的值。 从我读过的例子中，leastsq似乎不允许输入数据，以获得我需要的输出。

谢谢您的帮助

Answer 1

这是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)) 

的residuals的平方的总和为​​我们试图功能最小化：

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

这里我生成一些随机数据。您想要在这里加载您的真实数据。

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) 

产生类似

3.49914274899 

这是optimize.leastsq发现​​最合适的值。

在这里，我们生成使用值​​的的PLP价值，我们才发现：

PLP_fit=func(kd,p0,l0) 

下面是PLP与p0一个阴谋。蓝线来自数据，红线是最合适的曲线。

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

Answer 2

另一种选择是使用lmfit。

他们提供了一个伟大的example让你开始：。

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