实践中的最小二乘法答案

【问题标题】：Least Squares method in practice实践中的最小二乘法
【发布时间】：2017-11-19 03:40:15
【问题描述】：

非常简单的回归任务。我有三个变量x1, x2, x3 带有一些随机噪声。我知道目标方程：y = q1*x1 + q2*x2 + q3*x3。现在我想找到目标系数：q1, q2, q3 评估使用平均相对平方误差 (RSE) (Prediction/Real - 1)^2 来评估我们的预测方法的性能。

在研究中，我发现这是普通的最小二乘问题。但是我无法从互联网上的示例中获得如何在 Python 中解决这个特定问题。假设我有数据：

import numpy as np

sourceData = np.random.rand(1000, 3)
koefs = np.array([1, 2, 3])
target = np.dot(sourceData, koefs)

（在现实生活中，数据是嘈杂的，不是正态分布。）如何在 python 中使用最小二乘法找到这个 koefs？ 任何 lib 使用情况。

【问题讨论】：

不，这不是普通的最小二乘问题。在普通最小二乘法中，您可以最小化平方误差之和。你的损失函数是不同的。对于平均绝对偏差，有quantile regression 但我不确定当你取百分比时估计器的行为是否相同。你的问题似乎更适合Cross Validated（统计部分 - 也许不是图书馆建议部分）。
@ayhan researchgate.net/publication/… 在多元线性 (ML) 模型部分重新制定后，这里有一些证据。 “这是一个可以通过奇异值分解解决的普通最小二乘 (OLS) 问题”
损失函数不是那里的平均绝对百分比误差。
你是对的，我已将平均绝对百分比误差更改为相对平方误差

标签： python numpy machine-learning scipy regression

【解决方案1】：

@ayhan 提出了宝贵意见。

您的代码存在问题：实际上您收集的数据中没有噪音。输入数据有噪声，但在相乘之后，您不会添加任何额外的噪声。

我在您的测量中添加了一些噪音，并使用最小二乘公式来拟合参数，这是我的代码：

data = np.random.rand(1000,3)

true_theta = np.array([1,2,3])
true_measurements = np.dot(data, true_theta)

noise = np.random.rand(1000) * 1

noisy_measurements = true_measurements + noise

estimated_theta = np.linalg.inv(data.T @ data) @ data.T @ noisy_measurements

estimated_theta 将接近于true_theta。如果您不在测量中添加噪声，它们将是相等的。

我使用了 python3 矩阵乘法语法。你可以用np.dot 代替@

这会使代码变长，所以我将公式拆分：

MTM_inv = np.linalg.inv(np.dot(data.T, data))
MTy = np.dot(data.T, noisy_measurements)
estimated_theta = np.dot(MTM_inv, MTy)

你可以在这里阅读最小二乘：https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#The_general_problem

更新：

或者你可以使用builtin least squares函数：

np.linalg.lstsq(data, noisy_measurements)

【讨论】：

非常酷，谢谢.. 和 numpy.linalg.lstsq 一样吗？不幸的是，我无法提供损失来处理数据中的异常值
很好，我不知道那个功能。看来您的评论将是完美的答案；）-是的，该功能是相同的。可能会更好。我不知道他们的内部实现

【解决方案2】：

除了@lhk 的答案，我还发现了很棒的 scipy Least Squares function。用它很容易得到请求的行为。

这样我们可以提供一个自定义函数来返回残差并形成相对平方误差而不是绝对平方差：

import numpy as np
from scipy.optimize import least_squares
data = np.random.rand(1000,3)

true_theta = np.array([1,2,3])
true_measurements = np.dot(data, true_theta)

noise = np.random.rand(1000) * 1

noisy_measurements = true_measurements + noise
#noisy_measurements[-1] = data[-1]  @ (1000 * true_theta) - uncoment this outliner to see how much Relative Squared Error esimator works better then default abs diff for this case.


def my_func(params, x, y):
     res = (x @ params) / y - 1 # If we change this line to: (x @ params) - y - we will got the same result as np.linalg.lstsq
     return res

res = least_squares(my_func, x0,  args=(data, noisy_measurements) ) 
estimated_theta = res.x

此外，我们可以使用 loss 参数函数提供自定义损失，该函数将处理残差并形成最终损失。

【讨论】：