多元问题的分段拟合

Question

我正在尝试使用 scipy curve_fit 中的多元拟合来拟合 3 维数据（即 2 个自变量和 1 个因变量）。我希望对同样的问题进行分段拟合。我试图在of this 的基础上继续进行，但没有成功。问题定义如下：

import numpy as np
from scipy.optimize import curve_fit
#..........................................................................................................
def F0(X, a, b, c, c0, y0):
    x, y = X
    value = []
    for i in range(0, len(x)):
        if y[i] < y0:
            lnZ = x[i] + c0*y[i]
        else:
            lnZ = x[i] + c*y[i]
        val = a + (b*lnZ)
        value.append(val)
    return value
#..........................................................................................................
def F1(X, a, b, c):
    x, y = X
    lnZ = x + c*y
    value = a + (b*lnZ)
    return value
#..........................................................................................................
x = [-2.302585093,
-2.302585093, 
-2.302585093, 
-2.302585093, 
-2.302585093,
-2.302585093,
-2.302585093,
0,
0,
0,
0,
0,
0,
0,
2.302585093,
2.302585093,
2.302585093,
2.302585093,
2.302585093,
2.302585093,
2.302585093
]

y = [7.55E-04,
7.85E-04,
8.17E-04,
8.52E-04,
8.90E-04,
9.32E-04,
9.77E-04,
7.55E-04,
7.85E-04,
8.17E-04,
8.52E-04,
8.90E-04,
9.32E-04,
9.77E-04,
7.55E-04,
7.85E-04,
8.17E-04,
8.52E-04,
8.90E-04,
9.32E-04,
9.77E-04
]

z = [4.077424497,
4.358253892, 
4.610475878, 
4.881769469,
5.153063061,
5.323277142,
5.462023074,
4.610475878,
4.840765517,
5.04864602,
5.235070966,
5.351407761,
5.440090728,
5.540693448,
4.960439843,
5.118257381,
5.266539115,
5.370479367,
5.440090728,
5.528296904,
5.5816974,
]

popt, pcov = curve_fit(F0, (x, y), z, method = 'lm')
print(popt)

popt, pcov = curve_fit(F1, (x, y), z, method = 'lm')
print(popt)

输出是：

[1.34957781e+00 1.05456428e-01 1.00000000e+00 4.14879613e+04
 1.00000000e+00]
[1.34957771e+00 1.05456434e-01 4.14879603e+04]

您可以看到分段拟合中的参数值保持为初始值。我知道我没有以正确的方式做这件事。请纠正我。

Answer 1

问题的主要来源是这种方法对定义从一个函数到另一个函数的切换的变量值不敏感（有关类似解释，请参阅this 响应）。而且启动参数的选择也不好。

由于没有提供起始值，curve_fit 为所有拟合参数选择值 1（请参阅herep0 的默认值）。由于拟合算法通过对参数进行微小变化来工作，y0 在 1 左右以小步长变化，这不会导致函数输出发生变化（所有y 值都远小于 1）。因为y[i] < y0 总是True 并且只有第一个分支被评估，并且函数的输出不依赖于c 的值。这就解释了为什么 y0 和 c 保持初始值。

人们可能期望将 y0 初始值设置在所评估的值范围内（即大约 8E-4）可能会解决问题。事实上，由于第二个分支被评估，c 的值现在已经优化。尽管如此，y0 值将保持不变。由于拟合算法可以测试值的非常小的变化，因此这些变化不足以从两个实验性y 值之间的间隔移动到另一个。在这种特殊情况下，如果选择 8E-4，微小的变化将永远不足以使其超过 8.17E-04 或低于 7.85E-4，即包含初始 y0 选择的值。

通常可以绕过这个问题，使函数显式依赖y0 的值。一个明智的选择是重新定义函数，使y0 处的值无论采用哪个分支都是相同的（即确保函数是连续的）。在这种情况下，函数定义并不能确保如此。一个合理的改变是：

def F2(X, a, b, c, c0, y0):
    x, y = X
    value = []
    for i in range(0, len(x)):
        lnZ = x[i] + c0 * y[i]
        if y[i] >= y0:
            lnZ += c * (y[i]-y0)
        val = a + (b*lnZ)
        value.append(val)
    return value

这会改变参数c 的含义，并将结果限制为仅连续函数。在这种情况下，y0 的值确实是函数的转折点。然而，它产生了预期的结果：

popt2, pcov = curve_fit(F2, (x, y), z, p0=(1, 1, 1E4, 1E4, 9.1E-4), method = 'lm')
print(popt2)

结果：

[-1.93417968e-01  1.05456433e-01 -3.65740192e+04  5.97890809e+04
  8.64354057e-04]

更好的 (pythonic) 函数定义避免了 for 循环：

def F3(X, a, b, c, c0, y0):
    x, y = X
    lnZ = x + c0 * y
    idx = np.where(y>=y0)
    lnZ[idx] += c * (y[idx] - y0)
    rv = a + (b * lnZ)
    return rv

对于较大的数据集，这可能会快得多。