使用 scipy.ODR 的线性回归失败（解决方案中的排名不全）答案

【问题标题】：Linear Regression using scipy.ODR fails (Not full rank at solution)使用 scipy.ODR 的线性回归失败（解决方案中的排名不全）
【发布时间】：2015-02-01 07:47:38
【问题描述】：

因此尝试使用 scipy.odr 进行线性回归。然而，它惨遭失败。 scipy.odr 以前为我工作过，我没有在我的代码中看到任何错误。我能想到的唯一原因是斜率可能太小，但我不明白这会如何打扰 scipy。感谢您的帮助。

代码：

#!/usr/bin/env python3 -i
# -*- coding: iso-8859-1 -*-

import matplotlib.pyplot as plt
import numpy as np
from scipy.odr import *


fig = plt.figure()
ax1 = fig.add_subplot(111)


x   = np.linspace(0,1E15,10)
y   = 1E-15*x-2



ax1.set_xlim(-0.05E15,1.1E15)
ax1.set_ylim(-2.1, -0.7)


ax1.plot(x, y, 'o')

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

linear = Model(f)
mydata = RealData(x, y)
myodr = ODR(mydata, linear, beta0=[1., 2.])
myoutput = myodr.run()
myoutput.pprint()

a, b = myoutput.beta
sa, sb = myoutput.sd_beta

xp = np.linspace(ax1.get_xlim()[0], ax1.get_xlim()[1], 1000)
yp = a*xp+b


ax1.plot(xp,yp)

plt.show()

这是最终的终端输出：

Beta: [ -4.84615388e-15   2.00000000e+00]
Beta Std Error: [  8.14077323e-16   0.00000000e+00]
Beta Covariance: [[  1.46153845e-31   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00]]
Residual Variance: 4.534412955465587
Inverse Condition #: 1.0
Reason(s) for Halting:
  Problem is not full rank at solution
  Parameter convergence

这是生成的图形：

编辑：我的 odr 回归代码来自 http://docs.scipy.org/doc/scipy/reference/odr.html

【问题讨论】：

你从myoutput.Beta得到的系数不正确，如果你愿意可以使用numpy.polyfit
如果想要简单的线性回归，可以使用numpy.polyfit
感谢您的回复基鲁巴哈兰。我知道他们错了（从图片中你可以清楚地看到）。（我之前在这组数据上尝试过 polyfit，它可以工作）。但是：我不想使用 polyfit，因为它不会返回任何系数错误。而且我也不想使用 scipy.optimize.curve_fit （这会）。我喜欢 scipy.ODR，因为它可以解释我的 x 和 y 数据中的错误。（这个数据集中不需要，那只是我的问题的简化版）
哦，我还尝试让使用正交距离回归（因此得名）的 scipy.ODR 使用最小二乘算法。它没有改变任何东西..（或者我没有告诉 scipy 这样做）

标签： python matplotlib scipy regression one-definition-rule

【解决方案1】：

import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure()
ax1 = fig.add_subplot(111)
x   = np.linspace(0,1E15,10)
y   = 1E-15*x-2    
ax1.set_xlim(-0.05E15,1.1E15)
ax1.set_ylim(-2.1, -0.7)    
ax1.plot(x, y, 'o')
# Fit using odr


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

sx = np.std(x)
sy = np.std(y)
linear = Model(f)
mydata = RealData(x=x,y=y, sx=sx, sy=sy)
myodr = ODR(mydata, linear, beta0=[1.00000000e-15, 2.])
myoutput = myodr.run()
myoutput.pprint()

a, b = myoutput.beta
sa, sb = myoutput.sd_beta

xp = np.linspace(min(x), max(x), 1000)
yp = a*xp+b
ax1.plot(xp,yp)
plt.show()

【讨论】：

最重要的是你在beta0中给出的初始值。首先使用 np.polyfit 估计系数，并在 beta0 中使用。