python - 线性回归 ODR 失败

Question

遵循 this answer 中的建议我为 beta0 使用了多种值组合，如此处所示，这些值来自 polyfit。

此示例已更新，以显示 X 值与 Y 值的相对比例的影响(X 范围是 Y 的 0.1 到 100 倍):

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 range(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 拟合相似的唯一方法是交换 X 和 Y，或者如此处所示增加 X 相对于 Y 的值范围，但仍然没有真正的解决方案:)

=== 编辑 ===

我不希望 ODR 与 polyfit 相同。我展示 polyfit 只是为了强调 ODR 拟合是错误的，这不是数据的问题。

=== 解决方案 ===

感谢 @norok2 在 Y 范围为 X 的 0.001 到 100000 倍时的回答:

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 range(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()

Answer 1

polyfit() 和正交距离回归 (ODR) 拟合之间的主要区别是 polyfit 在 x 上的误差可以忽略不计的假设下工作。如果违反此假设，就像在您的数据中一样，您不能指望这两种方法会产生相似的结果。特别是，ODR() 对您指定的错误非常敏感。如果您没有指定任何错误/权重，它将为 x 和 y 分配一个值 1，这意味着两者之间的任何比例差异x 和 y 会影响结果(即所谓的数值调节)。

相反，polyfit() 在计算拟合之前，对数据应用某种预白化(参见其 source code 的第 577 行左右)以获得更好的数值条件。

因此，如果您希望 ODR() 匹配 polyfit()，您可以简单地微调 Y 上的错误以更改你的数值条件。我测试过这适用于 1e-10 和 Y 的 1e10 之间的任何数值条件(它是 /10。 或 1e-1 在你的例子中)。

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-10 和 Y 的 1e10 之间的任何数值条件(它是 /10 . 或 1e-1 在你的例子中)。

请注意，这仅对条件良好的合身有意义。