如何在Python3中对携带非对称误差棒的数据进行拟合

核心原理 #

对于带有非对称误差的拟合，不能直接使用默认对称误差的最小二乘方法，需要根据残差的正负动态选择对应的上下误差作为权重计算依据，本质是基于半高斯分布假设的加权最小二乘，也可以根据需求选择极大似然估计、贝叶斯拟合等方案。

方案1：基于scipy的自定义实现（无需额外安装第三方库） #

实现步骤 #

1. 首先将原始的[y_best, y_min, y_max]格式转换为上下误差值：
   • 下误差：y_err_low = y_best - y_min
   • 上误差：y_err_high = y_max - y_best
2. 自定义非对称残差计算逻辑，残差为正时用上误差做权重，残差为负时用下误差做权重
3. 调用scipy.optimize.least_squares执行拟合

完整代码示例 #

import numpy as np
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

# ---------------------- 1. 准备数据 ----------------------
# 自变量x数组
x = np.array([1, 2, 3, 4, 5, 6])
# 每个点对应格式[y_best, y_min, y_max]
raw_data = np.array([
    [2.1, 1.8, 2.5],
    [3.5, 3.2, 4.0],
    [4.8, 4.2, 5.3],
    [6.3, 5.8, 6.9],
    [7.7, 7.1, 8.3],
    [9.2, 8.6, 9.8]
])
y_best = raw_data[:, 0]
y_err_low = y_best - raw_data[:, 1]
y_err_high = raw_data[:, 2] - y_best

# ---------------------- 2. 定义拟合模型 ----------------------
# 示例为一阶线性模型，可替换为任意你需要的拟合函数（如二次函数、指数函数等）
def fit_model(params, x):
    slope, intercept = params
    return slope * x + intercept

# ---------------------- 3. 定义非对称残差函数 ----------------------
def asymmetric_residual(params, x, y_best, y_err_low, y_err_high):
    y_pred = fit_model(params, x)
    residual = y_pred - y_best
    # 动态选择误差作为权重
    sigma = np.where(residual > 0, y_err_high, y_err_low)
    return residual / sigma

# ---------------------- 4. 执行拟合 ----------------------
# 初始参数猜测，根据你的模型调整
init_params = [1, 0]
fit_result = least_squares(asymmetric_residual, init_params, args=(x, y_best, y_err_low, y_err_high))
slope_fit, intercept_fit = fit_result.x
print(f"拟合结果：斜率={slope_fit:.3f}, 截距={intercept_fit:.3f}")

# ---------------------- 5. 结果可视化（可选） ----------------------
plt.errorbar(x, y_best, yerr=[y_err_low, y_err_high], fmt='o', capsize=5, label='原始数据')
plt.plot(x, fit_model(fit_result.x, x), 'r-', label=f'拟合曲线: y={slope_fit:.2f}x+{intercept_fit:.2f}')
plt.legend()
plt.show()

方案2：基于lmfit库的简化实现 #

如果需要更完善的拟合参数误差评估、拟合优度计算等功能，可以使用专门的拟合库lmfit，首先执行pip install lmfit安装依赖，实现逻辑更简洁，也内置了非对称误差的支持。

import numpy as np
from lmfit import Model
import matplotlib.pyplot as plt

# 数据准备同上
x = np.array([1, 2, 3, 4, 5, 6])
raw_data = np.array([
    [2.1, 1.8, 2.5],
    [3.5, 3.2, 4.0],
    [4.8, 4.2, 5.3],
    [6.3, 5.8, 6.9],
    [7.7, 7.1, 8.3],
    [9.2, 8.6, 9.8]
])
y_best = raw_data[:, 0]
y_err_low = y_best - raw_data[:, 1]
y_err_high = raw_data[:, 2] - y_best

# 定义拟合模型
def linear_model(x, slope, intercept):
    return slope * x + intercept

model = Model(linear_model)
params = model.make_params(slope=1, intercept=0)

# 执行非对称误差拟合
result = model.fit(y_best, params, x=x, weights=lambda res, *_: np.where(res>0, 1/y_err_high, 1/y_err_low))
print(result.fit_report())

# 可视化
result.plot()
plt.show()

注意事项 #

• 上述方案默认假设误差服从半高斯分布，符合大部分科研场景的需求，如果你的误差属于泊松、二项等特殊分布，建议改用对应分布的极大似然估计实现拟合
• 拟合前请先校验数据合法性，避免出现y_min > y_best或y_max < y_best的异常值
• 可以通过计算校正卡方值（残差平方和除以自由度，即数据点数减拟合参数个数）评估拟合优度，理想值接近1

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