基于Python的线性拟合误差传播及未测点y值求解

嘿，这个问题我熟！你手里的每个y值都带了测量不确定度dy，所以得用加权线性最小二乘拟合（给不确定度小的数据点更大权重，让拟合结果更靠谱），之后再基于拟合出的线性模型预测未测量点，同时通过误差传播计算预测值的不确定度。下面分两种方法来实现：

方法一：手动实现（理解底层原理） #

先手动推导计算，帮你搞懂每个步骤的逻辑：

import numpy as np

# 你的原始数据
x = np.array([1, 2, 3, 4])
y = np.array([4.1, 5.8, 8.1, 9.7])
dy = np.array([0.2, 0.3, 0.2, 0.4])

# 计算权重：不确定度越小，权重越大，这里取 1/(dy²)
weights = 1 / dy**2

# 计算加权求和项，用于求解拟合参数
sum_w = np.sum(weights)
sum_wx = np.sum(weights * x)
sum_wy = np.sum(weights * y)
sum_wx2 = np.sum(weights * x**2)
sum_wxy = np.sum(weights * x * y)

# 求解斜率m和截距b
denominator = sum_w * sum_wx2 - sum_wx**2
m = (sum_w * sum_wxy - sum_wx * sum_wy) / denominator
b = (sum_wx2 * sum_wy - sum_wx * sum_wxy) / denominator

print(f"拟合得到线性模型：y = {m:.4f}x + {b:.4f}")

# 计算参数m和b的协方差矩阵，用于误差传播
cov_matrix = np.array([[sum_wx2, sum_wx], [sum_wx, sum_w]]) / denominator
var_m = cov_matrix[0, 0]  # m的方差
var_b = cov_matrix[1, 1]  # b的方差
cov_mb = cov_matrix[0, 1] # m和b的协方差

# 预测未测量点x_unm对应的y值（这里举个例子，你可以替换成自己需要的点）
x_unm = 2.5
y_unm = m * x_unm + b

# 误差传播计算预测值的不确定度
# 公式：var(y_unm) = x_unm²*var(m) + var(b) + 2*x_unm*cov(m,b)
var_y_unm = x_unm**2 * var_m + var_b + 2 * x_unm * cov_mb
dy_unm = np.sqrt(var_y_unm)

print(f"当x={x_unm}时，预测y值为 {y_unm:.4f}，不确定度为 ±{dy_unm:.4f}")

方法二：用Scipy简化实现（高效便捷） #

从Scipy 1.7.0版本开始，linregress支持传入权重参数，能直接给出拟合参数及其不确定度，代码更简洁：

from scipy.stats import linregress
import numpy as np

# 原始数据
x = np.array([1, 2, 3, 4])
y = np.array([4.1, 5.8, 8.1, 9.7])
dy = np.array([0.2, 0.3, 0.2, 0.4])

# 权重设置和手动方法一致
weights = 1 / dy**2

# 执行加权线性拟合
fit_result = linregress(x, y, weights=weights)

# 提取拟合参数和不确定度
m = fit_result.slope
b = fit_result.intercept
m_err = fit_result.stderr       # 斜率m的标准差
b_err = fit_result.intercept_stderr  # 截距b的标准差
cov_mb = fit_result.rvalue * m_err * b_err  # m和b的协方差（通过相关系数推导）

print(f"拟合得到线性模型：y = {m:.4f}x + {b:.4f}")
print(f"斜率m的不确定度：±{m_err:.4f}，截距b的不确定度：±{b_err:.4f}")

# 预测未测量点
x_unm = 2.5
y_unm = m * x_unm + b

# 计算预测值的不确定度
var_y_unm = x_unm**2 * m_err**2 + b_err**2 + 2 * x_unm * cov_mb
dy_unm = np.sqrt(var_y_unm)

print(f"当x={x_unm}时，预测y值为 {y_unm:.4f}，不确定度为 ±{dy_unm:.4f}")

关键注意点 #

• 加权的意义：如果忽略dy直接做普通线性拟合，会给不确定度大的数据点不合理的权重，导致拟合结果偏差。
• 误差传播不能忽略协方差：预测值的不确定度不仅和m、b各自的方差有关，还和两者的协方差有关，漏掉协方差项会导致不确定度计算不准确。
• 如果你需要预测多个未测量点，可以把x_unm改成数组，用向量运算批量计算，效率更高。

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