Python中基于中值的线性回归

Question

我想通过最小化绝对误差中位数来执行一维线性回归.

虽然最初假设它应该是一个相当标准的用例,但快速搜索令人惊讶地发现所有回归和插值函数都使用均方误差.

因此我的问题是:是否有一个函数可以对一个维度执行基于中值误差的线性回归？

Answer 1

正如评论中已经指出的那样,即使您要求的内容定义明确,其解决方案的正确方法也将取决于模型的属性.让我们看看为什么,让我们看看通才优化方法能给你带来多远,让我们看看一些数学可以如何简化问题.底部包含可复制的解决方案.

首先,在专业算法适用的意义上,最小二乘拟合比你试图做的更"容易"; 例如,SciPy leastsq使用Levenberg - Marquardt算法,该算法假设您的优化目标是平方和.当然,在线性回归的特殊情况下,问题也可以通过分析解决.

除了实际优势之外,最小二乘线性回归在理论上也是合理的:如果你的观测的残差是独立的和正态分布的(如果你发现中心极限定理适用于你的模型,你可以证明这是正确的),那么最大似然估计您的模型参数将是通过最小二乘法获得的参数.类似地,最小化平均绝对误差优化目标的参数将是拉普拉斯分布式残差的最大似然估计.如果你事先知道你的数据是如此肮脏以至于对残差正常性的假设会失败,那么你现在要做的事情将优于普通的最小二乘法,但即便如此,你也可以证明其他假设会影响到选择目标函数,所以我很好奇你是如何在这种情况下结束的？

使用数值方法

除此之外,一些一般性评论适用.首先,请注意SciPy提供了大量通用算法,您可以直接在您的案例中应用.作为一个例子,让我们看看如何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在评论中指出的那样,目标的不平滑性很快成为一个问题,但是,再次取决于你的模型到底是什么样的,你可能会发现自己看到的东西足够凸显可以拯救你.

如果您的参数空间是低维的,只需绘制优化格局即可直观了解优化的稳健性.

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的许多算法都受益于提供目标的雅可比行列式,即使您的目标不可微分,再次依赖于您要优化的内容,您的残差也可能是,因此,您的因为你能够提供衍生物,所以几乎在任何地方都可以区分目标(例如,中位数的导数成为函数的导数,其值为中位数).

在我们的例子中,提供雅可比行列似乎并不特别有用,如下例所示; 在这里,我们将残差的方差提高到足以使整个事物分崩离析.

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 - 最小平方回归中位数,其第二部分包括一个算法,用于将上面的二维优化问题简化为可能更容易解决的一维优化问题.假设如上所述我们有奇数个数据点,所以我们不必担心中位数定义的模糊性.

首先要注意的是,如果您只有一个变量(在您第二次阅读您的问题时,实际上可能是您感兴趣的情况),那么很容易证明以下函数提供了最小的分析.

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*xsg

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--')