Основанная на медиане линейная регрессия в Python

Question

Я хотел бы выполнить одномерную линейную регрессию путем минимизации средней абсолютной ошибки.

Первоначально предполагая, что это должен быть достаточно стандартный вариант использования, быстрый поиск неожиданно выявил, что все функции регрессии и интерполяции используют среднеквадратичную ошибку.

Поэтому мой вопрос: есть ли функция, которая может выполнять линейную регрессию на основе медианной ошибки для одного измерения?

Answer 1

Как уже указывалось в комментариях, даже если то, что вы запрашиваете, само по себе четко определено, правильный подход к его решению будет зависеть от свойств вашей модели. Давайте посмотрим, почему, давайте посмотрим, насколько далеко вас продвигает универсальный подход к оптимизации, и посмотрим, как математика может упростить задачу. Копируемое решение включено внизу.

Прежде всего, подгонка наименьших квадратов "проще", чем то, что вы пытаетесь сделать, в том смысле, что применяются специализированные алгоритмы; например, SciPy's leastsqиспользует алгоритм Левенберга - Марквардта, который предполагает, что ваша цель оптимизации представляет собой сумму квадратов. Конечно, в частном случае линейной регрессии проблема также может быть решена аналитически.

Помимо практических преимуществ, линейная регрессия наименьших квадратов также может быть теоретически обоснована: если остатки ваших наблюдений независимы и нормально распределены (что вы можете обосновать, если обнаружите, что теорема о центральном пределе применима в вашей модели), тогда оценка максимального правдоподобия из ваших параметров модели будут те, которые получены через наименьших квадратов. Аналогично, параметры, минимизирующие среднюю цель оптимизации абсолютной ошибки, будут оценками максимального правдоподобия для распределенных остатков Лапласа. Теперь то, что вы пытаетесь сделать, будет иметь преимущество перед обычными наименьшими квадратами, если вы заранее знаете, что ваши данные настолько грязны, что допущения относительно нормальности невязок будут ошибочными, но даже тогда вы сможете обосновать другие предположения, которые повлияют на Выбор целевой функции, поэтому мне интересно, как вы оказались в такой ситуации?

Используя численные методы

При этом некоторые общие замечания применимы. Прежде всего, обратите внимание, что 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*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--')