曲线拟合参数范围
Curve fit with parameter bounds
我有实验数据：
xdata = [85,86,87,88,89,90,91,91.75,93,96,100,101,102,103,104,105,106,107.25,108.25,109,109.75,111,112,112.75,114,115.25,116,116.75,118,119.25,120,121,122,122.5,123.5,125.25,126,126.75,127.75,129.25,130.25,131,131.75,133,134.25,135,136,137,138,139,140,141,142,143,144,144.75,146,146.75,148,149.25,150,150.5,152,153.25,154,155,156.75,158,159,159.75,161,162,162.5,164,165,166]
ydata = [0.2,0.21,0.18,0.21,0.19,0.2,0.21,0.2,0.18,0.204,0.208,0.2,0.21,0.25,0.2,0.19,0.216,0.22,0.224,0.26,0.229,0.237,0.22,0.246,0.25,0.264,0.29,0.274,0.29,0.3,0.27,0.32,0.38,0.348,0.372,0.398,0.35,0.42,0.444,0.48,0.496,0.55,0.51,0.54,0.57,0.51,0.605,0.57,0.65,0.642,0.6,0.66,0.7,0.688,0.69,0.705,0.67,0.717,0.69,0.728,0.75,0.736,0.73,0.744,0.72,0.76,0.752,0.74,0.76,0.7546,0.77,0.74,0.758,0.74,0.78,0.76]
和公式f(x) = m1 + m2 / (1 + e ^ (-m3*(x - m4)))。我需要用最小二乘法找到 m1,
m2, m3, m4，其中
0.05 < m1 < 0.3
0.3 < 平方米 < 0.8
0.05 < 立方米 < 0.5
100 < 立方米 < 200.
我使用 curve_fit，我的函数是：
def f(xdata, m1, m2, m3, m4):
if m1 > 0.05 and m1 < 0.3 and \
m2 > 0.3 and m2 < 0.8 and \
m3 > 0.05 and m3 < 0.5 and \
m4 > 100 and m4 < 200:
return m1 + (m2 * 1. / (1 + e ** (-m3 * (x - m4))))
return (abs(m1) + abs(m2) + abs(m3) + abs(m4)) * 1e14 # some large number
但是程序return报错：RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 1000.
怎么办？
import numpy as np
from scipy.optimize import curve_fit
from math import e
xdata = np.array([85,86,87,88,89,90,91,91.75,93,96,100,101,102,103,104,105,106,107.25,108.25,109,109.75,111,112,112.75,114,115.25,116,116.75,118,119.25,120,121,122,122.5,123.5,125.25,126,126.75,127.75,129.25,130.25,131,131.75,133,134.25,135,136,137,138,139,140,141,142,143,144,144.75,146,146.75,148,149.25,150,150.5,152,153.25,154,155,156.75,158,159,159.75,161,162,162.5,164,165,166])`
ydata = np.array([0.2,0.21,0.18,0.21,0.19,0.2,0.21,0.2,0.18,0.204,0.208,0.2,0.21,0.25,0.2,0.19,0.216,0.22,0.224,0.26,0.229,0.237,0.22,0.246,0.25,0.264,0.29,0.274,0.29,0.3,0.27,0.32,0.38,0.348,0.372,0.398,0.35,0.42,0.444,0.48,0.496,0.55,0.51,0.54,0.57,0.51,0.605,0.57,0.65,0.642,0.6,0.66,0.7,0.688,0.69,0.705,0.67,0.717,0.69,0.728,0.75,0.736,0.73,0.744,0.72,0.76,0.752,0.74,0.76,0.7546,0.77,0.74,0.758,0.74,0.78,0.76])
def f(xdata, m1, m2, m3, m4):
if m1 > 0.05 and m1 < 0.3 and \
m2 > 0.3 and m2 < 0.8 and \
m3 > 0.05 and m3 < 0.5 and \
m4 > 100 and m4 < 200:
return m1 + (m2 * 1. / (1 + e ** (-m3 * (x - m4))))
return (abs(m1) + abs(m2) + abs(m3) + abs(m4)) * 1e14
print curve_fit(f, xdata, ydata)
将初始参数设置为有用的值：
curve_fit(f, xdata, ydata, p0=(0.1, 0.5, 0.1, 150)))
此外，在您的函数中使用 xdata 而不是 x f:
return m1 + (m2 * 1. / (1 + e ** (-m3 * (xdata - m4))))
这是我修改的程序：
def f(xdata, m1, m2, m3, m4):
if (0.05 < m1 < 0.3 and
0.3 < m2 < 0.8 and
0.05 < m3 < 0.5 and
100 < m4 < 200):
return m1 + (m2 * 1. / (1 + e ** (-m3 * (xdata - m4))))
return 1e38
print(curve_fit(f, xdata, ydata, p0=(0.1, 0.5, 0.1, 150)))
结果：
(array([ 0.19567035, 0.56792559, 0.13434829, 129.98915877]),
array([[ 2.94622909e-05, -3.96126279e-05, 1.99236054e-05,
7.48438125e-04],
[ -3.96126279e-05, 9.24145662e-05, -4.62302643e-05,
5.04671621e-04],
[ 1.99236054e-05, -4.62302643e-05, 3.77364832e-05,
-2.43866126e-04],
[ 7.48438125e-04, 5.04671621e-04, -2.43866126e-04,
1.34700612e-01]]))
或者，您也可以使用 lmfit，它允许您轻松设置边界并避免在函数中使用 "ugly" if 语句。您获得的参数如下：
m1: 0.19567033 +/- 0.005427 (2.77%) (init= 0.1)
m2: 0.56792558 +/- 0.009613 (1.69%) (init= 0.5)
m3: 0.13434829 +/- 0.006143 (4.57%) (init= 0.2)
m4: 129.989156 +/- 0.367009 (0.28%) (init= 150)
您获得的输出如下所示：
这是完整的代码，有几条注释；如果您还有其他问题，请告诉我：
from lmfit import minimize, Parameters, Parameter, report_fit
import numpy as np
xdata = np.array([85,86,87,88,89,90,91,91.75,93,96,100,101,102,103,104,105,106,107.25,108.25,109,109.75,111,112,112.75,114,115.25,116,116.75,118,119.25,120,121,122,122.5,123.5,125.25,126,126.75,127.75,129.25,130.25,131,131.75,133,134.25,135,136,137,138,139,140,141,142,143,144,144.75,146,146.75,148,149.25,150,150.5,152,153.25,154,155,156.75,158,159,159.75,161,162,162.5,164,165,166])
ydata = np.array([0.2,0.21,0.18,0.21,0.19,0.2,0.21,0.2,0.18,0.204,0.208,0.2,0.21,0.25,0.2,0.19,0.216,0.22,0.224,0.26,0.229,0.237,0.22,0.246,0.25,0.264,0.29,0.274,0.29,0.3,0.27,0.32,0.38,0.348,0.372,0.398,0.35,0.42,0.444,0.48,0.496,0.55,0.51,0.54,0.57,0.51,0.605,0.57,0.65,0.642,0.6,0.66,0.7,0.688,0.69,0.705,0.67,0.717,0.69,0.728,0.75,0.736,0.73,0.744,0.72,0.76,0.752,0.74,0.76,0.7546,0.77,0.74,0.758,0.74,0.78,0.76])
def fit_fc(params, x, data):
m1 = params['m1'].value
m2 = params['m2'].value
m3 = params['m3'].value
m4 = params['m4'].value
model = m1 + (m2 * 1. / (1 + np.exp(-m3 * (x - m4))))
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('m1', value= 0.1, min=0.05, max=0.3)
params.add('m2', value= 0.5, min=0.3, max=0.8)
params.add('m3', value= 0.2, min=0.05, max=0.5)
params.add('m4', value= 150.0, min=100, max=200)
# do fit, here with leastsq model
result = minimize(fit_fc, params, args=(xdata, ydata))
# calculate final result
final = ydata + result.residual
# write error report
report_fit(params)
#plot results
try:
import pylab
pylab.plot(xdata, ydata, 'k+')
pylab.plot(xdata, final, 'r')
pylab.show()
except:
pass
我有实验数据：
xdata = [85,86,87,88,89,90,91,91.75,93,96,100,101,102,103,104,105,106,107.25,108.25,109,109.75,111,112,112.75,114,115.25,116,116.75,118,119.25,120,121,122,122.5,123.5,125.25,126,126.75,127.75,129.25,130.25,131,131.75,133,134.25,135,136,137,138,139,140,141,142,143,144,144.75,146,146.75,148,149.25,150,150.5,152,153.25,154,155,156.75,158,159,159.75,161,162,162.5,164,165,166]
ydata = [0.2,0.21,0.18,0.21,0.19,0.2,0.21,0.2,0.18,0.204,0.208,0.2,0.21,0.25,0.2,0.19,0.216,0.22,0.224,0.26,0.229,0.237,0.22,0.246,0.25,0.264,0.29,0.274,0.29,0.3,0.27,0.32,0.38,0.348,0.372,0.398,0.35,0.42,0.444,0.48,0.496,0.55,0.51,0.54,0.57,0.51,0.605,0.57,0.65,0.642,0.6,0.66,0.7,0.688,0.69,0.705,0.67,0.717,0.69,0.728,0.75,0.736,0.73,0.744,0.72,0.76,0.752,0.74,0.76,0.7546,0.77,0.74,0.758,0.74,0.78,0.76]
和公式f(x) = m1 + m2 / (1 + e ^ (-m3*(x - m4)))。我需要用最小二乘法找到 m1,
m2, m3, m4，其中
0.05 < m1 < 0.3
0.3 < 平方米 < 0.8
0.05 < 立方米 < 0.5
100 < 立方米 < 200.
我使用 curve_fit，我的函数是：
def f(xdata, m1, m2, m3, m4):
if m1 > 0.05 and m1 < 0.3 and \
m2 > 0.3 and m2 < 0.8 and \
m3 > 0.05 and m3 < 0.5 and \
m4 > 100 and m4 < 200:
return m1 + (m2 * 1. / (1 + e ** (-m3 * (x - m4))))
return (abs(m1) + abs(m2) + abs(m3) + abs(m4)) * 1e14 # some large number
但是程序return报错：RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 1000.
怎么办？
import numpy as np
from scipy.optimize import curve_fit
from math import e
xdata = np.array([85,86,87,88,89,90,91,91.75,93,96,100,101,102,103,104,105,106,107.25,108.25,109,109.75,111,112,112.75,114,115.25,116,116.75,118,119.25,120,121,122,122.5,123.5,125.25,126,126.75,127.75,129.25,130.25,131,131.75,133,134.25,135,136,137,138,139,140,141,142,143,144,144.75,146,146.75,148,149.25,150,150.5,152,153.25,154,155,156.75,158,159,159.75,161,162,162.5,164,165,166])`
ydata = np.array([0.2,0.21,0.18,0.21,0.19,0.2,0.21,0.2,0.18,0.204,0.208,0.2,0.21,0.25,0.2,0.19,0.216,0.22,0.224,0.26,0.229,0.237,0.22,0.246,0.25,0.264,0.29,0.274,0.29,0.3,0.27,0.32,0.38,0.348,0.372,0.398,0.35,0.42,0.444,0.48,0.496,0.55,0.51,0.54,0.57,0.51,0.605,0.57,0.65,0.642,0.6,0.66,0.7,0.688,0.69,0.705,0.67,0.717,0.69,0.728,0.75,0.736,0.73,0.744,0.72,0.76,0.752,0.74,0.76,0.7546,0.77,0.74,0.758,0.74,0.78,0.76])
def f(xdata, m1, m2, m3, m4):
if m1 > 0.05 and m1 < 0.3 and \
m2 > 0.3 and m2 < 0.8 and \
m3 > 0.05 and m3 < 0.5 and \
m4 > 100 and m4 < 200:
return m1 + (m2 * 1. / (1 + e ** (-m3 * (x - m4))))
return (abs(m1) + abs(m2) + abs(m3) + abs(m4)) * 1e14
print curve_fit(f, xdata, ydata)
将初始参数设置为有用的值：
curve_fit(f, xdata, ydata, p0=(0.1, 0.5, 0.1, 150)))
此外，在您的函数中使用 xdata 而不是 x f:
return m1 + (m2 * 1. / (1 + e ** (-m3 * (xdata - m4))))
这是我修改的程序：
def f(xdata, m1, m2, m3, m4):
if (0.05 < m1 < 0.3 and
0.3 < m2 < 0.8 and
0.05 < m3 < 0.5 and
100 < m4 < 200):
return m1 + (m2 * 1. / (1 + e ** (-m3 * (xdata - m4))))
return 1e38
print(curve_fit(f, xdata, ydata, p0=(0.1, 0.5, 0.1, 150)))
结果：
(array([ 0.19567035, 0.56792559, 0.13434829, 129.98915877]),
array([[ 2.94622909e-05, -3.96126279e-05, 1.99236054e-05,
7.48438125e-04],
[ -3.96126279e-05, 9.24145662e-05, -4.62302643e-05,
5.04671621e-04],
[ 1.99236054e-05, -4.62302643e-05, 3.77364832e-05,
-2.43866126e-04],
[ 7.48438125e-04, 5.04671621e-04, -2.43866126e-04,
1.34700612e-01]]))
或者，您也可以使用 lmfit，它允许您轻松设置边界并避免在函数中使用 "ugly" if 语句。您获得的参数如下：
m1: 0.19567033 +/- 0.005427 (2.77%) (init= 0.1)
m2: 0.56792558 +/- 0.009613 (1.69%) (init= 0.5)
m3: 0.13434829 +/- 0.006143 (4.57%) (init= 0.2)
m4: 129.989156 +/- 0.367009 (0.28%) (init= 150)
您获得的输出如下所示：
这是完整的代码，有几条注释；如果您还有其他问题，请告诉我：
from lmfit import minimize, Parameters, Parameter, report_fit
import numpy as np
xdata = np.array([85,86,87,88,89,90,91,91.75,93,96,100,101,102,103,104,105,106,107.25,108.25,109,109.75,111,112,112.75,114,115.25,116,116.75,118,119.25,120,121,122,122.5,123.5,125.25,126,126.75,127.75,129.25,130.25,131,131.75,133,134.25,135,136,137,138,139,140,141,142,143,144,144.75,146,146.75,148,149.25,150,150.5,152,153.25,154,155,156.75,158,159,159.75,161,162,162.5,164,165,166])
ydata = np.array([0.2,0.21,0.18,0.21,0.19,0.2,0.21,0.2,0.18,0.204,0.208,0.2,0.21,0.25,0.2,0.19,0.216,0.22,0.224,0.26,0.229,0.237,0.22,0.246,0.25,0.264,0.29,0.274,0.29,0.3,0.27,0.32,0.38,0.348,0.372,0.398,0.35,0.42,0.444,0.48,0.496,0.55,0.51,0.54,0.57,0.51,0.605,0.57,0.65,0.642,0.6,0.66,0.7,0.688,0.69,0.705,0.67,0.717,0.69,0.728,0.75,0.736,0.73,0.744,0.72,0.76,0.752,0.74,0.76,0.7546,0.77,0.74,0.758,0.74,0.78,0.76])
def fit_fc(params, x, data):
m1 = params['m1'].value
m2 = params['m2'].value
m3 = params['m3'].value
m4 = params['m4'].value
model = m1 + (m2 * 1. / (1 + np.exp(-m3 * (x - m4))))
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('m1', value= 0.1, min=0.05, max=0.3)
params.add('m2', value= 0.5, min=0.3, max=0.8)
params.add('m3', value= 0.2, min=0.05, max=0.5)
params.add('m4', value= 150.0, min=100, max=200)
# do fit, here with leastsq model
result = minimize(fit_fc, params, args=(xdata, ydata))
# calculate final result
final = ydata + result.residual
# write error report
report_fit(params)
#plot results
try:
import pylab
pylab.plot(xdata, ydata, 'k+')
pylab.plot(xdata, final, 'r')
pylab.show()
except:
pass