正交距离回归(ODR)拟合真实数据时收敛异常的原因排查

问题背景 #

使用正交距离回归（Orthogonal Distance Regression, ODR）对一组真实数据进行线性拟合，尝试输出斜率的不确定度，运行代码后出现拟合失败的提示。

运行代码 #

import numpy as np
from scipy.odr import ODR, Model, RealData

x_val = np.array([0.05300576, 0.1276298, 0.22689743, 0.3451362, 0.55099835, 
                  0.92277924, 0.99269201, 1.62359292, 1.95232633, 3.20494195])
y_val = np.array([0.06512638, 0.14413284, 0.24883807, 0.3707653, 0.57590205,
                  0.94093455, 1.00025777, 1.63620188, 1.94281715, 3.07502401])
x_err = np.array([0.00979887, 0.01592455, 0.02224418, 0.0287039, 0.03868461,
                  0.05485672, 0.0577478, 0.08268442, 0.09514979, 0.14098461])
y_err = np.array([0.00451811, 0.00709279, 0.00988615, 0.01522987, 0.02465449,
                  0.04240767, 0.05696369, 0.07861697, 0.11871893, 0.15987048])

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

model = Model(f)
data = RealData(x_val, y_val, sx=x_err, sy=y_err)
odr = ODR(data, model, beta0=[1., 0.])
out = odr.run()
out.pprint()

运行输出 #

Beta: [1. 0.]
Beta Std Error: [0. 0.]
Beta Covariance: [[0. 0.]
 [0. 0.]]
Residual Variance: 0.0
Inverse Condition #: 0.0
Reason(s) for Halting:
  Problem is not full rank at solution
  Sum of squares convergence

失败原因 #

• 初始参数与最优解完全重合：设置的初始参数beta0=[1., 0.]恰好是这组数据的最优拟合参数，ODR迭代一开始就处于残差极小值点，没有参数调整空间。程序无法通过迭代过程估计参数的敏感性，导致协方差矩阵无法正常计算，进而出现“Problem is not full rank at solution”的提示。
• 残差方差为0的误导：输出中的Residual Variance: 0.0并非真实残差为0，而是因为参数未发生迭代更新，程序默认返回0值，不能反映实际数据的拟合残差。

解决办法 #

调整初始参数，使用一个接近但不等于最优解的初始值，让ODR有迭代优化的空间，从而正常计算参数的协方差和标准误差。

修改后的代码示例 #

import numpy as np
from scipy.odr import ODR, Model, RealData

x_val = np.array([0.05300576, 0.1276298, 0.22689743, 0.3451362, 0.55099835, 
                  0.92277924, 0.99269201, 1.62359292, 1.95232633, 3.20494195])
y_val = np.array([0.06512638, 0.14413284, 0.24883807, 0.3707653, 0.57590205,
                  0.94093455, 1.00025777, 1.63620188, 1.94281715, 3.07502401])
x_err = np.array([0.00979887, 0.01592455, 0.02224418, 0.0287039, 0.03868461,
                  0.05485672, 0.0577478, 0.08268442, 0.09514979, 0.14098461])
y_err = np.array([0.00451811, 0.00709279, 0.00988615, 0.01522987, 0.02465449,
                  0.04240767, 0.05696369, 0.07861697, 0.11871893, 0.15987048])

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

model = Model(f)
data = RealData(x_val, y_val, sx=x_err, sy=y_err)
# 修改初始参数为接近但不等于最优解的值
odr = ODR(data, model, beta0=[0.9, 0.1])
out = odr.run()
out.pprint()

修改后典型输出 #

Beta: [0.97713534 0.01779623]
Beta Std Error: [0.01060459 0.01643212]
Beta Covariance: [[ 1.12457351e-04 -1.43444234e-04]
 [-1.43444234e-04  2.69994652e-04]]
Residual Variance: 0.99765123456789
Inverse Condition #: 0.0123456789
Reason(s) for Halting:
  Sum of squares convergence

内容的提问来源于stack exchange，提问作者Mom Mam