python - 使用 scipy Optimize curve_fit 拟合步骤位置变化的步骤函数

Question

我正在尝试拟合看起来像这样的 x y 数据

x = np.linspace(-2, 2, 1000)
a = 0.5
yl = np.ones_like(x[x < a]) * -0.4 + np.random.normal(0, 0.05, x[x < a].shape[0])
yr = np.ones_like(x[x >= a]) * 0.4 + np.random.normal(0, 0.05, x[x >= a].shape[0])
y = np.concatenate((yl, yr))
plt.scatter(x, y, s=2, color='k')

我正在使用 Heaviside 阶跃函数的变体

def f(x, a, b): return 0.5 * b * (np.sign(x - a))

并拟合

popt, pcov = curve_fit(f, x, y, p0=p)

其中 p 是一些初始猜测。对于任何 p curve_fit 只适合 b 而不是 a例如:

popt, pcov = curve_fit(f, x, y, p0=[-1.0, 0])我们得到 popt 是 [-1., 0.20117665]

popt, pcov = curve_fit(f, x, y, p0=[.5, 2])我们得到 popt 是 [.5, 0.79902]

popt, pcov = curve_fit(f, x, y, p0=[1.5, -2])我们得到 taht popt 是 [1.5, 0.40128229]

为什么 curve_fit 不拟合 a？

Answer 1

正如其他人所提到的，curve_fit(以及scipy.optimize中的所有其他求解器)可以很好地优化连续变量，但不能优化离散变量。它们的工作原理都是对参数值进行小的(例如，在 1.e-7 级别)更改，并查看结果中发生了什么(如果有)更改，并使用该更改来细化这些值，直到最小残差为成立。使用 np.sign 的模型函数:

def f(x, a, b): return 0.5 * b * (np.sign(x - a))

a 值的如此小的变化根本不会改变模型或拟合结果。也就是说，首先拟合将尝试起始值，例如 a=-1.0 或 a=0.5，然后尝试 a=-0.999999995 或 a=0.500000005。对于 np.sign(x-a)，它们都会给出相同的结果。拟合不知道需要将 a 更改 1 才能对结果产生任何影响。它不可能知道这一点。 np.sign() 和 np.sin() 相差一个字母，但在这方面的行为却截然不同。

真实数据采取一步是很常见的，但采样足够精细，以便该一步不会完全在一步中发生。在这种情况下，您将能够使用各种函数形式(线性斜坡、误差函数、反​​正切、逻辑等)对步骤进行建模。 @JamesPhilipps 的彻底回答给出了一种方法。我可能会使用 lmfit (作为其主要作者之一)并愿意通过查看数据来猜测参数的起始值，也许:

import numpy as np

x = np.linspace(-2, 2, 1000)
a = 0.5
yl = np.ones_like(x[x < a]) * -0.4 + np.random.normal(0, 0.05, x[x < a].shape[0])
yr = np.ones_like(x[x >= a]) * 0.4 + np.random.normal(0, 0.05, x[x >= a].shape[0])
y = np.concatenate((yl, yr))

from lmfit.models import StepModel, ConstantModel

model = StepModel() + ConstantModel()
params = model.make_params(center=0, sigma=1, amplitude=1., c=-0.5)

result = model.fit(y, params, x=x)

print(result.fit_report())

import matplotlib.pyplot as plt
plt.scatter(x, y, label='data')
plt.plot(x, result.best_fit, marker='o', color='r', label='fit')
plt.show()

这将提供良好的拟合并打印出结果

[[Model]]
    (Model(step, form='linear') + Model(constant))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 50
    # data points      = 1000
    # variables        = 4
    chi-square         = 2.32729556
    reduced chi-square = 0.00233664
    Akaike info crit   = -6055.04839
    Bayesian info crit = -6035.41737
##  Warning: uncertainties could not be estimated:
[[Variables]]
    amplitude:  0.80013762 (init = 1)
    center:     0.50083312 (init = 0)
    sigma:      4.6009e-04 (init = 1)
    c:         -0.40006255 (init = -0.5)

请注意，它会找到步骤的中心，因为它假设步骤有一些有限的宽度(sigma)，但随后发现宽度更小比 x 中的步长大。但还要注意，它无法计算参数中的不确定性，因为如上所述，解附近的 center(您的 a)的微小变化不会改变最终的拟合。 FWIW StepModel 可以使用线性、误差函数、反​​正切或逻辑作为阶跃函数。

如果您构建的测试数据具有较小的步长宽度，请说类似的东西

from scipy.special import erf    
y = 0.638  * erf((x-0.574)/0.005)  + np.random.normal(0, 0.05, len(x))

然后拟合就能够找到最佳解决方案并评估不确定性。

我希望这能解释为什么与模型函数的拟合无法细化 a 的值，以及可以采取哪些措施。