python - 拟合阶跃函数

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

其中第一个数字是最小平方找到的，第二个是它应该找到的实际最优函数的值。

Answer 1

事实证明，如果我将 epsfcn= 参数添加到 minimumsq，拟合效果会更好：

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

结果是

[ 248.00000146   -8.8273455     0.40818216]

我的基本理解是，第一个自由参数必须移动得超过相邻点之间的间距才能影响残差的平方，而 epsfcn 与使用多大的步骤来找到梯度或类似的东西有关。

Answer 2

我建议逼近阶跃函数。而不是在“变化点”处的无限斜率使其线性超过一个 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 3

我不认为最小二乘拟合是提出一步近似值的方法。我不相信它会给你一个令人满意的不连续性描述。在解决这个问题时，我不会首先想到最小二乘法。

为什么不使用傅里叶级数近似呢？在不连续处，您总是会遇到 Gibbs 现象，但函数的其余部分可以近似为您和您的 CPU 可以承受的范围。

你到底打算用这个做什么？一些上下文可能会有所帮助。

Answer 4

很可能您的优化陷入了局部最小值。我不知道 leastsq 真正的工作原理是什么，但如果你给它一个初始估计 (0, 0, 0)，它也会卡在那里。

您可以在数字上检查初始估计的梯度（在 +/- epsilon 处评估非常小的 epsilon 并将 bei 2*epsilon 除以差值），我敢打赌它会在 0 左右。

Answer 5

使用 statsmodel ols。ols 使用普通最小二乘法进行曲线拟合