scipy曲线拟合无法拟合tophat函数

Question

我试图用一个top-hat函数来拟合某些数据，即f（x）对于整个实线是常数，除了一个有限长度的线段等于另一个常数。我的参数是tophat函数的两个常数，中点和宽度，我尝试使用scipy.optimize.curve U拟合为了得到这些。不幸的是，曲线拟合无法获得帽子的宽度。不管我做什么，它都拒绝测试宽度的任何值，除了我开始的那个值之外，它与其他数据的拟合度非常差。以下代码片段说明了该问题：

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()

为什么曲线拟合在正确处理这一点上非常糟糕，我如何才能修复它？在

Answer 1

就其性质而言，非线性最小二乘拟合（如curve_fit()）适用于实数浮点数，不擅长处理离散变量。在拟合过程中，对每个变量进行小的更改（例如，在1e-7级别），并使用该小更改对拟合结果的影响来决定如何更改该变量以改进拟合。对于离散采样的数据，hat_mid和/或{}的微小变化很容易小于数据点的间距，因此对拟合没有任何影响。这就是为什么{}在这个问题上“极其糟糕”。在

您可能会发现，给这些步骤指定一个有限的宽度（也就是说，与离散数据的步长相当）有助于更好地找到边缘所在。在

Answer 2

首先，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不能使用它。所以我建立了一个目标函数（偏差的绝对值之和）并最小化：

^{pr2}$

结果更好，因为从宽度为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]

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

你可以尝试其他的优化方法/目标函数组合，但我还没有找到一个能使hat可靠扩展的方法。在

可以尝试的一件事是不要使用接近真实水平的太的参数；这有时可能会造成伤害。与

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]

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