基于scipy.curvefit的不确定度传播与置信区间计算问询

问题 #

使用scipy.optimize.curve_fit实现3个高斯函数之和的数据拟合，已完成拟合代码并得到参数结果：

拟合代码 #

# Initial data for fitting
x_array = np.array(sep_df.E)
y_array_3gauss = np.array(sep_df.exp_cs)

def _1gaussian(x, amp1,cen1,sigma1,offset):
    return amp1*(1/(sigma1*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x-cen1)/sigma1)**2)))+offset

def _3gaussian(x, amp1,cen1,sigma1,amp2,cen2,sigma2,amp3,cen3,sigma3,offset):
    return _1gaussian(x, amp1,cen1,sigma1,offset=0) + \
        _1gaussian(x, amp2,cen2,sigma2,offset=0) + \
           _1gaussian(x, amp3,cen3,sigma3,offset=0) + offset

#initial_guesses for Gaussians
amp1 = 100 #max value without an offset (!)
cen1 = 140 # position of a center
sigma1 = 1 # sd of a gaussian, can be calculated approx. as  HWHM / 2.355 

amp2 = 32
cen2 = 157
sigma2 = 1

amp3 = 17.5
cen3 = 171.5
sigma3 = 1

offset_initial_guess = y_array_3gauss.mean()

p0=[amp1, cen1, sigma1, 
amp2, cen2, sigma2, 
amp3, cen3, sigma3, 
offset_initial_guess]

# using a scipy.optimize.curve_fit for parameters Estimation
popt_3gauss, pcov_3gauss = scipy.optimize.curve_fit(_3gaussian, x_array, y_array_3gauss, p0=p0)
perr_3gauss = np.sqrt(np.diag(pcov_3gauss)) # errors (??)

print('Popt_3 gauss')
print(popt_3gauss)

i=0
for param in popt_3gauss:
    print(f'Guess: {p0[i]} -> value: {param} (+/-) {perr_3gauss[i]}')
    i+=1

pars_1 = np.append(popt_3gauss[0:3], popt_3gauss[9])
pars_2 = np.append(popt_3gauss[3:6], popt_3gauss[9])
pars_3 = np.append(popt_3gauss[6:9], popt_3gauss[9])

#calculating the separate Gaussians
gauss_peak_1 = _1gaussian(x_array, *pars_1)
gauss_peak_2 = _1gaussian(x_array, *pars_2)
gauss_peak_3 = _1gaussian(x_array, *pars_3)

参数输出结果 #

Guess: 100 -> value: 19.921886501569567 (+/-) 0.18211089486661997
Guess: 140 -> value: 140.8226385680359 (+/-) 0.0009978640529532633
Guess: 1 -> value: 0.07977753969265024 (+/-) 0.0008591843799752477
Guess: 32 -> value: 5.8061836613068865 (+/-) 0.21223980806115864
Guess: 157 -> value: 157.24985139555835 (+/-) 0.005092072398387486
Guess: 1 -> value: 0.08218041022663795 (+/-) 0.0034588647851462877
Guess: 17.5 -> value: 4.183133300983996 (+/-) 0.2522036049333162
Guess: 171.5 -> value: 171.47025791173272 (+/-) 0.008818904590601183
Guess: 1 -> value: 0.11713718144344663 (+/-) 0.008042004990244404
Guess: 4.743919339218138 -> value: 4.016878986311514 (+/-) 0.04473028381895628

需要解决的问题：

• 能否通过np.sqrt(np.diag(pcov_3gauss))计算参数误差？
• 应采用何种方法传播参数不确定度并计算拟合的不确定带？

回答 #

1. 参数误差计算：np.sqrt(np.diag(pcov_3gauss))是有效的 #

scipy.optimize.curve_fit返回的pcov是参数的协方差矩阵，其对角线元素对应每个参数的方差。对这些对角线元素取平方根，得到的就是参数的标准误差（即参数估计值的标准差），这个计算完全正确。

注意前提条件：

• 结果基于拟合的线性近似假设，即模型在最优参数popt附近可近似为线性函数；
• curve_fit默认假设数据残差服从正态分布且方差恒定（同方差性）。若数据不符合这些假设，标准误差的可靠性会下降，但仍是行业通用的参数不确定度估计方法。

2. 不确定度传播与拟合置信区间计算 #

方法一：Delta方法（线性不确定度传播）

这是最常用的方法，基于参数协方差矩阵，通过计算模型对参数的偏导数，将参数的不确定度传播到拟合值上，步骤如下：

1. 对每个x值，计算模型_3gaussian(x, *p)对所有参数的偏导数（即雅可比矩阵的行向量）；
2. 用该向量与协方差矩阵pcov做矩阵乘法，再与向量的转置相乘，得到该x处拟合值的方差；
3. 对方差取平方根得到拟合值的标准误差，再乘以置信水平对应的分位数（比如95%置信区间用1.96），得到置信区间上下限。

示例代码：

import numpy as np
from scipy.optimize import curve_fit

# 定义计算雅可比矩阵的函数（手动计算偏导数）
def jac_3gaussian(x, amp1,cen1,sigma1,amp2,cen2,sigma2,amp3,cen3,sigma3,offset):
    # 对amp1的偏导
    d_amp1 = (1/(sigma1*np.sqrt(2*np.pi))) * np.exp(-0.5*((x-cen1)/sigma1)**2)
    # 对cen1的偏导
    d_cen1 = amp1 * (x - cen1) / (sigma1**3 * np.sqrt(2*np.pi)) * np.exp(-0.5*((x-cen1)/sigma1)**2)
    # 对sigma1的偏导
    term1 = (x-cen1)**2 / sigma1**4
    term2 = 1 / sigma1**2
    d_sigma1 = amp1 / np.sqrt(2*np.pi) * (term1 - term2) * np.exp(-0.5*((x-cen1)/sigma1)**2)
    # 对amp2、cen2、sigma2的偏导
    d_amp2 = (1/(sigma2*np.sqrt(2*np.pi))) * np.exp(-0.5*((x-cen2)/sigma2)**2)
    d_cen2 = amp2 * (x - cen2) / (sigma2**3 * np.sqrt(2*np.pi)) * np.exp(-0.5*((x-cen2)/sigma2)**2)
    term1_2 = (x-cen2)**2 / sigma2**4
    term2_2 = 1 / sigma2**2
    d_sigma2 = amp2 / np.sqrt(2*np.pi) * (term1_2 - term2_2) * np.exp(-0.5*((x-cen2)/sigma2)**2)
    # 对amp3、cen3、sigma3的偏导
    d_amp3 = (1/(sigma3*np.sqrt(2*np.pi))) * np.exp(-0.5*((x-cen3)/sigma3)**2)
    d_cen3 = amp3 * (x - cen3) / (sigma3**3 * np.sqrt(2*np.pi)) * np.exp(-0.5*((x-cen3)/sigma3)**2)
    term1_3 = (x-cen3)**2 / sigma3**4
    term2_3 = 1 / sigma3**2
    d_sigma3 = amp3 / np.sqrt(2*np.pi) * (term1_3 - term2_3) * np.exp(-0.5*((x-cen3)/sigma3)**2)
    # 对offset的偏导
    d_offset = np.ones_like(x)
    # 组合成雅可比矩阵（每行对应一个x的偏导数）
    jac = np.column_stack([d_amp1, d_cen1, d_sigma1, d_amp2, d_cen2, d_sigma2, d_amp3, d_cen3, d_sigma3, d_offset])
    return jac

# 计算拟合值的标准误差
jac = jac_3gaussian(x_array, *popt_3gauss)
y_var = np.einsum('ij,jk,ik->i', jac, pcov_3gauss, jac)
y_std = np.sqrt(y_var)

# 计算95%置信区间
confidence_level = 0.95
z_score = 1.96  # 正态分布分位数
y_upper = _3gaussian(x_array, *popt_3gauss) + z_score * y_std
y_lower = _3gaussian(x_array, *popt_3gauss) - z_score * y_std

方法二：蒙特卡洛采样（非线性不确定度传播）

如果模型非线性较强，Delta方法的线性近似误差较大，蒙特卡洛采样更准确：

1. 从参数的多元正态分布中采样（均值为popt，协方差矩阵为pcov）；
2. 对每一组采样参数，计算对应的拟合曲线；
3. 对每个x值，取所有采样拟合值的分位数（比如2.5%和97.5%分位数）作为置信区间上下限。

示例代码：

import numpy as np
from scipy.stats import multivariate_normal

# 采样次数（次数越多结果越稳定）
n_samples = 10000

# 从参数的多元正态分布中采样
params_samples = multivariate_normal.rvs(mean=popt_3gauss, cov=pcov_3gauss, size=n_samples)

# 计算所有采样的拟合曲线
y_samples = np.array([_3gaussian(x_array, *params) for params in params_samples])

# 计算95%置信区间（取2.5%和97.5%分位数）
y_lower = np.percentile(y_samples, 2.5, axis=0)
y_upper = np.percentile(y_samples, 97.5, axis=0)

3. 单个高斯峰的不确定度 #

若需单独计算每个高斯峰的不确定带，只需将上述方法中的_3gaussian替换为对应的单个高斯函数（_1gaussian），并传入对应参数的协方差子矩阵即可。例如，第一个高斯峰的参数是popt_3gauss[0:3] + [popt_3gauss[9]]，对应的协方差子矩阵是pcov_3gauss[np.ix_([0,1,2,9], [0,1,2,9])]。

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