Подгонка ортогональной регрессии в методе наименьших квадратов - PullRequest
Новая
       59

Подгонка ортогональной регрессии в методе наименьших квадратов

12 голосов
/ 21 февраля 2012

Метод leastsq в scipy lib соответствует кривой для некоторых данных. И этот метод подразумевает, что в этих данных значения Y зависят от некоторого аргумента X. И вычисляет минимальное расстояние между кривой и точкой данных по оси Y (dy)

Но что, если мне нужно рассчитать минимальное расстояние по обеим осям (dy и dx)

Есть ли способы реализовать этот расчет?

Вот пример кода при использовании вычисления по одной оси: 

import numpy as np
from scipy.optimize import leastsq

xData = [some data...]
yData = [some data...]

def mFunc(p, x, y):
    return y - (p[0]*x**p[1])  # is takes into account only y axis

plsq, pcov = leastsq(mFunc, [1,1], args=(xData,yData))
print plsq

Я недавно попробовал библиотеку scipy.odr, и она возвращает правильные результаты только для линейной функции. Для других функций, таких как y = a * x ^ b, он возвращает неверные результаты. Вот как я это использую: 

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

myModel = Model(f)
myData = Data(xData, yData)
myOdr = ODR(myData, myModel , beta0=[1,1])
myOdr.set_job(fit_type=0) #if set fit_type=2, returns the same as leastsq
out = myOdr.run()
out.pprint()

Возвращает неверные результаты, которые не нужны, а в некоторых входных данных даже близко к реальным. Может быть, есть какие-то особые способы его использования, что я делаю не так?

Ответы [ 3 ]

7 голосов
/ 21 февраля 2012

scipy.odr реализует ортогональную регрессию расстояния.См. Инструкции по базовому использованию в строке документации:

https://github.com/scipy/scipy/blob/master/scipy/odr/odrpack.py#L27

6 голосов
/ 24 февраля 2012

Я нашел решение.Scipy Odrpack работает норально, но для правильных результатов ему нужно хорошее начальное предположение.Поэтому я разделил процесс на два этапа.

Первый шаг: найдите начальное предположение, используя метод наименьших квадратов ординат.

Второй шаг: замените это первоначальное предположение в ODR как параметр бета0.

И он работает очень хорошо с приемлемой скоростью.

Спасибо, ребята, ваш совет направил меня к правильному решению

0 голосов
/ 21 февраля 2012

Если / когда вы можете инвертировать функцию, описанную p, вы можете просто включить x-pinverted (y) в mFunc, я полагаю, как sqrt (a ^ 2 + b ^ 2), так (псевдокод) 

return sqrt( (y - (p[0]*x**p[1]))^2 + (x - (pinverted(y))^2)

например для 

y=kx+m   p=[m,k]    
pinv=[-m/k,1/k]

return sqrt( (y - (p[0]+x*p[1]))^2 + (x - (pinv[0]+y*pinv[1]))^2)

Но то, что вы просите, в некоторых случаях проблематично. Например, если полиномиальная (или ваша x ^ j) кривая имеет минимальное значение ym при y (m) и у вас есть точка x, y ниже, чем ym, какое значение вы хотите вернуть? Там не всегда решение.

...