曲线拟合 - 单调递增的导数

Question

我正在尝试为一些实验数据找到一个合适的曲线拟合。我的数据有两个特点：首先，y值会随着x值不断增加；其次，y对x的变化率（也就是dy/dx）也会不断增加。我试过很多种拟合方法，比如多项式拟合和单变量样条，但都没有得到我想要的效果。

所以，我在寻找一个可以帮助我定义最终曲线已知约束的曲线拟合函数（可能在scipy里有）。下面是我的一些数据示例，附带了一条拟合线，但这条线的变化率并没有表现出单调增加的特性。

import numpy as np
import matplotlib.pyplot as plt

data =  np.array([[  6.30991828, -10.22329935],
                  [  6.30991828, -10.2127338 ],
                  [  6.47697236, -10.01359361],
                  [  6.47697236,  -9.89353722],
                  [  6.47697236,  -9.81708052],
                  [  6.55108034,  -9.42113403],
                  [  6.55108034,  -9.21932801],
                  [  6.58617165,  -8.40428977],
                  [  6.62007321,  -7.6500927 ]])

interp = np.linspace(min(data[:,0]), max(data[:,0]), 20)
f = np.polyfit(data[:,0], data[:,-1], 3)
data_interp = np.polyval(f, interp)
plt.plot(data[:,0], data[:,1], 'x', interp, data_interp, '-')

补充：我相信在MATLAB中可以用slmengine来实现这个。

Answer 1

你的数据很有意思：你有三个不连续点，分别在6.30991828、6.47697236和6.55108034。这些点是真的吗？是你想要捕捉的内容吗？

没有任何连续的函数能够正确地表示这些不连续点。你唯一的办法是在不连续点的两侧进行分段拟合。这样你会得到三个不同的拟合：

1. x 小于 6.47697236
2. 6.47697236 小于 x 小于 6.55108034
3. x 大于 6.55108034

当然，在不连续点上，这个函数会有多个值。

如果这些不连续点对你来说没有意义，我建议你可以使用任何一个三次多项式，这样会得到一个连续的，并且一阶导数是递增的函数。

Answer 2

编辑：再试一次。我之前发的回答不太完整。而且我在阅读方面也没做好。我希望这次能更好。

from scipy.optimize import minimize
import numpy as np
import matplotlib.pyplot as plt

data = np.array([[  6.30991828, -10.22329935],
                  [  6.30991828, -10.2127338 ],
                  [  6.47697236, -10.01359361],
                  [  6.47697236,  -9.89353722],
                  [  6.47697236,  -9.81708052],
                  [  6.55108034,  -9.42113403],
                  [  6.55108034,  -9.21932801],
                  [  6.58617165,  -8.40428977],
                  [  6.62007321,  -7.6500927 ]])

x = data[:, 0]

def polynomial(p, x):
    return p[0]+p[1]*x+p[2]*x**2+p[3]*x**3

def constraint_2nd_der(p):
    return 2*p[2]+6*p[3]*x

def constraint_1st_der(p):
    return p[1]+2*p[2]*x+3*p[3]*x**2

def objective(p):
    return ((polynomial(p, x)-data[:, 1])**2).sum()


cons = (dict(type='ineq', fun=constraint_1st_der), dict(type='ineq', fun=constraint_2nd_der))
res = minimize(objective, x0=np.array([0., 0., 0., 0.]), method='SLSQP', constraints=cons)
if res.success:
    pars = res.x
    x = np.linspace(data[:, 0].min(), data[:, 0].max(), 100)
    pol = polynomial(pars, x)
    plt.plot(data[:, 0], data[:, 1], 'x', x, pol, '-')
    plt.show()
else:
    print 'Failed'