如何在Python中获取含相关参数的拟合结果的不确定度？

不用自己写大量代码处理相关参数的不确定度，Scipy里有现成工具可以简化流程，下面是具体方案：

1. 用minimize替代fmin获取协方差矩阵 #

fmin是较基础的优化函数，scipy.optimize.minimize更灵活，能直接返回Hessian矩阵的逆（即参数的协方差矩阵），而协方差矩阵正好对应卡方最小值+1的置信区间（1σ水平）。修改后的代码如下：

import numpy as np
from scipy.optimize import minimize

def polynomial(x_value, parameters):
    a, b = parameters
    return (a + b) * x_value**2 + (a - b) * x_value - 2.6

def calculate_chi_squared(parameters, data):
    prediction = polynomial(data[:, 0], parameters)
    return np.sum((prediction - data[:, 1])**2 / data[:, 2]**2)

data = np.genfromtxt('polynomial_data_2.csv', delimiter=',')
INITIAL_GUESS_2 = [0, 0]  # 替换为你的初始猜测值

# 执行优化，选择支持返回Hessian的方法（比如L-BFGS-B）
result = minimize(calculate_chi_squared, INITIAL_GUESS_2, args=(data,), method='L-BFGS-B')
optimised_parameters = result.x
chi2_min = result.fun

2. 从协方差矩阵提取不确定度和相关性 #

通过result.hess_inv拿到协方差矩阵后，直接计算参数的不确定度（对角元平方根）和相关系数：

# 转换为稠密矩阵（L-BFGS-B返回的是稀疏矩阵）
cov_matrix = result.hess_inv.todense()

# 参数a、b的1σ不确定度
a_uncert = np.sqrt(cov_matrix[0, 0])
b_uncert = np.sqrt(cov_matrix[1, 1])

# 参数a和b的相关系数
corr_ab = cov_matrix[0, 1] / (a_uncert * b_uncert)

print(f"最优参数a: {optimised_parameters[0]:.4f} ± {a_uncert:.4f}")
print(f"最优参数b: {optimised_parameters[1]:.4f} ± {b_uncert:.4f}")
print(f"参数相关系数: {corr_ab:.4f}")

3. 快速绘制卡方=卡方_min+1的等高线 #

如果需要可视化置信区间，用网格生成+matplotlib.contour就能实现，不用手动搜索：

import matplotlib.pyplot as plt
from scipy.stats import chi2

# 2个参数的1σ置信区间对应卡方阈值（68.3%置信度）
chi2_threshold = chi2_min + chi2.ppf(0.683, df=2)

# 生成参数网格（以最优参数为中心，覆盖±2倍不确定度范围）
a_grid = np.linspace(optimised_parameters[0]-2*a_uncert, optimised_parameters[0]+2*a_uncert, 100)
b_grid = np.linspace(optimised_parameters[1]-2*b_uncert, optimised_parameters[1]+2*b_uncert, 100)
A, B = np.meshgrid(a_grid, b_grid)

# 批量计算网格点的卡方值
chi2_grid = np.array([calculate_chi_squared([a, b], data) for a, b in zip(A.flatten(), B.flatten())]).reshape(A.shape)

# 绘制等高线和最优参数点
plt.contour(A, B, chi2_grid, levels=[chi2_threshold], colors='red', linewidths=2)
plt.scatter(optimised_parameters[0], optimised_parameters[1], marker='x', color='blue', s=100, label='最优参数')
plt.xlabel('参数a')
plt.ylabel('参数b')
plt.legend()
plt.show()

核心思路 #

卡方最小值+1的等高线对应的是参数的1σ置信区间，而通过优化函数返回的协方差矩阵可以直接得到这个区间的信息，不需要手动遍历参数空间搜索等高线，大幅减少代码量。

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