заказная регрессия с использованием среднего абсолютного отклонения - PullRequest
Новая
       94

заказная регрессия с использованием среднего абсолютного отклонения

0 голосов
/ 25 февраля 2020

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

Рассмотрим эти случайно сгенерированные ndarrays x и y: 

import numpy as np

np.random.seed(42)

x = np.random.rand(30) * 10
y = 1.5 * x + 0.3 + (np.random.rand(30) - 0.5) * 3.5

Теперь я могу определить среднее / среднее абсолютное отклонение любого набора точек данных с помощью: 

def aad(X, Y, a, b): # assumes X and Y are of the identical shape/size
    n = X.size # highly unsafe!
    U = (a * X + Y - b) / 2 / a
    V = (a * X + Y + b) / 2
    E = np.sqrt(np.power((X - U), 2) + np.power((Y - V), 2))
    return E.sum() / n

, что, на мой взгляд, является наилучшим способом количественного определения соответствия линии y = a * x + b паре точек данных. Функция просто находит ближайшую точку предполагаемой линии к любой точке данных и затем вычисляет перпендикулярное расстояние между точкой и линией.

Теперь мне нужно иметь функцию, скажем: 

linearFit(X, Y)

, который с учетом ndarrays идентичной формы X и Y, находит a и b, которые составляют aad(X, Y, a, b) минимум. Важно, чтобы результат был абсолютным минимумом, а не только локальным.

Конечно, в духе лучших практик SO, я уже попробовал scipy.optimize функции fmin и brute, как вы можете видеть в вышеупомянутом посте и также здесь . Однако, похоже, я не могу понять правильный синтаксис для этих функций. Буду признателен, если вы поможете мне найти каноническую и производительную реализацию для предполагаемой функции linearFit. Спасибо за вашу поддержку заранее.

PS Временное решение проблемы здесь : 

from scipy.optimize import minimize

aad_ = lambda P: aad(P[0], P[1], x1, y1)
minimize(aad_, x0=[X0, Y0])

однако результаты, которые я получаю, не столь многообещающие! Решение не удается, и я получаю сообщение:

Требуемая ошибка не обязательно достигается из-за потери точности

1 Ответ

0 голосов
/ 28 февраля 2020

Прежде всего, благодаря этому посту Я понял, что это не обычная регрессия наименьших квадратов (OLS), как обсуждалось в комментариях выше. На самом деле его называют многими именами, среди которых регрессия Деминга, регрессия ортогонального расстояния (ODR) и наименьшие квадраты (TLS). Также есть, конечно , Python пакет scipy.odr для этого! Его синтаксис немного странный, и документация не очень помогает, но хороший учебник можно найти здесь .

Nex Я обнаружил небольшую ошибку в определении aad, переименовал и исправил ее на: 

def aaod(a, b, X, Y): # assumes X and Y are of the identical shape/size
    n = X.size # still highly unsafe! don't use it in real production
    U = (a * X + Y - b) / 2 / a
    V = (a * X + Y + b) / 2
    E = np.sqrt(np.power((X - U), 2) + np.power((Y - V), 2))
    return E.sum() / n

, обозначая среднее абсолютное ортогональное расстояние. Теперь определяем нашу функцию подбора как: 

from scipy.optimize import minimize
from scipy.stats import linregress

def odrFit(X, Y):
    X0 = linregress(X, Y) # wait this is cheating!
    aaod_ = lambda P: aaod(P[0], P[1], X, Y)
    res = minimize(aaod_, x0=X0[:2], method = 'Nelder-Mead')
    res_list = res.x.tolist()
    res_list.append(aaod_(res_list))
    return res_list

, что не обязательно является наиболее эффективной и канонической реализацией. Обходной путь с временной функцией lambda, которую я узнал из здесь и method = 'Nelder-Mead' из здесь . Реализацию scipy.odr можно также выполнить следующим образом: 

from scipy.odr import Model, ODR, RealData

def f(B, x):
    return B[0]*x + B[1]

linear = Model(f)
mydata = RealData(x, y)
myodr = ODR(mydata, linear, beta0=[1., 2.])
myoutput = myodr.run()

Теперь сравниваем результат между нашей заказной функцией odrFit() и scipy.stats.linregress(): 

slope, intercept, r_value, p_value, std_err = linregress(x,y)

print(*odrFit(x, y)) 
# --> 1.4804181575739097, 0.47304584702448255, 0.6008218016339527

print(slope, intercept, aaod(slope, intercept, x, y))
# --> 1.434483032725671 0.5747705643012724 0.608021569291401

print(*myoutput.beta, aaod(*myoutput.beta, x, y))
# --> 1.5118079563432785 0.23562547897245803 0.6055838996104704

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

...