Scipy的最佳曲线拟合极限
问题内容:
有什么方法可以为Scipy的“优化曲线拟合”提供限制吗?
我的例子:
def optimized_formula(x, m_1, m_2, y_1, y_2, ratio_2):
return (log(x[0]) * m_1 + m_2)*((1 - x[1]/max_age)*(1-ratio_2)) + ((log(x[1]) * y_1 + y_2)*(x[1]/max_age)*ratio_2)
popt, pcov = optimize.curve_fit(optimized_formula, usage_and_age, prices)
x [0]是age,max_age是常数。考虑到这一点,当x [0]接近最大值时,x [1] / max_age接近1。
是否可以提供一个约束/限制,从而使x [1] / max_age> 0.3和x [1] / max_age <0.7以及其他约束,例如m_1 <0,m_2>
0,依此类推。
问题答案:
如另一个答案中所建议,您可以将lmfit用于此类问题。因此,我还添加了一个有关如何使用它的示例,以防有人对此主题感兴趣。
假设您有一个数据集,如下所示:
xdata = np.array([177.,180.,183.,187.,189.,190.,196.,197.,201.,202.,203.,204.,206.,218.,225.,231.,234.,
252.,262.,266.,267.,268.,277.,286.,303.])
ydata = np.array([0.81,0.74,0.78,0.75,0.77,0.81,0.73,0.76,0.71,0.74,0.81,0.71,0.74,0.71,
0.72,0.69,0.75,0.59,0.61,0.63,0.64,0.63,0.35,0.27,0.26])
并且您想要使模型适合数据,如下所示:
model = n1 + (n2 * x + n3) * 1./ (1. + np.exp(n4 * (n5 - x)))
受到限制
0.2 < n1 < 0.8
-0.3 < n2 < 0
lmfit
然后使用(版本0.8.3)获得以下输出:
n1: 0.26564921 +/- 0.024765 (9.32%) (init= 0.2)
n2: -0.00195398 +/- 0.000311 (15.93%) (init=-0.005)
n3: 0.87261892 +/- 0.068601 (7.86%) (init= 1.0766)
n4: -1.43507072 +/- 1.223086 (85.23%) (init=-0.36379)
n5: 277.684530 +/- 3.768676 (1.36%) (init= 274)
如您所见,拟合非常好地再现了数据,并且参数在请求的范围内。
这是完整的代码,它重现了该图并带有一些其他注释:
from lmfit import minimize, Parameters, Parameter, report_fit
import numpy as np
xdata = np.array([177.,180.,183.,187.,189.,190.,196.,197.,201.,202.,203.,204.,206.,218.,225.,231.,234.,
252.,262.,266.,267.,268.,277.,286.,303.])
ydata = np.array([0.81,0.74,0.78,0.75,0.77,0.81,0.73,0.76,0.71,0.74,0.81,0.71,0.74,0.71,
0.72,0.69,0.75,0.59,0.61,0.63,0.64,0.63,0.35,0.27,0.26])
def fit_fc(params, x, data):
n1 = params['n1'].value
n2 = params['n2'].value
n3 = params['n3'].value
n4 = params['n4'].value
n5 = params['n5'].value
model = n1 + (n2 * x + n3) * 1./ (1. + np.exp(n4 * (n5 - x)))
return model - data #that's what you want to minimize
# create a set of Parameters
# 'value' is the initial condition
# 'min' and 'max' define your boundaries
params = Parameters()
params.add('n1', value= 0.2, min=0.2, max=0.8)
params.add('n2', value= -0.005, min=-0.3, max=10**(-10))
params.add('n3', value= 1.0766, min=-1000., max=1000.)
params.add('n4', value= -0.36379, min=-1000., max=1000.)
params.add('n5', value= 274.0, min=0., max=1000.)
# do fit, here with leastsq model
result = minimize(fit_fc, params, args=(xdata, ydata))
# write error report
report_fit(params)
xplot = np.linspace(min(xdata), max(xdata), 1000)
yplot = result.values['n1'] + (result.values['n2'] * xplot + result.values['n3']) * \
1./ (1. + np.exp(result.values['n4'] * (result.values['n5'] - xplot)))
#plot results
try:
import pylab
pylab.plot(xdata, ydata, 'k+')
pylab.plot(xplot, yplot, 'r')
pylab.show()
except:
pass
编辑:
如果您使用的是0.9.x版,则需要相应地调整代码。在此处检查从0.8.3到0.9.x进行了哪些更改。