如何用Scipy ODR拟合幂律并获取点与拟合曲线的垂直残差？

Scipy ODR的delta属性返回的是模型预测值与观测值的纵向残差（|y_model - y_data|），确实没有内置方法直接输出数据点到拟合幂律曲线的垂直（最短）距离。要手动计算，步骤如下：

1. 明确幂律模型形式 #

拟合的幂律模型通常为：

y = a * x^b

其中a和b是ODR拟合得到的最优参数。

2. 推导最小距离的求解方程 #

对于单个数据点(x_i, y_i)，曲线上距离它最近的点(x*, y*)满足：两点连线与曲线在(x*, y*)处的切线垂直。结合幂律模型的导数（切线斜率），可得到关于x*的方程：

(x_i - x*) + a*b*x*^(b-1)*(y_i - a*x*^b) = 0

这个方程没有解析解，需要用数值方法求解x*。

3. 数值求解并计算距离 #

使用Scipy的数值优化工具（如scipy.optimize.root或scipy.optimize.minimize_scalar）求解x*，再计算对应的y* = a*x*^b，最后用欧氏距离公式计算垂直距离。

代码示例 #

import numpy as np
from scipy.odr import ODR, Model, Data
from scipy.optimize import root

# 定义幂律模型
def power_law(p, x):
    a, b = p
    return a * (x ** b)

# 替换为你的实际数据
x_data = np.linspace(1, 10, 20)
y_data = 2 * (x_data ** 1.5) + np.random.normal(0, 2, size=len(x_data))

# ODR拟合流程
data = Data(x_data, y_data)
model = Model(power_law)
odr = ODR(data, model, beta0=[1, 1])  # 初始参数猜测
output = odr.run()

# 提取拟合参数
a_fit, b_fit = output.beta

# 定义求解x*的方程
def equation(x_star, x_i, y_i, a, b):
    return (x_i - x_star) + a*b*(x_star ** (b-1))*(y_i - a*(x_star ** b))

# 逐个计算每个数据点的垂直距离
vertical_distances = []
for x_i, y_i in zip(x_data, y_data):
    # 用数据点x值作为初始猜测，提升求解稳定性
    guess = x_i
    sol = root(equation, guess, args=(x_i, y_i, a_fit, b_fit))
    x_star = sol.x[0]
    y_star = a_fit * (x_star ** b_fit)
    dist = np.sqrt((x_i - x_star)**2 + (y_i - y_star)**2)
    vertical_distances.append(dist)

vertical_distances = np.array(vertical_distances)
print("垂直距离数组：", vertical_distances)

注意事项 #

• 如果你的幂律模型带有常数项或其他变形，需要同步调整模型函数和求解方程。
• 若数值求解出现不收敛情况，可改用scipy.optimize.minimize_scalar直接最小化平方距离函数，有时稳定性更好。
• 初始猜测值建议贴合数据范围，比如用模型在x_i处的对应趋势值，避免求解陷入局部最优。

内容的提问来源于stack exchange，提问作者Chryssi Koukouraki