我想通过最小化中位数绝对误差来执行一维线性回归。
最初假设它应该是一个相当标准的用例，但快速搜索令人惊讶地发现，所有回归和内插函数都使用均方误差。
因此，我的问题是：是否存在可以对一维进行基于中值误差的线性回归的函数？
请您参考如下方法:
正如评论中已经指出的那样，即使您自己要求的是定义明确的，但解决方案的正确方法也将取决于模型的属性。让我们看看为什么，让我们看看通才优化方法能使您走多远，让我们看看一些数学如何简化问题。底部包含可复制复制的解决方案。
首先，在应用专门算法的意义上，最小二乘拟合比您尝试做的“容易”。例如，SciPy的leastsq使用Levenberg--Marquardt algorithm假定您的optimization objective是平方和。当然，在线性回归的特殊情况下，问题也可能是solved analytically。
除了实际优势之外，最小二乘线性回归理论上也可以证明是合理的：如果观测值的residuals是independent和normally distributed（如果您发现central limit theorem适用于模型，则可以证明是合理的），那么模型参数的maximum likelihood estimate将是通过最小二乘法获得的参数。同样，最小化mean absolute error优化目标的参数将是Laplace distributed残差的最大似然估计。现在，如果您事先知道数据非常脏，以至于对残差的正态性的假设将失败，那么您要尝试的操作将比普通最小二乘法具有优势，但是即使那样，您仍然可以证明会影响到残差正态性的其他假设。选择目标函数，所以我很好奇您在这种情况下最终会如何？
使用数值方法
顺便说一句，一些一般性的评论适用。首先，请注意，SciPy带有large selection of general purpose algorithms，您可以直接将其应用于您的案例。例如，让我们看看如何在单变量情况下应用minimize。
# Generate some data
np.random.seed(0)
n = 200
xs = np.arange(n)
ys = 2*xs + 3 + np.random.normal(0, 30, n)
# Define the optimization objective
def f(theta):
return np.median(np.abs(theta[1]*xs + theta[0] - ys))
# Provide a poor, but not terrible, initial guess to challenge SciPy a bit
initial_theta = [10, 5]
res = minimize(f, initial_theta)
# Plot the results
plt.scatter(xs, ys, s=1)
plt.plot(res.x[1]*xs + res.x[0])
这样肯定会更糟。正如@sascha在评论中指出的那样，目标的不平滑很快成为一个问题，但是，再次取决于模型的模样，您可能会发现自己正在寻找的 convex足以节省您的钱。
如果您的参数空间是低维的，则只需绘制优化范围即可直观地了解优化的稳定性。
theta0s = np.linspace(-100, 100, 200)
theta1s = np.linspace(-5, 5, 200)
costs = [[f([theta0, theta1]) for theta0 in theta0s] for theta1 in theta1s]
plt.contour(theta0s, theta1s, costs, 50)
plt.xlabel('$\\theta_0$')
plt.ylabel('$\\theta_1$')
plt.colorbar()
在上面的特定示例中，如果初始猜测无效，则通用优化算法将失败。
initial_theta = [10, 10000]
res = minimize(f, initial_theta)
plt.scatter(xs, ys, s=1)
plt.plot(res.x[1]*xs + res.x[0])
还请注意，SciPy的许多算法都受益于目标的 Jacobian，即使您的目标是不可区分的，再次取决于您要优化的内容，残差也很可能会因此而产生，您的目标可能是可区分的 almost everywhere，因为您能够提供导数（例如，中位数的导数变成值为中位数的函数的导数）。
在我们的例子中，提供雅可比行列似乎没有什么特别的帮助，如以下示例所示；在这里，我们增加了残差的方差，足以使整个事物崩溃。
np.random.seed(0)
n = 201
xs = np.arange(n)
ys = 2*xs + 3 + np.random.normal(0, 50, n)
initial_theta = [10, 5]
res = minimize(f, initial_theta)
plt.scatter(xs, ys, s=1)
plt.plot(res.x[1]*xs + res.x[0])
def fder(theta):
"""Calculates the gradient of f."""
residuals = theta[1]*xs + theta[0] - ys
absresiduals = np.abs(residuals)
# Note that np.median potentially interpolates, in which case the np.where below
# would be empty. Luckily, we chose n to be odd.
argmedian = np.where(absresiduals == np.median(absresiduals))[0][0]
residual = residuals[argmedian]
sign = np.sign(residual)
return np.array([sign, sign * xs[argmedian]])
res = minimize(f, initial_theta, jac=fder)
plt.scatter(xs, ys, s=1)
plt.plot(res.x[1]*xs + res.x[0])
在此示例中，我们发现自己陷入了奇异之处。
theta = res.x
delta = 0.01
theta0s = np.linspace(theta[0]-delta, theta[0]+delta, 200)
theta1s = np.linspace(theta[1]-delta, theta[1]+delta, 200)
costs = [[f([theta0, theta1]) for theta0 in theta0s] for theta1 in theta1s]
plt.contour(theta0s, theta1s, costs, 100)
plt.xlabel('$\\theta_0$')
plt.ylabel('$\\theta_1$')
plt.colorbar()
而且，这是您在最低限度附近会发现的一团糟：
theta0s = np.linspace(-20, 30, 300)
theta1s = np.linspace(1, 3, 300)
costs = [[f([theta0, theta1]) for theta0 in theta0s] for theta1 in theta1s]
plt.contour(theta0s, theta1s, costs, 50)
plt.xlabel('$\\theta_0$')
plt.ylabel('$\\theta_1$')
plt.colorbar()
如果您发现自己在这里，则可能需要使用其他方法。如@sascha所述，仍然使用通用优化方法的示例包括用更简单的方法替换目标。另一个简单的示例是使用各种不同的初始输入来运行优化：
min_f = float('inf')
for _ in range(100):
initial_theta = np.random.uniform(-10, 10, 2)
res = minimize(f, initial_theta, jac=fder)
if res.fun < min_f:
min_f = res.fun
theta = res.x
plt.scatter(xs, ys, s=1)
plt.plot(theta[1]*xs + theta[0])
部分分析方法
请注意，将
theta最小化
f的值还将最小化残差平方的中值。搜索“最小中位数平方”可能会为您提供有关此特定问题的更多相关信息。
在这里，我们遵循 Rousseeuw -- Least Median of Squares Regression，其第二部分包括一种算法，该算法可将上述二维优化问题简化为可能更易于解决的一维问题。如上所述，假设我们有奇数个数据点，因此我们不必担心中位数定义的歧义。
首先要注意的是，如果您只有一个变量（对问题的二读实际上可能是您感兴趣的情况），那么很容易证明以下函数可以提供最小的分析能力。 。
def least_median_abs_1d(x: np.ndarray):
X = np.sort(x) # For performance, precompute this one.
h = len(X)//2
diffs = X[h:] - X[:h+1]
min_i = np.argmin(diffs)
return diffs[min_i]/2 + X[min_i]
现在，诀窍在于，对于固定的
theta1，通过将以上内容应用于
theta0可获得
f(theta0, theta1)最小化
ys - theta0*xs的值。换句话说，我们已将问题简化为单个变量的函数（以下称为
g）的最小化。
def best_theta0(theta1):
# Here we use the data points defined above
rs = ys - theta1*xs
return least_median_abs_1d(rs)
def g(theta1):
return f([best_theta0(theta1), theta1])
尽管这可能比上面的二维优化问题更容易受到攻击，但是我们还没有完全脱离森林，因为此新功能具有其自身的局部最小值：
theta1s = np.linspace(0, 3, 500)
plt.plot(theta1s, [g(theta1) for theta1 in theta1s])
theta1s = np.linspace(1.5, 2.5, 500)
plt.plot(theta1s, [g(theta1) for theta1 in theta1s])
在我有限的测试中，
basinhopping似乎能够一致地确定最小值。
from scipy.optimize import basinhopping
res = basinhopping(g, -10)
print(res.x) # prints [ 1.72529806]
此时，我们可以将所有内容包装起来，并检查结果看起来是否合理：
def least_median(xs, ys, guess_theta1):
def least_median_abs_1d(x: np.ndarray):
X = np.sort(x)
h = len(X)//2
diffs = X[h:] - X[:h+1]
min_i = np.argmin(diffs)
return diffs[min_i]/2 + X[min_i]
def best_median(theta1):
rs = ys - theta1*xs
theta0 = least_median_abs_1d(rs)
return np.median(np.abs(rs - theta0))
res = basinhopping(best_median, guess_theta1)
theta1 = res.x[0]
theta0 = least_median_abs_1d(ys - theta1*xs)
return np.array([theta0, theta1]), res.fun
theta, med = least_median(xs, ys, 10)
# Use different colors for the sets of points within and outside the median error
active = ((ys < theta[1]*xs + theta[0] + med) & (ys > theta[1]*xs + theta[0] - med))
not_active = np.logical_not(active)
plt.plot(xs[not_active], ys[not_active], 'g.')
plt.plot(xs[active], ys[active], 'r.')
plt.plot(xs, theta[1]*xs + theta[0], 'b')
plt.plot(xs, theta[1]*xs + theta[0] + med, 'b--')
plt.plot(xs, theta[1]*xs + theta[0] - med, 'b--')
