蟒蛇非線性最小二乘擬合

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>