如何在lmfit中提取伪Voigt峰中心不确定度并绘制误差棒图

• 使用Python的lmfit包进行数据拟合：对多组数据采用PseudoVoigt+常数模型拟合，提取峰中心位置（约100个），再用正弦波模型对这些位置拟合并绘图。
• 遇到的问题：无法获取PseudoVoigt拟合得到的峰中心位置的不确定度，尝试eval_uncertainty函数仅能获取峰的振幅和sigma值的不确定度，无法提取峰中心的置信区间。
• 需求：将第二张图的散点图替换为带有PseudoVoigt峰中心置信区间的误差棒图。

lmfit的ModelResult拟合结果对象中，参数的不确定度可通过以下方式获取：

• 标准误差：直接访问out.params['参数名'].stderr
• 置信区间：调用out.conf_interval()方法计算（默认95%置信水平）

修改步骤：

1. 更新Voigt_app函数，返回峰中心的标准误差和置信区间
2. 收集所有峰中心的位置、标准误差/置信区间
3. 使用plt.errorbar()替代plt.scatter()绘制带误差棒的图

import numpy as np
import matplotlib.pyplot as plt
from lmfit.models import PseudoVoigtModel, ConstantModel, SineModel

def Voigt_app(x,y,a,b):
    psvo_mod = PseudoVoigtModel(prefix='psvo_')
    pars = psvo_mod.guess(y,x=x)
    
    con_mod = ConstantModel(prefix='con_')
    pars = pars.update(con_mod.make_params(c = (y[0]+y[1])/2))
    
    mod = psvo_mod + con_mod
    out = mod.fit(y, pars, x=x)
    
    # 获取峰中心的标准误差
    center_stderr = out.params['psvo_center'].stderr
    # 计算95%置信区间
    ci = out.conf_interval()
    center_ci_low, center_ci_high = ci['psvo_center'][1][1], ci['psvo_center'][2][1]
    
    details = out.fit_report(min_correl=0.25)
    return out, details, center_stderr, (center_ci_low, center_ci_high)

def Sin_app(x,y):
    sin_mod = SineModel(prefix='sin_')
    pars = sin_mod.guess(y,x=x)
    
    con_mod = ConstantModel(prefix='con_')
    pars = pars.update(con_mod.make_params(c = np.mean(y)))
    
    mod = sin_mod + con_mod
    out = mod.fit(y,pars,x=x)
    details = out.fit_report(min_correl=0.25)
                       
    return out, details

data = np.asarray(
      [[14407.0 , 26989.0 , 31026.0 , 31194.0 , 31302.0 , 91996.0 , 112250.0 , 112204.0 , 113431.0 , 75097.0 , 55942.0 , 55826.0 , 56181.0 , 28266.0 , 14687.0],
       [12842.0 , 13275.0 , 13581.0 , 13943.0 , 14914.0 , 20463.0 , 21132.0 , 21457.0 , 23308.0 , 32017.0 , 30808.0 , 30927.0 , 30337.0 , 21025.0 , 15216.0],
       [15770.0 , 17677.0 , 19008.0 , 20911.0 , 25958.0 , 27984.0 , 28164.0 , 31161.0 , 33517.0 , 29430.0 , 28673.0 , 32725.0 , 28495.0 , 24527.0 , 25173.0],
       [16299.0 , 20067.0 , 25102.0 , 34968.0 , 37757.0 , 44871.0 , 51347.0 , 60154.0 , 54985.0 , 53383.0 , 45776.0 , 40836.0 , 30816.0 , 27922.0 , 26066.0]])
a,b = 1932, 1947
centers = []
center_stderrs = []
center_cis = []  # 存储(下限, 上限)

colors = ['red','black','green','blue']
fig, ax = plt.subplots()
x = np.arange(15)  # 补充原代码缺失的x定义
for i in range(4):
    ax.scatter(x, data[i,:], label=f'data {i}', color=colors[i])
    out, details, stderr, ci = Voigt_app(x, data[i,:], a, b)
    ax.plot(x, out.best_fit, label=f'fit {i}', color=colors[i])
    centers.append(out.values['psvo_center'])
    center_stderrs.append(stderr)
    center_cis.append(ci)
ax.set_title('Data points and pseudo-Voigt fits')
ax.legend()


fig, ax = plt.subplots()
# 选项1：用标准误差绘制误差棒
ax.errorbar(np.arange(4), centers, yerr=center_stderrs, fmt='ro', capsize=5, label='data with stderr')
# 选项2：用95%置信区间绘制误差棒（取消注释即可）
# ax.errorbar(np.arange(4), centers, yerr=np.array(center_cis).T, fmt='ro', capsize=5, label='data with 95% CI')

out, details = Sin_app(np.arange(4), np.asarray(centers))
ax.plot(np.arange(4), out.best_fit, color='black', label='fit')
ax.set_title('Peak positions and sine-wave fit')
ax.legend()

plt.show()

• 原代码遗漏了x的定义，已补充x = np.arange(15)。
• out.params['psvo_center'].stderr直接返回峰中心参数的标准误差，是拟合后参数的标准差估计值。
• out.conf_interval()返回参数的置信区间字典，其中ci['psvo_center']的第2和第3个元素分别对应95%置信区间的下限和上限。
• errorbar的yerr参数支持两种输入：一维数组表示对称误差（标准误差），二维数组表示上下不对称的置信区间（需转置后传入）。

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