Scipy的curve_fit无法拟合tophat函数。

Question

我正在尝试将一个顶帽函数拟合到一些数据上，即 f(x) 在整个实数线上都是常数，除了一个有限长度的区间，它等于另一个常数。我的参数是顶帽函数的两个常数、中点和宽度，我正在尝试使用scipy.optimize.curve_fit来得到所有这些参数。不幸的是，curve_fit在获取帽子的宽度时遇到了麻烦。无论我做什么，它都拒绝测试除我开始的宽度之外的任何值，并且其余数据的拟合非常糟糕。以下代码片段说明了这个问题：

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

def tophat(x, base_level, hat_level, hat_mid, hat_width):
    ret=[]
    for xx in x:
        if hat_mid-hat_width/2. < xx < hat_mid+hat_width/2.:
            ret.append(hat_level)
        else:
            ret.append(base_level)
    return np.array(ret)

x = np.arange(-10., 10., 0.01)
y = tophat(x, 1.0, 5.0, 0.0, 1.0)+np.random.rand(len(x))*0.2-0.1

guesses = [ [1.0, 5.0, 0.0, 1.0],
            [1.0, 5.0, 0.0, 0.1],
            [1.0, 5.0, 0.0, 2.0] ]

plt.plot(x,y)

for guess in guesses:
    popt, pcov = curve_fit( tophat, x, y, p0=guess )
    print popt
    plt.plot( x, tophat(x, popt[0], popt[1], popt[2], popt[3]) )

plt.show()

为什么curve_fit在处理这个问题时表现如此糟糕，我该如何修复它？

- TheBigH

---

顶帽导数不连续，因此我使用 sp.curve_fit 可能会失败。 您可以定义一个残差函数，作为您的顶帽和数据之间的差异，并使用类似于 scipy.optimize.brute 的方法将其最小化。 - Brenlla

Answer 1

首先，tophat 的定义可以使用 numpy.where 而不是循环：

def tophat(x, base_level, hat_level, hat_mid, hat_width):
    return np.where((hat_mid-hat_width/2. < x) & (x < hat_mid+hat_width/2.), hat_level, base_level)

其次，棘手的不连续目标函数阻碍了`curve_fit`调用的优化算法。Nelder-Mead方法通常适用于粗糙的函数，但看起来`curve_fit`不能使用它。因此，我设置了一个目标函数（仅为偏差的绝对值之和）并最小化它：

def objective(params, x, y):
    return np.sum(np.abs(tophat(x, *params) - y))

plt.plot(x,y)

for guess in guesses:
    res = minimize(objective, guess, args=(x, y), method='Nelder-Mead')
    print(res.x)
    plt.plot(x, tophat(x, *(res.x)))

结果更好了，因为开始使用宽度为2的帽子，它会缩小到正确大小（请参见三次猜测中的最后一次）。

[9.96041297e-01 5.00035502e+00 2.39462103e-04 9.99759984e-01]
[ 1.00115808e+00  4.94088711e+00 -2.21340843e-05  1.04924153e-01]
[9.95947108e-01 4.99871040e+00 1.26575116e-03 9.97908018e-01]

不幸的是，当起始猜测是一个太窄的帽子时，优化器仍然会卡住。

您可以尝试其他的优化方法/目标函数组合，但我还没有找到一种可以使帽子可靠扩大的方法。

一个值得尝试的方法是不要使用那些离真实水平太近的参数；有时这可能会产生负面影响。

guesses = [ [1.0, 1.0, 0.0, 1.0],
            [1.0, 1.0, 0.0, 0.1],
            [1.0, 1.0, 0.0, 2.0] ]

我曾经设法得到

[ 1.00131181  4.99156649 -0.01109271  0.96822019]
[ 1.00137925  4.97879423 -0.05091561  1.096166  ]
[ 1.00130568  4.98679988 -0.01133717  0.99339777]

这是正确的三种宽度。然而，这只是在几次尝试中的一部分（优化过程的初始化存在一定的随机性）。一些具有相同初始点的其他尝试失败了；该过程不够健壮。

- user6655984

Answer 2

由于非线性最小二乘拟合（例如使用curve_fit()）的本质，它处理实数浮点数时表现良好，但在处理离散变量时效果不佳。在拟合过程中，每个变量都会发生微小变化（例如1e-7级别），并且该微小变化对拟合结果的影响被用来决定如何改变该变量以改善拟合。对于离散采样数据，你的hat_mid和/或hat_width的微小变化可能比数据点的间距还要小，因此对拟合没有任何影响。这就是为什么curve_fit在这个问题上表现“极其糟糕”的原因。

你可能会发现，给步骤一个有限的宽度（与离散数据的步长相当）有助于更好地找到帽子边缘的位置。

- M Newville

Answer 3

你也可以尝试拟合以下公式：

f = A0 (-Erf[A4*(u - A1) - A3] + Erf[A4*(u + A1) - A3]) + A2

其中，A0 与阶跃高度成比例，A1 与阶跃宽度成比例，A2 是基线的垂直偏移量，A3 应该是阶跃中心点距离 0 的水平距离，阶跃的斜率与 A4 成比例。

- Steve Boege