python - 线性回归ODR失败

按照this answer中的建议，我为beta0使用了几个值的组合，如图所示，polyfit中的值。
更新此示例是为了显示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()

似乎beta0的组合对。。。获得polyfit和ODR fit相似的唯一方法是交换X和Y，或者如图所示，增加X相对于Y的值范围，但仍然不是真正的解决方案：）

==编辑===
我不希望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()

最佳答案

polyfit()与正交距离回归（ODR）拟合的关键区别在于，polyfit的工作假设是x上的误差可以忽略不计。如果违反了这个假设，就像在数据中一样，就不能期望这两个方法产生类似的结果。
尤其是，ODR()对您指定的错误非常敏感。
如果您没有指定任何误差/权重，它将为1和x指定值y，这意味着x和y之间的任何刻度差都将影响结果（所谓的数值调节）。
相反，polyfit()在计算拟合之前，会对数据应用某种预白化（参见其source code的577行周围）以获得更好的数值条件。
因此，如果希望ODR()与polyfit()匹配，只需微调Y上的错误即可更改数值条件。
我测试了这对于1e-10和1e10之间的任何数值条件作用（在您的示例中是Y或/ 10.）。

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

致：

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

（编辑：注意上面一行有错别字）
我测试了这对于1e-1和1e-10之间的任何数值条件作用（在您的示例中是1e10或Y）。
注意，这只适用于条件良好的配合。

关于python - 线性回归ODR失败，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/52723212/