用条件方程拟合曲线

Question

问题

我已经创建了一个曲线拟合练习(见下面的函数代码)，但是我想添加到这个功能中。

我需要能够定义以下条件：slope at min(xdata) = 0。

(换句话说:我希望拟合的曲线从水平梯度开始)

我已经尝试过的

我花了相当多的时间研究适合，并评估了其他选项(lmfit包装和bit函数scipy.optimize.fmin_slsqp、scipy.optimize.minimize等)。lmfit只允许我对参数设置一个静态条件，例如p1 =2* p2 + 3，但是它不允许我动态地寻址min(xdata)，并且不能在约束中使用派生。

Scipy只允许我最小化函数(找到一个最优的x，但参数p已经知道)。也可以用来定义参数的特定范围。我无法定义第二个函数，可以用来约束曲线拟合过程中的参数。

我需要能够将条件直接传递给曲线拟合算法(而不是通过将条件引入cubic_fit()方程来解决这个问题-似乎可以消除p3并将其定义为其他参数和min的组合(Xdata))。我的实际拟合函数要复杂得多，我需要在一批数据上迭代地运行这个脚本(varying (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()

以上代码使用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的约束。

Answer 1

您的条件是多项式=0在xmin处的导数可以表示为一个简单的约束，这意味着变量p2、p3和p4实际上不是独立的。衍生条件是

p2 + 2*p3*xmin + 3*p4*xmin**2 = 0

其中xmin是xdata的最小值。此外，xmin将在fit之前被知道(如果在编写脚本时不一定)，您可以使用它来约束这三个参数之一。由于xmin可能为零(实际上，它适用于您的情况)，因此约束应该是

p2 = - 2*p3*xmin - 3*p4*xmin**2

使用lmfit时，原始的、无约束的fit看起来如下(我稍微清理了一下)：

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的值和不确定性就不会为零。

Answer 2

正如我在评论中提到的那样，您只需符合正确的功能即可。不过，我忘了常量。所以函数是a*(x-xmin)**2*(x-xn)+c

由于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]

和

​

​

Answer 3

此解决方案使用scipy.optimize.leastsq。使用自制的residuals函数，实际上不需要将xmin作为附加参数传递给fit。fit功能和另一个职位一样，因此没有约束的必要。这看起来像是：

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]

图像上传目前不起作用..。无论出于什么原因，结果都与另一职位相同。