回归模型中成本函数的 L1 范数而不是 L2 范数答案

【问题标题】：L1 norm instead of L2 norm for cost function in regression model回归模型中成本函数的 L1 范数而不是 L2 范数
【发布时间】：2019-01-23 18:26:44
【问题描述】：

我想知道 Python 中是否有一个函数可以与scipy.linalg.lstsq 完成相同的工作，但使用“最小绝对偏差”回归而不是“最小二乘”回归 (OLS)。我想使用L1 规范，而不是L2 规范。

事实上，我有 3d 点，我想要它们中最适合的平面。常见的方法是使用最小二乘法，例如 Github link。但是众所周知，这并不总是提供最佳拟合，尤其是当我们的数据集中有闯入者时。最好计算最小的绝对偏差。两种方法的区别更多解释here。

它不能通过诸如 MAD 之类的函数来求解，因为它是一个 Ax = b 矩阵方程，并且需要循环来最小化结果。我想知道是否有人知道 Python 中的相关函数 - 可能在线性代数包中 - 可以计算“最小绝对偏差”回归？

【问题讨论】：

Where can I find mad (mean absolute deviation) in scipy?的可能重复
@Joel 这似乎不是链接问题的副本。虽然两者都处理 MAD，但这个问题更进一步将其用作优化目标，这不是链接问题的内容。
@N.Wouda 这两个问题都询问如何在 Python 中使用 MAD。我链接的问题如何在处理优化目标时“更进一步”？
能否请您发布数据链接？
@Joel OP 写道“如果 Python 中有一个函数可以与 scipy.linalg.lstsq 进行相同的工作，但最小化绝对偏差而不是最小二乘偏差”。这不是sm.robust.mad 所做的：它只是计算偏差，而不是优化参数。这就是为什么我觉得这个问题更进一步，因此不应该作为重复关闭。

标签： python machine-learning regression least-squares

【解决方案1】：

使用scipy.optimize.minimize 和自定义cost_function 自己滚动并不难。

让我们先进口必需品，

from scipy.optimize import minimize
import numpy as np

并定义一个自定义成本函数（以及一个用于获取拟合值的便利包装器），

def fit(X, params):
    return X.dot(params)


def cost_function(params, X, y):
    return np.sum(np.abs(y - fit(X, params)))

那么，如果你有一些X（设计矩阵）和y（观察），我们可以做以下，

output = minimize(cost_function, x0, args=(X, y))

y_hat = fit(X, output.x)

x0 是最佳参数的一些合适的初始猜测（您可以在这里接受@JamesPhillips 的建议，并使用 OLS 方法中的拟合参数）。

无论如何，当用一个有点做作的例子进行测试时，

X = np.asarray([np.ones((100,)), np.arange(0, 100)]).T
y = 10 + 5 * np.arange(0, 100) + 25 * np.random.random((100,))

我发现，

      fun: 629.4950595335436
 hess_inv: array([[  9.35213468e-03,  -1.66803210e-04],
       [ -1.66803210e-04,   1.24831279e-05]])
      jac: array([  0.00000000e+00,  -1.52587891e-05])
  message: 'Optimization terminated successfully.'
     nfev: 144
      nit: 11
     njev: 36
   status: 0
  success: True
        x: array([ 19.71326758,   5.07035192])

还有，

fig = plt.figure()
ax = plt.axes()

ax.plot(y, 'o', color='black')
ax.plot(y_hat, 'o', color='blue')

plt.show()

蓝色为拟合值，黑色为数据。

【讨论】：

你用过分位数回归function吗？我四处询问，我的前辈说“LAD 模型是分位数回归的一个特例，其中 q=0.5。”我想知道这两个函数有多大不同，它们的结果会有多大不同。
@Sara.Eft 我听说过它，但我以前从未使用过它。我怀疑你会想改用statsmodels.regression.quantile_regression.QuantReg，主要是因为它有更多的功能，使用q = 0.5（在手册中也有）。我没有意识到也可以使用分位数回归。

【解决方案2】：

您可以使用 scipy.minimize 函数解决您的问题。您必须设置要最小化的函数（在我们的例子中是 Z= aX + bY + c 形式的平面）和误差函数（L1 范数），然后使用一些起始值运行最小化器。

import numpy as np
import scipy.linalg
from scipy.optimize import minimize
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

def fit(X, params):
    # 3d Plane Z = aX + bY + c
    return X.dot(params[:2]) + params[2]

def cost_function(params, X, y):
    # L1- norm
    return np.sum(np.abs(y - fit(X, params)))

我们生成 3d 点

# Generating  3-dim points
mean = np.array([0.0,0.0,0.0])
cov = np.array([[1.0,-0.5,0.8], [-0.5,1.1,0.0], [0.8,0.0,1.0]])
data = np.random.multivariate_normal(mean, cov, 50)

最后我们运行最小化器

output = minimize(cost_function, [0.5,0.5,0.5], args=(np.c_[data[:,0], data[:,1]], data[:, 2]))
y_hat = fit(np.c_[data[:,0], data[:,1]], output.x)

X,Y = np.meshgrid(np.arange(min(data[:,0]), max(data[:,0]), 0.5),    np.arange(min(data[:,1]), max(data[:,1]), 0.5))
XX = X.flatten()
YY = Y.flatten()


# # evaluate it on grid
Z = output.x[0]*X + output.x[1]*Y + output.x[2]
fig = plt.figure(figsize=(10,10))
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.2)
ax.scatter(data[:,0], data[:,1], data[:,2], c='r')
plt.show()

注意：我使用了之前的响应代码和来自 github 的 code 作为开始

【讨论】：