在Python中拟合单调曲线

Question

我正在尝试将单调曲线拟合到一些近乎单调的数据上。 （X值是单调的，Y值应该是单调的，但噪声通常比点之间的基础值变化还要大。）这是我目前在做的事情的摘要：

def goodness_of_fit(Xfit):
    assert(is_sorted(Xfit))
    # ( Calculate the area between the fit line and the join-the-dots line from the data )

scipy.optimize.minimize(goodness_of_fit, x0=numpy.linspace(xmin, xmax))

我找不到一种方法来使优化算法保持Xfit数组排序 - 有人有什么建议吗？（由于数组的大小过大，创建N-1个单独的排序约束并使用受限制的优化函数变得不可行。）如果最佳方案只能在其他语言中实现，我愿意使用除Python之外的其他语言。

（注意：我确实拟合了X值而不是Y值 - 这是因为我最终想要绘制dX/dY曲线，并且使其不像从原始数据绘制时那样吹气成荒谬的价值。然而，如果固定X值上拟合Y值更容易，我可以选择这样做。）

- user1915639

---

我计划尝试的一个想法是使用LBFGSB（有界优化），并在每次迭代后更新边界。 - user1915639

Answer 1

如何从Xfit创建一个严格单调子集的新数组，然后将曲线拟合到该子集上呢？ 可以尝试以下方法：

Xfit      = np.hstack(((-1,-1,0,1,1),np.arange(1,10),(-7,9,10,10)))
len_mono  = 0
Xfit_mono = zeros(Xfit.size)
Xfit_mono_ind = zeros(Xfit.size)
Xfit_mono[len_mono] = Xfit[0]

for(iX, x) in enumerate(Xfit):
    if iX > 0:
        if x > Xfit_mono[len_mono]:
            len_mono = len_mono + 1
            Xfit_mono[len_mono] = x
            Xfit_mono_ind[len_mono] = iX 
Xfit_mono_ind = Xfit_mono_ind[:len_mono+1]
Xfit_mono = Xfit_mono[:len_mono+1]
print(Xfit_mono_ind)
print(Xfit_mono)

那么在进行曲线拟合时，您可以使用 Xfit_mono 替代原有的方法，并且当您需要选择相关联的 y 值时，您可以使用 Xfit_mono_ind 的值作为 y 的索引。

更新内容如下：

if x > Xfit[iX-1]:

使用：

if x > Xfit_mono[len_mono]:

对于递增的情况，旧解法是可行的。然而，如果出现递减的情况，旧结果将会给出错误的答案。这个更新后的解法在两种情况下都能得到期望的结果。

- BKay

---

你是指： 如果 x > Xfit_mono[len_mono] 而不是： 如果 x > Xfit[iX-1] ？ - user1915639

没错，我只测试了递增的集合，我的逻辑是正确的。然而，如果集合是递减的，你的建议仍然能提供正确的结果，而我的则不行。我已经进行了更新。 - BKay

这看起来对某些人可能有效；然而在我的情况下，X最有可能变得非单调的区域也是我最关心获得良好拟合的区域 - 我不想简单地跳过Y值。 - user1915639

Answer 2

最终，我通过适应连续X值之间的差异而不是值本身来完成它 - 然后我可以使用简单的下限为零。

def fg(X):
    # return the objective function and its Jacobian

def adjusted_fg(X_diff):
    X = cumsum(X_diff)
    score, jac = fg(X)
    jac[1:] = np.diff(jac)
    return score, jac

X0_diff[1:] = [X0[0]] + np.diff(X0)
bounds = [(None, None)] + [(0, None) for i in range(len(X0)-1)]

scipy.optimize.minimize(adjusted_fg, X0_diff, method='L-BFGS-B', bounds=bounds)

对于大的N，如果没有在接近解决方案的地方开始，很容易变得不稳定或陷入局部最小值，因此我将首先尝试适合较小的N（例如N/10），然后进行插值以获得更好的X0以适配大的N。

（注：对于我的问题，我实际上想要严格的排序，所以我使用了正下限而不是0）

- user1915639