使用scipy curve_fit求解加权最小二乘直线的斜率及拟合误差

初始实现与问题 #

实现加权线性拟合的基础代码如下：

from scipy.optimize import curve_fit

# 定义线性拟合函数，a为斜率，b为截距
def func(x, a, b):
    return a * x + b

# 测试数据
x = [0, 1, 2, 3]
y = [160, 317, 3302, 16002]
# 每个y值的测量误差，用于加权最小二乘
yerr = [0.0791, 0.0562, 0.0174, 0.0079]

最初尝试提取斜率的写法存在问题，且不清楚如何获取拟合误差：

# 错误写法：仅接收一个返回值
popt = curve_fit(func, x, y, sigma=yerr)
slope = popt[0]

正确实现方案 #

curve_fit调用后会返回两个结果：

• 一维数组popt：存储拟合得到的最优参数，和自定义函数的参数顺序一一对应，此处popt[0]为斜率，popt[1]为截距
• 二维数组pcov：参数的协方差矩阵，对角线元素对应各参数的方差，参数的拟合误差为对应方差的平方根

完整正确代码示例 #

import numpy as np
from scipy.optimize import curve_fit

def func(x, a, b):
    return a * x + b

x = [0, 1, 2, 3]
y = [160, 317, 3302, 16002]
yerr = [0.0791, 0.0562, 0.0174, 0.0079]

# 同时接收两个返回值
popt, pcov = curve_fit(func, x, y, sigma=yerr)
# 提取斜率、截距
slope, intercept = popt
# 提取斜率误差、截距误差
slope_err = np.sqrt(pcov[0,0])
intercept_err = np.sqrt(pcov[1,1])

原理说明 #

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