线性回归 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 范围为 0.001 到 100000 倍 X 时的回答：

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 上的错误以更改您的数值条件。 我测试了这适用于Y 的1e-10 和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))

（编辑：请注意上面一行有错字）

我测试了这适用于Y 的1e-10 和1e10 之间的任何数值条件（在您的示例中为/ 10. 或1e-1）。

请注意，这只对条件良好的拟合有意义。

Answer 2

我无法在注释中格式化源代码，所以把它放在这里。此代码使用 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()