Линейная регрессия ODR не работает - PullRequest
Новая
       36

Линейная регрессия ODR не работает

0 голосов
/ 09 октября 2018

Следуя рекомендациям в этот ответ Я использовал несколько комбинаций значений для бета0, и, как показано здесь, значения из полифита.

Этот пример ОБНОВЛЕН, чтобы показатьЭффект относительных шкал значений X по сравнению с Y: 

from random import random, seed
from scipy import polyfit
from scipy import odr
import numpy as np
from matplotlib import pyplot as plt

seed(1)
X = np.array([random() for i in range(1000)])
Y = np.array([i + random()**2 for i in range(1000)])

for num in xrange(1, 5):
    plt.subplot(2, 2, num)
    plt.title('X range is %.1f times Y' % (float(100 / max(X))))
    X *= 10
    z = np.polyfit(X, Y, 1)
    plt.plot(X, Y, 'k.', alpha=0.1)

    # Fit using odr
    def f(B, X):
        return B[0]*X + B[1]    

    linear = odr.Model(f)
    mydata = odr.RealData(X, Y)
    myodr = odr.ODR(mydata, linear, beta0=z)
    myodr.set_job(fit_type=0)
    myoutput = myodr.run()
    a, b = myoutput.beta
    sa, sb = myoutput.sd_beta
    xp = np.linspace(plt.xlim()[0], plt.xlim()[1], 1000)
    yp = a*xp+b
    plt.plot(xp, yp, label='ODR')
    yp2 = z[0]*xp+z[1]
    plt.plot(xp, yp2, label='polyfit')
    plt.legend()
    plt.ylim(-1000, 2000)
plt.show()

Кажется, что никакая комбинация бета0 не помогает ... Единственный способ получить подобие polyfit и ODR - поменять местами X и Y, OR какпоказанный здесь, чтобы увеличить диапазон значений X относительно Y, все еще не действительно решение:)

new example

=== РЕДАКТИРОВАТЬ ===

Я не хочу, чтобы ODR был таким же, как polyfit.Я показываю polyfit только для того, чтобы подчеркнуть, что подгонка ODR неверна и это не проблема данных.

=== РЕШЕНИЕ ===

благодаря ответу @ norok2: 

from random import random, seed
from scipy import polyfit
from scipy import odr
import numpy as np
from matplotlib import pyplot as plt
seed(1)
X = np.array([random() / 1000 for i in range(1000)])
Y = np.array([i + random()**2 for i in range(1000)])
plt.figure(figsize=(12, 12))
for num in xrange(1, 10):
    plt.subplot(3, 3, num)
    plt.title('Y range is %.1f times X' % (float(100 / max(X))))
    X *= 10
    z = np.polyfit(X, Y, 1)
    plt.plot(X, Y, 'k.', alpha=0.1)
    # Fit using odr
    def f(B, X):
        return B[0]*X + B[1]    
    linear = odr.Model(f)
    mydata = odr.RealData(X, Y, 
                          sy=min(1/np.var(Y), 1/np.var(X)))  # here the trick!! :)
    myodr = odr.ODR(mydata, linear, beta0=z)
    myodr.set_job(fit_type=0)
    myoutput = myodr.run()
    a, b = myoutput.beta
    sa, sb = myoutput.sd_beta
    xp = np.linspace(plt.xlim()[0], plt.xlim()[1], 1000)
    yp = a*xp+b
    plt.plot(xp, yp, label='ODR')
    yp2 = z[0]*xp+z[1]
    plt.plot(xp, yp2, label='polyfit')

    plt.legend()
    plt.ylim(-1000, 2000)
plt.show()

example3

Ответы [ 2 ]

0 голосов
/ 09 октября 2018

Ключевая разница между polyfit() и подбором ортогональной регрессии (ODR) заключается в том, что полифит работает в предположении, что ошибка на x незначительна.Если это предположение нарушается, как это происходит в ваших данных, вы не можете ожидать, что оба метода приведут к схожим результатам.В частности, ODR() очень чувствителен к указанным вами ошибкам.Если вы не укажете какую-либо ошибку / взвешивание, ему будет присвоено значение 1 для x и y, означающее, что любая разница в шкале между x и y повлияет на результаты (поэтомуназываемое числовым условием).

Напротив, polyfit() перед вычислением подбора применяет какое-то предварительное отбеливание к данным (см. строку 577 его исходного кода )для лучшей обработки чисел.

Поэтому, если вы хотите, чтобы ODR() соответствовал polyfit(), вы можете просто настроить ошибку на Y, чтобы изменить вашу обработку чисел.Я проверил, что это работает для любого числового условия между 1e-10 и 1e10 вашего Y (в вашем примере это / 10. или 1e-1). 

mydata = odr.RealData(X, Y)
# equivalent to: odr.RealData(X, Y, sx=1, sy=1)

to: 

mydata = odr.RealData(X, Y, sx=1, sy=1/np.var(Y))

(РЕДАКТИРОВАТЬ: обратите внимание, что в строке выше была опечатка)

Я проверил, что это работает для любого числового условия между 1e-10 и 1e10 вашего Y (этов вашем примере это / 10. или 1e-1).

Обратите внимание, что это имеет смысл только для хорошо подготовленных припадков.

0 голосов
/ 09 октября 2018

Я не могу отформатировать исходный код в комментарии, поэтому разместите его здесь.Этот код использует ODR для расчета статистики соответствия, обратите внимание на строку с «порядком параметров для odr», так что я использую функцию-обертку для вызова ODR моей «фактической» функции. 

from scipy.optimize import curve_fit
import numpy as np
import scipy.odr
import scipy.stats

x = np.array([5.357, 5.797, 5.936, 6.161, 6.697, 6.731, 6.775, 8.442, 9.861])
y = np.array([0.376, 0.874, 1.049, 1.327, 2.054, 2.077, 2.138, 4.744, 7.104])

def f(x,b0,b1):
    return b0 + (b1 * x)


def f_wrapper_for_odr(beta, x): # parameter order for odr
    return f(x, *beta)

parameters, cov= curve_fit(f, x, y)

model = scipy.odr.odrpack.Model(f_wrapper_for_odr)
data = scipy.odr.odrpack.Data(x,y)
myodr = scipy.odr.odrpack.ODR(data, model, beta0=parameters,  maxit=0)
myodr.set_job(fit_type=2)
parameterStatistics = myodr.run()
df_e = len(x) - len(parameters) # degrees of freedom, error
cov_beta = parameterStatistics.cov_beta # parameter covariance matrix from ODR
sd_beta = parameterStatistics.sd_beta * parameterStatistics.sd_beta
ci = []
t_df = scipy.stats.t.ppf(0.975, df_e)
ci = []
for i in range(len(parameters)):
    ci.append([parameters[i] - t_df * parameterStatistics.sd_beta[i], parameters[i] + t_df * parameterStatistics.sd_beta[i]])

tstat_beta = parameters / parameterStatistics.sd_beta # coeff t-statistics
pstat_beta = (1.0 - scipy.stats.t.cdf(np.abs(tstat_beta), df_e)) * 2.0    # coef. p-values

for i in range(len(parameters)):
    print('parameter:', parameters[i])
    print('   conf interval:', ci[i][0], ci[i][1])
    print('   tstat:', tstat_beta[i])
    print('   pstat:', pstat_beta[i])
    print()
...