使用平均绝对偏差定制回归答案

【问题标题】：custom made regression using average absolute deviation使用平均绝对偏差定制回归
【发布时间】：2020-06-08 16:51:18
【问题描述】：

在this post 之后，我现在严重怀疑R-squared 或F-test 是否很好地表明了对某些具有随机噪声的数据的良好线性拟合。因此，我想开发一个定制的回归函数，这样我既可以了解它是如何工作的，也可以改进现有的工具。

考虑这些随机生成的 ndarrays x 和 y：

import numpy as np

np.random.seed(42)

x = np.random.rand(30) * 10
y = 1.5 * x + 0.3 + (np.random.rand(30) - 0.5) * 3.5

现在我可以定义任何一组数据点的平均/平均绝对偏差：

def aad(X, Y, a, b): # assumes X and Y are of the identical shape/size
    n = X.size # highly unsafe!
    U = (a * X + Y - b) / 2 / a
    V = (a * X + Y + b) / 2
    E = np.sqrt(np.power((X - U), 2) + np.power((Y - V), 2))
    return E.sum() / n

在我看来，这是将y = a * x + b 线的适合度量化为一对数据点的最佳方式。该函数简单地找到假设线到任何数据点的最近点，然后计算该点与线之间的垂直距离。

现在我需要一个函数，比如说：

linearFit(X, Y)

给出X 和Y 的形状相同的ndarray，找到a 和b，这使得aad(X, Y, a, b) 最小。重要的是，结果必须是绝对最小值，而不仅仅是局部结果。

当然，本着 SO 最佳实践的精神，我已经尝试了 scipy.optimize 函数 fmin 和 brute，正如您在 above-mentioned post 和 here 中看到的那样。但是，我似乎无法理解这些函数的正确语法。如果您能帮我找到假定的linearFit 函数的规范和高性能实现，我将不胜感激。提前感谢您的支持。

P.S.提供了一个临时解决方法here：

from scipy.optimize import minimize

aad_ = lambda P: aad(P[0], P[1], x1, y1)
minimize(aad_, x0=[X0, Y0])

但是，我得到的结果并不那么有希望！求解器不成功，我收到消息：

由于精度损失，不一定能达到预期误差

【问题讨论】：

计算最小绝对偏差回归的方法有很多，google一些算法就可以了。虽然这是一个迭代问题，但与 ols 相比有一些缺点
@bryan60 好的，这就是它的名称。我不知道。我也不知道常用的方法叫做Ordinary least squares (OLS)回归。谢谢。 scipy.optimize 函数有什么简洁的方法吗？
我真的从未遇到过需要 L1 回归的用例。去除异常值会给你大致相同的结果。你不太适合随机数据。这就是随机数据的意义所在。任何相关性都是虚假的。
@bryan60 编辑了帖子。它不是真正的随机数据，而是带有随机噪声的数据。
如果噪声是随机的，那么绝对误差与平方误差不会有太大区别。当异常值的权重较小时，绝对值最强。如果您可以测量噪声的原因并将其用于多元回归，您将获得更好的拟合。或者，如果噪声确实像信号噪声一样是随机的，还有其他技术。

标签： python numpy scipy regression scipy-optimize

【解决方案1】：

首先，感谢this post，我意识到这不是上面 cmets 中讨论的普通最小二乘 (OLS) 回归。它实际上有许多名称，其中包括戴明回归、正交距离回归 (ODR) 和总最小二乘法 (TLS)。还有，of course，a Python packagescipy.odr 也是如此！它的语法有点奇怪，文档帮助不大，但是可以在here 找到一个很好的教程。

接下来我在aad 定义中发现了一个小错误，并将其重命名并修复为：

def aaod(a, b, X, Y): # assumes X and Y are of the identical shape/size
    n = X.size # still highly unsafe! don't use it in real production
    U = (a * X + Y - b) / 2 / a
    V = (a * X + Y + b) / 2
    E = np.sqrt(np.power((X - U), 2) + np.power((Y - V), 2))
    return E.sum() / n

代表平均绝对正交距离。现在将我们的拟合函数定义为：

from scipy.optimize import minimize
from scipy.stats import linregress

def odrFit(X, Y):
    X0 = linregress(X, Y) # wait this is cheating!
    aaod_ = lambda P: aaod(P[0], P[1], X, Y)
    res = minimize(aaod_, x0=X0[:2], method = 'Nelder-Mead')
    res_list = res.x.tolist()
    res_list.append(aaod_(res_list))
    return res_list

这不一定是最高效和最规范的实现。我从 here 和 method = 'Nelder-Mead' 从 here 学到的临时 lambda 函数的解决方法。 scipy.odr 的实现也可以这样完成：

from scipy.odr import Model, ODR, RealData

def f(B, x):
    return B[0]*x + B[1]

linear = Model(f)
mydata = RealData(x, y)
myodr = ODR(mydata, linear, beta0=[1., 2.])
myoutput = myodr.run()

现在比较我们定制的odrFit()函数和scipy.stats.linregress()的结果：

slope, intercept, r_value, p_value, std_err = linregress(x,y)

print(*odrFit(x, y)) 
# --> 1.4804181575739097, 0.47304584702448255, 0.6008218016339527

print(slope, intercept, aaod(slope, intercept, x, y))
# --> 1.434483032725671 0.5747705643012724 0.608021569291401

print(*myoutput.beta, aaod(*myoutput.beta, x, y))
# --> 1.5118079563432785 0.23562547897245803 0.6055838996104704

这表明我们的函数实际上比 Scipy 的最小绝对偏差回归方法更准确。这实际上可能只是纯粹的运气，需要进行更多测试才能得出可靠的结论。完整代码见here。

【讨论】：