用条件方程拟合曲线
Curve fitting with conditional equation
问题
我创建了一个曲线拟合练习（见下面的功能代码），但我想添加功能。
我需要能够定义以下条件：slope at min(xdata) = 0.
（换句话说：我希望拟合曲线从水平梯度开始）
我试过的
我花了很多时间研究 scipy.optimize.curve_fit and evaluated other options (lmfit package 和 scipy 函数 scipy.optimize.fmin_slsqp, scipy.optimize.minimize，等等）。 lmfit 只允许我对参数设置一个静态条件，例如 p1 = 2 * p2 + 3。但它不允许我动态地处理 min(xdata)，我不能在约束中使用派生。
Scipy 只允许我最小化函数（找到一个最优的 x，但参数 p 是已知的）。或者它可以用于定义参数的特定范围。我无法定义可用于在曲线拟合期间约束参数的第二个函数。
我需要能够将条件直接传递给曲线拟合算法（而不是通过将条件带入 cubic_fit() 方程来解决问题 - 似乎可以消除例如 p3 并将其定义为其他参数和 min(xdata)) 的组合。我的实际拟合函数要复杂得多，我需要 运行 这个脚本在一批数据上迭代（变化 min(xdata)）。我不能每次都手动更改拟合函数...
我很感激任何建议，也许还有其他软件包可以对曲线拟合问题进行更复杂的定义？
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
import scipy.optimize
# generate dummy data - on which I will run a curve fit below
def cubic_fit_with_noise(x, p1, p2, p3, p4):
return p1 + p2*x + p3*x**2 + p4*x**3 + np.random.rand()
xdata = [x * 0.1 for x in range(0, 100)]
ydata = np.array( [cubic_fit_with_noise (x, 2, 0.4, -.2,0.02) for x in xdata] )
# now, run the curve-fit
# set up the fitting function:
def cubic_fit(x, p1, p2, p3, p4):
return p1 + p2*x + p3*x**2 + p4*x**3
# define starting point:
s1 = 2.5
s2 = 0.2
s3 = -.2
s4 = 0.02
# scipy curve fitting:
popt, pcov = scipy.optimize.curve_fit(cubic_fit, xdata, ydata, p0=(s1,s2,s3,s4))
y_modelled = np.array([cubic_fit(x, popt[0], popt[1], popt[2], popt[3]) for x in xdata])
print(popt) # prints out the 4 parameters p1,p2,p3,p4 defined in curve-fitting
plt.plot(xdata, ydata, 'bo')
plt.plot(xdata, y_modelled, 'r-')
plt.show()
上面的代码 运行s 和 Python3（如果你有 Python2，修正打印语句）。
作为补充，我想引入导数：
def cubic_fit_derivative(x, p1, p2, p3, p4):
return p2 + 2.0 * p3 * x + 3 * p4 * x**2
和 cubic_fit_derivative(min(xdata), p1,p2,p3,p4) = 0.
的约束
你的多项式导数在 xmin 处的条件 = 0 可以表示为一个简单的约束，这意味着变量 p2、p3 和 p4 实际上并不独立。求导条件为
p2 + 2*p3*xmin + 3*p4*xmin**2 = 0
其中 xmin 是 xdata 的最小值。此外，xmin 将在适合之前已知（如果在编写脚本时不一定），您可以使用它来约束三个参数之一。由于 xmin 可能为零（实际上，这是针对您的情况），因此约束应该是
p2 = - 2*p3*xmin - 3*p4*xmin**2
使用 lmfit，原始的、无约束的拟合看起来像这样（我稍微清理了一下）：
import numpy as np
from lmfit import Model
import matplotlib.pylab as plt
# the model function:
def cubic_poly(x, p1, p2, p3, p4):
return p1 + p2*x + p3*x**2 + p4*x**3
xdata = np.arange(100) * 0.1
ydata = cubic_poly(xdata, 2, 0.4, -.2, 0.02)
ydata = ydata + np.random.normal(size=len(xdata), scale=0.05)
# make Model, create parameters, run fit, print results
model = Model(cubic_poly)
params = model.make_params(p1=2.5, p2=0.2, p3=-0.0, p4=0.0)
result = model.fit(ydata, params, x=xdata)
print(result.fit_report())
plt.plot(xdata, ydata, 'bo')
plt.plot(xdata, result.best_fit, 'r-')
plt.show()
打印：
[[Model]]
Model(cubic_poly)
[[Fit Statistics]]
# function evals = 13
# data points = 100
# variables = 4
chi-square = 0.218
reduced chi-square = 0.002
Akaike info crit = -604.767
Bayesian info crit = -594.347
[[Variables]]
p1: 2.00924432 +/- 0.018375 (0.91%) (init= 2.5)
p2: 0.39427207 +/- 0.016155 (4.10%) (init= 0.2)
p3: -0.19902928 +/- 0.003802 (1.91%) (init=-0)
p4: 0.01993319 +/- 0.000252 (1.27%) (init= 0)
[[Correlations]] (unreported correlations are < 0.100)
C(p3, p4) = -0.986
C(p2, p3) = -0.967
C(p2, p4) = 0.914
C(p1, p2) = -0.857
C(p1, p3) = 0.732
C(p1, p4) = -0.646
并生成一个图
现在，要添加你的约束条件，我们将添加xmin作为固定参数，并像上面那样约束p2，将上面的替换为：
params = model.make_params(p1=2.5, p2=0.2, p3=-0.0, p4=0.0)
# add an extra parameter for `xmin`
params.add('xmin', min(xdata), vary=False)
# constrain p2 so that the derivative is 0 at xmin
params['p2'].expr = '-2*p3*xmin - 3*p4*xmin**2'
result = model.fit(ydata, params, x=xdata)
print(result.fit_report())
plt.plot(xdata, ydata, 'bo')
plt.plot(xdata, result.best_fit, 'r-')
plt.show()
现在打印
[[Model]]
Model(cubic_poly)
[[Fit Statistics]]
# function evals = 10
# data points = 100
# variables = 3
chi-square = 1.329
reduced chi-square = 0.014
Akaike info crit = -426.056
Bayesian info crit = -418.241
[[Variables]]
p1: 2.39001759 +/- 0.023239 (0.97%) (init= 2.5)
p2: 0 +/- 0 (nan%) == '-2*p3*xmin - 3*p4*xmin**2'
p3: -0.10858258 +/- 0.002372 (2.19%) (init=-0)
p4: 0.01424411 +/- 0.000251 (1.76%) (init= 0)
xmin: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
C(p3, p4) = -0.986
C(p1, p3) = -0.742
C(p1, p4) = 0.658
和
这样的情节
如果xmin不为零（比如xdata = np.linspace(-10, 10, 101)，p2的值和不确定性就不会为零。
正如我在评论中提到的，您只需要安装正确的功能即可。不过我忘记了常量。所以函数将是 a*(x-xmin)**2*(x-xn)+c
As curvefit 不采用其他参数，例如leatssq，唯一的技巧就是通过xmin。我通过一个全局变量来做到这一点（也许不是最好的方法，但它有效。欢迎就如何做得更好发表评论）。
最后，您只需在代码中添加以下几行：
def cubic_zero(x,a,xn,const):
global xmin
return (a*(x-xmin)**2*(x-xn)+const)
和
xmin=xdata[0]
popt2, pcov2 = scipy.optimize.curve_fit(cubic_zero, xdata, ydata)
y_modelled2 = np.array([cubic_zero(x, *popt2) for x in xdata])
print(popt2)
plt.plot(xdata, y_modelled2, color='#ee9900',linestyle="--")
提供
>>>[ 0.01429367 7.63190327 2.92604132]
和
此解决方案使用 scipy.optimize.leastsq。使用自制的 residuals 函数，实际上不需要将 xmin 作为附加参数传递给拟合。 fit 函数与另一个 post 中一样，因此不需要约束。这看起来像：
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
def cubic_fit_with_noise(x, p1, p2, p3, p4):
return p1 + p2*x + p3*x**2 + p4*x**3 + .2*(1-2*np.random.rand())
def cubic_zero(x,a,xn,const, xmin):
return (a*(x-xmin)**2*(x-xn)+const)
def residuals(params, dataX,dataY):
a,xn,const=params
xmin=dataX[0]
dist=np.fromiter( (y-cubic_zero(x,a,xn,const, xmin) for x,y in zip(dataX,dataY)), np.float)
return dist
xdata = np.linspace(.5,10.5,100)
ydata = np.fromiter( (cubic_fit_with_noise (x, 2, 0.4, -.2,0.02) for x in xdata), np.float )
# scipy curve fitting with leastsq:
initialGuess=[.3,.3,.3]
popt2, pcov2, info2, msg2, ier2 = leastsq(residuals,initialGuess, args=(xdata, ydata), full_output=True)
fullparams=np.append(popt2,xdata[0])
y_modelled2 = np.array([cubic_zero(x, *fullparams) for x in xdata])
print(popt2)
print(pcov2)
print np.array([ -popt2[0]*xdata[0]**2*popt2[1]+popt2[2],popt2[0]*(xdata[0]**2+2*xdata[0]*popt2[1]),-popt2[0]*(2*xdata[0]+popt2[1]),popt2[0] ])
plt.plot(xdata, ydata, 'bo')
plt.plot(xdata, y_modelled2, 'r-')
plt.show()
并提供：
>>>[ 0.01710749 7.69369653 2.38986378]
>>>[[ 4.33308441e-06 5.61402017e-04 2.71819763e-04]
[ 5.61402017e-04 1.10367937e-01 5.67852980e-02]
[ 2.71819763e-04 5.67852980e-02 3.94127702e-02]]
>>>[ 2.35695882 0.13589672 -0.14872733 0.01710749]
目前无法上传图片...无论出于何种原因，但结果与另一个相同post
问题
我创建了一个曲线拟合练习（见下面的功能代码），但我想添加功能。
我需要能够定义以下条件：slope at min(xdata) = 0.
（换句话说：我希望拟合曲线从水平梯度开始）
我试过的
我花了很多时间研究 scipy.optimize.curve_fit and evaluated other options (lmfit package 和 scipy 函数 scipy.optimize.fmin_slsqp, scipy.optimize.minimize，等等）。 lmfit 只允许我对参数设置一个静态条件，例如 p1 = 2 * p2 + 3。但它不允许我动态地处理 min(xdata)，我不能在约束中使用派生。
Scipy 只允许我最小化函数（找到一个最优的 x，但参数 p 是已知的）。或者它可以用于定义参数的特定范围。我无法定义可用于在曲线拟合期间约束参数的第二个函数。
我需要能够将条件直接传递给曲线拟合算法（而不是通过将条件带入 cubic_fit() 方程来解决问题 - 似乎可以消除例如 p3 并将其定义为其他参数和 min(xdata)) 的组合。我的实际拟合函数要复杂得多，我需要 运行 这个脚本在一批数据上迭代（变化 min(xdata)）。我不能每次都手动更改拟合函数...
我很感激任何建议，也许还有其他软件包可以对曲线拟合问题进行更复杂的定义？
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
import scipy.optimize
# generate dummy data - on which I will run a curve fit below
def cubic_fit_with_noise(x, p1, p2, p3, p4):
return p1 + p2*x + p3*x**2 + p4*x**3 + np.random.rand()
xdata = [x * 0.1 for x in range(0, 100)]
ydata = np.array( [cubic_fit_with_noise (x, 2, 0.4, -.2,0.02) for x in xdata] )
# now, run the curve-fit
# set up the fitting function:
def cubic_fit(x, p1, p2, p3, p4):
return p1 + p2*x + p3*x**2 + p4*x**3
# define starting point:
s1 = 2.5
s2 = 0.2
s3 = -.2
s4 = 0.02
# scipy curve fitting:
popt, pcov = scipy.optimize.curve_fit(cubic_fit, xdata, ydata, p0=(s1,s2,s3,s4))
y_modelled = np.array([cubic_fit(x, popt[0], popt[1], popt[2], popt[3]) for x in xdata])
print(popt) # prints out the 4 parameters p1,p2,p3,p4 defined in curve-fitting
plt.plot(xdata, ydata, 'bo')
plt.plot(xdata, y_modelled, 'r-')
plt.show()
上面的代码 运行s 和 Python3（如果你有 Python2，修正打印语句）。 作为补充，我想引入导数：
def cubic_fit_derivative(x, p1, p2, p3, p4):
return p2 + 2.0 * p3 * x + 3 * p4 * x**2
和 cubic_fit_derivative(min(xdata), p1,p2,p3,p4) = 0.
你的多项式导数在 xmin 处的条件 = 0 可以表示为一个简单的约束，这意味着变量 p2、p3 和 p4 实际上并不独立。求导条件为
p2 + 2*p3*xmin + 3*p4*xmin**2 = 0
其中 xmin 是 xdata 的最小值。此外，xmin 将在适合之前已知（如果在编写脚本时不一定），您可以使用它来约束三个参数之一。由于 xmin 可能为零（实际上，这是针对您的情况），因此约束应该是
p2 = - 2*p3*xmin - 3*p4*xmin**2
使用 lmfit，原始的、无约束的拟合看起来像这样（我稍微清理了一下）：
import numpy as np
from lmfit import Model
import matplotlib.pylab as plt
# the model function:
def cubic_poly(x, p1, p2, p3, p4):
return p1 + p2*x + p3*x**2 + p4*x**3
xdata = np.arange(100) * 0.1
ydata = cubic_poly(xdata, 2, 0.4, -.2, 0.02)
ydata = ydata + np.random.normal(size=len(xdata), scale=0.05)
# make Model, create parameters, run fit, print results
model = Model(cubic_poly)
params = model.make_params(p1=2.5, p2=0.2, p3=-0.0, p4=0.0)
result = model.fit(ydata, params, x=xdata)
print(result.fit_report())
plt.plot(xdata, ydata, 'bo')
plt.plot(xdata, result.best_fit, 'r-')
plt.show()
打印：
[[Model]]
Model(cubic_poly)
[[Fit Statistics]]
# function evals = 13
# data points = 100
# variables = 4
chi-square = 0.218
reduced chi-square = 0.002
Akaike info crit = -604.767
Bayesian info crit = -594.347
[[Variables]]
p1: 2.00924432 +/- 0.018375 (0.91%) (init= 2.5)
p2: 0.39427207 +/- 0.016155 (4.10%) (init= 0.2)
p3: -0.19902928 +/- 0.003802 (1.91%) (init=-0)
p4: 0.01993319 +/- 0.000252 (1.27%) (init= 0)
[[Correlations]] (unreported correlations are < 0.100)
C(p3, p4) = -0.986
C(p2, p3) = -0.967
C(p2, p4) = 0.914
C(p1, p2) = -0.857
C(p1, p3) = 0.732
C(p1, p4) = -0.646
并生成一个图
现在，要添加你的约束条件，我们将添加xmin作为固定参数，并像上面那样约束p2，将上面的替换为：
params = model.make_params(p1=2.5, p2=0.2, p3=-0.0, p4=0.0)
# add an extra parameter for `xmin`
params.add('xmin', min(xdata), vary=False)
# constrain p2 so that the derivative is 0 at xmin
params['p2'].expr = '-2*p3*xmin - 3*p4*xmin**2'
result = model.fit(ydata, params, x=xdata)
print(result.fit_report())
plt.plot(xdata, ydata, 'bo')
plt.plot(xdata, result.best_fit, 'r-')
plt.show()
现在打印
[[Model]]
Model(cubic_poly)
[[Fit Statistics]]
# function evals = 10
# data points = 100
# variables = 3
chi-square = 1.329
reduced chi-square = 0.014
Akaike info crit = -426.056
Bayesian info crit = -418.241
[[Variables]]
p1: 2.39001759 +/- 0.023239 (0.97%) (init= 2.5)
p2: 0 +/- 0 (nan%) == '-2*p3*xmin - 3*p4*xmin**2'
p3: -0.10858258 +/- 0.002372 (2.19%) (init=-0)
p4: 0.01424411 +/- 0.000251 (1.76%) (init= 0)
xmin: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
C(p3, p4) = -0.986
C(p1, p3) = -0.742
C(p1, p4) = 0.658
和
这样的情节如果xmin不为零（比如xdata = np.linspace(-10, 10, 101)，p2的值和不确定性就不会为零。
正如我在评论中提到的，您只需要安装正确的功能即可。不过我忘记了常量。所以函数将是 a*(x-xmin)**2*(x-xn)+c
As curvefit 不采用其他参数，例如leatssq，唯一的技巧就是通过xmin。我通过一个全局变量来做到这一点（也许不是最好的方法，但它有效。欢迎就如何做得更好发表评论）。
最后，您只需在代码中添加以下几行：
def cubic_zero(x,a,xn,const):
global xmin
return (a*(x-xmin)**2*(x-xn)+const)
和
xmin=xdata[0]
popt2, pcov2 = scipy.optimize.curve_fit(cubic_zero, xdata, ydata)
y_modelled2 = np.array([cubic_zero(x, *popt2) for x in xdata])
print(popt2)
plt.plot(xdata, y_modelled2, color='#ee9900',linestyle="--")
提供
>>>[ 0.01429367 7.63190327 2.92604132]
和
此解决方案使用 scipy.optimize.leastsq。使用自制的 residuals 函数，实际上不需要将 xmin 作为附加参数传递给拟合。 fit 函数与另一个 post 中一样，因此不需要约束。这看起来像：
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
def cubic_fit_with_noise(x, p1, p2, p3, p4):
return p1 + p2*x + p3*x**2 + p4*x**3 + .2*(1-2*np.random.rand())
def cubic_zero(x,a,xn,const, xmin):
return (a*(x-xmin)**2*(x-xn)+const)
def residuals(params, dataX,dataY):
a,xn,const=params
xmin=dataX[0]
dist=np.fromiter( (y-cubic_zero(x,a,xn,const, xmin) for x,y in zip(dataX,dataY)), np.float)
return dist
xdata = np.linspace(.5,10.5,100)
ydata = np.fromiter( (cubic_fit_with_noise (x, 2, 0.4, -.2,0.02) for x in xdata), np.float )
# scipy curve fitting with leastsq:
initialGuess=[.3,.3,.3]
popt2, pcov2, info2, msg2, ier2 = leastsq(residuals,initialGuess, args=(xdata, ydata), full_output=True)
fullparams=np.append(popt2,xdata[0])
y_modelled2 = np.array([cubic_zero(x, *fullparams) for x in xdata])
print(popt2)
print(pcov2)
print np.array([ -popt2[0]*xdata[0]**2*popt2[1]+popt2[2],popt2[0]*(xdata[0]**2+2*xdata[0]*popt2[1]),-popt2[0]*(2*xdata[0]+popt2[1]),popt2[0] ])
plt.plot(xdata, ydata, 'bo')
plt.plot(xdata, y_modelled2, 'r-')
plt.show()
并提供：
>>>[ 0.01710749 7.69369653 2.38986378]
>>>[[ 4.33308441e-06 5.61402017e-04 2.71819763e-04]
[ 5.61402017e-04 1.10367937e-01 5.67852980e-02]
[ 2.71819763e-04 5.67852980e-02 3.94127702e-02]]
>>>[ 2.35695882 0.13589672 -0.14872733 0.01710749]
目前无法上传图片...无论出于何种原因，但结果与另一个相同post