拟合阶跃函数

Question

我正在尝试使用scipy.optimize.leastsq来拟合一个阶跃函数。考虑以下例子：

import numpy as np
from scipy.optimize import leastsq

def fitfunc(p, x):
    y = np.zeros(x.shape)
    y[x < p[0]] = p[1]
    y[p[0] < x] = p[2]
    return y

errfunc = lambda p, x, y: fitfunc(p, x) - y # Distance to the target function

x = np.arange(1000)
y = np.random.random(1000)

y[x < 250.] -= 10

p0 = [500.,0.,0.]
p1, success = leastsq(errfunc, p0, args=(x, y))

print p1

这里的参数是阶跃的位置和两侧的水平值。奇怪的是，第一个自由参数一直没有变化，如果你运行这个，scipy会给出

[  5.00000000e+02  -4.49410173e+00   4.88624449e-01]

当第一个参数设置为250，第二个参数设置为-10时，它是最优的。

有没有人能告诉我为什么这可能不工作，以及如何让它正常工作呢？

如果我运行

print np.sum(errfunc(p1, x, y)**2.)
print np.sum(errfunc([250.,-10.,0.], x, y)**2.)

我发现：

12547.1054663
320.679545235

第一个数字是leastsq找到的值，第二个是它应该找到的实际最优函数的值。

Answer 1

我建议用一种方法来近似阶跃函数。与其在“变化点”处有无限大的斜率，不如在一个x的距离内让它呈线性变化（在这个例子中是1.0）。比如，如果函数的x参数xp被定义为这条线的中点，那么在xp-0.5时的值就是较低的y值，而在xp+0.5时的值就是较高的y值。在区间[xp-0.5; xp+0.5]内，函数的中间值就是这两个点之间的线性插值。

如果可以假设阶跃函数（或它的近似值）是从一个较低的值变到一个较高的值，那么我认为最后两个参数的初始猜测应该分别是最低的y值和最高的y值，而不是0.0和0.0。

---

我有两个更正：

1) np.random.random()返回的随机数范围是0.0到1.0。因此，平均值是+0.5，这也是第三个参数的值（而不是0.0）。第二个参数则是-9.5（+0.5 - 10.0），而不是-10.0。

所以

print np.sum(errfunc([250.,-10.,0.], x, y)**2.)

应该是

print np.sum(errfunc([250.,-9.5,0.5], x, y)**2.)

2) 在原始的fitfunc()中，如果x恰好等于p[0]，那么y的一个值会变成0.0。因此在这种情况下，它就不是一个阶跃函数了（更像是两个阶跃函数的和）。例如，当第一个参数的起始值是500时，就会发生这种情况。

Answer 2

我觉得最小二乘法拟合并不是处理阶跃函数的好方法。我认为它不能很好地描述不连续性。面对这个问题，我不会首先想到最小二乘法。

那为什么不考虑用傅里叶级数来近似呢？虽然在不连续的地方你总是会遇到吉布斯现象，但其他部分的函数可以根据你的需求和计算机的能力来进行很好的近似。

你到底打算用这个做什么呢？提供一些背景信息可能会更有帮助。

Answer 3

结果发现，如果我在leastsq函数中加上epsfcn这个参数，拟合效果会好很多：

p1, success = leastsq(errfunc, p0, args=(x, y), epsfcn=10.)

得到的结果是

[ 248.00000146   -8.8273455     0.40818216]

我基本理解是，首个自由参数需要移动的距离要大于相邻点之间的间距，这样才能影响残差的平方。而epsfcn这个参数与寻找梯度时使用的步长大小有关，或者说是类似的东西。