使用lmfit进行多峰高斯拟合效果不佳的优化求助

我尝试使用lmfit的GaussianModel对光谱数据进行多峰高斯拟合，但部分峰的拟合效果不理想。以下是原始代码：

import matplotlib.pyplot as plt
import numpy as np
from lmfit.models import Model, GaussianModel

fig, axs = plt.subplots(figsize=(7, 5))  
datax= [ 29.5,  30.5,  31.5,  32.5,  33.5,  34.5,  35.5,  36.5,  37.5,
        38.5,  39.5,  40.5,  41.5,  42.5,  43.5,  44.5,  45.5,  46.5,
        47.5,  48.5,  49.5,  50.5,  51.5,  52.5,  53.5,  54.5,  55.5,
        56.5,  57.5,  58.5,  59.5,  60.5,  61.5,  62.5,  63.5,  64.5,
        65.5,  66.5,  67.5,  68.5,  69.5,  70.5,  71.5,  72.5,  73.5,
        74.5,  75.5,  76.5,  77.5,  78.5,  79.5,  80.5,  81.5,  82.5,
        83.5,  84.5,  85.5,  86.5,  87.5,  88.5,  89.5,  91.5,  92.5,
        93.5,  95.5,  96.5, 100.5]

datay  = [     2,  81057, 586959,  15141,  11441,  10208,  16115,  23656,
        55111,  33352,  14368,   9032,   9810,   9344,  15003,  19679,
        16745,   6492,   5387,   4905,   4644,   5202,   7294,   4828,
         2961,   1900,   1920,   1823,   1737,   1521,   1452,    819,
          481,    438,    505,    409,    299,    230,    199,    113,
          109,     78,     91,     84,     46,     39,     33,     27,
           19,     18,     18,     16,      7,      4,      4,      5,
            2,      3,      2,      4,      1,      2,      1,      1,
            1,      2,      1]
 
plt.plot(datax, datay, '-s', markersize=5)
peak_indices   = [ 2,  8, 15, 22]        # indices of detected peaks

peakCentroids = [31.5, 37.5, 44.5, 51.5] # detected peaks

model = None 
for i, cen in enumerate(peakCentroids):
    peak = GaussianModel(prefix='p%d_' %(i))
    if model is None:
        model = peak
    else:
        model = model + peak

pars = model.make_params()
stdv = [0.23, 0.75, 0.68]
for i, cen in enumerate(peakCentroids):   # 0.23, 0.75, 0.68
    pars['p%d_center' % (i)].value = cen
    pars['p%d_amplitude' % (i)].value = datay[peak_indices[i]]            # is that a sensible initial value?
    pars['p%d_sigma' % (i)].value = 0.7+i/10 # sigmas are close to 0.7 

output = model.fit(datay, pars, x=datax)
#print(output.fit_report(min_correl=0.25))
plt.plot(datax,  output.best_fit, 'r-', lw=2)
plt.yscale("log")
plt.xlim(25,  80)
plt.ylim(200, 200000)
plt.show()

优化建议： #

• 修正振幅初始值：lmfit高斯模型的amplitude不是峰高，峰高计算公式为amplitude/(sigma*sqrt(2π))。所以初始振幅应设置为：

  pars['p%d_amplitude' % (i)].value = datay[peak_indices[i]] * pars['p%d_sigma' % (i)].value * np.sqrt(2*np.pi)
  

  这样初始值更贴合模型的物理意义，避免拟合方向偏离。

• 合理设置sigma初始值：你定义了stdv = [0.23, 0.75, 0.68]但未使用，直接将对应峰的已知sigma赋值给初始值更准确，第4个峰可补充合理值（如0.5）：

  stdv = [0.23, 0.75, 0.68, 0.5]
  for i, cen in enumerate(peakCentroids):
      pars['p%d_sigma' % (i)].value = stdv[i]
  

• 添加参数约束：给参数设置合理边界，防止拟合时参数异常漂移：

  pars['p%d_sigma' % (i)].min = 0.1  # sigma不能过小
  pars['p%d_center' % (i)].min = cen - 2
  pars['p%d_center' % (i)].max = cen + 2  # 中心值限制在峰位附近
  pars['p%d_amplitude' % (i)].min = 0  # 振幅非负
  

• 引入基线模型：数据低信号区存在明显基线，当前模型未考虑，导致后续峰拟合受干扰。添加常数基线：

  from lmfit.models import ConstantModel
  model += ConstantModel(prefix='baseline_')
  pars['baseline_c'].value = np.min(datay)  # 用数据最小值作为初始基线
  pars['baseline_c'].min = 0
  

• 使用加权拟合：数据动态范围大（对数刻度），用1/y作为权重让拟合更关注高信号峰：

  weights = 1 / np.array(datay)
  weights[datay == 0] = 1e-6  # 避免除以0
  output = model.fit(datay, pars, x=datax, weights=weights)
  

• 检查峰的数量：x=30.5处（datay=81057）信号强度很高，疑似被遗漏的峰，建议加入模型。也可使用自动峰检测工具：

  from scipy.signal import find_peaks
  peaks, _ = find_peaks(datay, height=10000, distance=3)  # 调整height和distance适配数据
  peakCentroids = [datax[p] for p in peaks]
  peak_indices = peaks.tolist()
  

• 分析拟合报告：取消注释print(output.fit_report(min_correl=0.25))，查看参数误差、相关性。若某参数误差过大，说明初始值或约束不合理，需针对性调整。

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