使用python scipy.curvefit计算不确定性传播和置信区间

igsr9ssn  于 2022-11-23  发布在  Python
关注(0)|答案(1)|浏览(244)

我正在尝试使用scipy.optimiz.curve_fit来拟合一个函数。我正在用3个高斯函数之和来实现数据的拟合。
这里是必须拟合的数据。

我需要估计每个高斯模型的参数以及这些参数的计算误差,因此我需要计算拟合过程的不确定性,并传播这些数据以显示建模数据的置信区间。
是否可以使用np.sqrt(np.diag(pcov_3gauss))计算出的参数误差进行计算?我在这里看到了:https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html
我应该使用什么方法来传播参数不确定性并计算拟合的不确定性带?
下面是我为拟合实现的代码。

# 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

和拟合结果:

对于第一个脉冲,使用缩放来缩放图像:

E轴重采样比例:

nszi6y05

nszi6y051#

我们可以使用强解析方法计算不确定性带(如果函数已知,就像我的例子一样),但是很难获得所有的导数,即使有可能,也必须进行大量的计算。
我尝试用简化的方法,用蒙特卡罗方法,反复随机抽样估计参数的误差值,计算总体拟合结果。
因此,我计算了所有随机噪声参数的函数结果,并获得了大量结果数据集。
然后,我只使用我得到的计算结果的最大值和最小值--来获得边界的数值。
结果如下所示。x1c 0d1x

相关问题