基于Scipy的双高斯拟合效果自动评估方法求助

第一步：修正拟合函数（关键错误） #

你当前的拟合函数存在逻辑错误，amp1 = amp1 - m*x - C会将标量参数转为数组，导致模型偏离“双高斯+线性背景”的目标。正确的函数应该让高斯峰的振幅是相对于背景的峰值高度，代码如下：

import numpy as np
from scipy.optimize import curve_fit

def gaussian2_same_wid(x, amp1, width1, amp2, m, C):
    # 线性背景
    background = m * x + C
    # 两个固定峰位的高斯峰（同宽度）
    f1 = amp1 * np.exp(-((x - 0.671644)**2) / (2 * width1**2))
    f2 = amp2 * np.exp(-((x - 0.673081)**2) / (2 * width1**2))
    # 总模型：峰+背景
    return f1 + f2 + background

对应的初始猜测和边界也需要调整，因为amp1和amp2现在是相对于背景的峰值，所以初始值可以用区间最大值减去背景估计，边界也应该设为正数（峰高不可能为负）：

# 估计背景（取x两端的均值）
bg_est = np.mean(y[(x < 0.6709) | (x > 0.675)])
max_SII_1 = np.max(y[(x >= 0.6709) & (x <= 0.67275)]) - bg_est
max_SII_2 = np.max(y[(x >= 0.67275) & (x <= 0.675)]) - bg_est

guess_same = (max_SII_1, 0.00045, max_SII_2, 0, bg_est)
bounds_same = [
    [0, 1e-5, 0, -5, 0],  # amp1/amp2下限为0，峰宽最小1e-5，背景参数范围
    [max_SII_1 + 0.2, 0.001, max_SII_2 + 0.2, 5, 3]  # 峰宽上限0.001（根据领域知识调整）
]

# 拟合
popt_same, pcov_same = curve_fit(
    gaussian2_same_wid, xdata=x, ydata=y, sigma=dy,
    p0=guess_same, bounds=bounds_same, maxfev=1000000
)

第二步：可靠的拟合评估指标 #

1. 缩减卡方（Reduced Chi-Squared） #

这是连续数据拟合的标准统计量，反映模型对数据的解释能力：

from scipy.stats import chi2

# 计算残差
residuals = y - gaussian2_same_wid(x, *popt_same)
# 计算卡方值
chi_sq = np.sum((residuals / dy)**2)
# 自由度：数据点数量 - 拟合参数数量（这里是5个）
dof = len(x) - len(popt_same)
red_chi_sq = chi_sq / dof

# 可选：计算p值（原假设：模型能解释数据）
p_value = 1 - chi2.cdf(chi_sq, dof)

• 判断标准：red_chi_sq ≈ 1说明拟合良好；red_chi_sq > 3通常表示拟合差或噪声估计不准；p_value < 0.05说明模型无法解释数据（需舍弃）。

2. 参数合理性校验 #

拟合参数必须符合物理意义，否则直接舍弃：

• 峰高（amp1, amp2）：必须为正（高斯峰是信号，高于背景），若为0或负数，说明峰被噪声淹没或拟合失败。
• 峰宽（width1）：根据SII峰的领域知识，设置合理范围（比如1e-5 < width1 < 0.001），超出范围说明拟合跑偏。
• 参数不确定性：计算参数的标准差（协方差矩阵对角线开根号），若某参数的标准差远大于参数本身（比如std(amp1) > amp1 * 0.5），说明该参数无法被可靠拟合，样本质量差。

# 计算参数标准差
param_std = np.sqrt(np.diag(pcov_same))
# 检查峰高是否为正且不确定性合理
valid_amp = (popt_same[0] > 0) and (popt_same[2] > 0) and (param_std[0] < popt_same[0]*0.5) and (param_std[2] < popt_same[2]*0.5)
# 检查峰宽是否在合理范围
valid_width = (popt_same[1] > 1e-5) and (popt_same[1] < 0.001)

3. 信噪比（SNR） #

直接评估目标峰的信号强度相对于噪声的比例，是筛选低信噪比样本的核心指标：

# 计算两个峰位的拟合信号值
peak1_y = gaussian2_same_wid(0.671644, *popt_same)
peak2_y = gaussian2_same_wid(0.673081, *popt_same)
# 计算对应峰位的背景值
peak1_bg = popt_same[3] * 0.671644 + popt_same[4]
peak2_bg = popt_same[3] * 0.673081 + popt_same[4]
# 峰的净信号高度
peak1_signal = peak1_y - peak1_bg
peak2_signal = peak2_y - peak2_bg
# 噪声水平：残差的标准差（或用输入的dy的均值）
noise_level = np.std(residuals)
# 计算SNR
snr1 = peak1_signal / noise_level
snr2 = peak2_signal / noise_level

# 阈值：SNR >= 3（统计上3σ为显著信号）
valid_snr = (snr1 >= 3) and (snr2 >= 3)

4. 残差分布检查 #

拟合良好的样本，残差应满足：

• 均值接近0（无系统偏差）
• 无明显趋势（比如残差随x增大而单调变化，说明模型漏了成分）

# 检查残差均值是否接近0（容忍度±0.1倍噪声水平）
valid_residual_mean = np.abs(np.mean(residuals)) < noise_level * 0.1
# 可选：检查残差是否随机（用线性拟合残差与x的R²，若R²<0.1则无明显趋势）
from scipy.stats import linregress
slope, _, r_value, _, _ = linregress(x, residuals)
valid_residual_trend = r_value**2 < 0.1

第三步：自动化筛选流程 #

将上述指标组合，编写批量处理函数：

def is_valid_fit(x, y, dy):
    # 1. 拟合模型（捕获curve_fit的异常）
    try:
        # 计算初始猜测
        bg_est = np.mean(y[(x < 0.6709) | (x > 0.675)])
        max_SII_1 = np.max(y[(x >= 0.6709) & (x <= 0.67275)]) - bg_est
        max_SII_2 = np.max(y[(x >= 0.67275) & (x <= 0.675)]) - bg_est
        guess_same = (max_SII_1, 0.00045, max_SII_2, 0, bg_est)
        bounds_same = [[0, 1e-5, 0, -5, 0], [max_SII_1+0.2, 0.001, max_SII_2+0.2, 5, 3]]
        
        popt, pcov = curve_fit(gaussian2_same_wid, x, y, sigma=dy, p0=guess_same, bounds=bounds_same, maxfev=1e6)
    except RuntimeError:
        # 拟合不收敛，直接标记为无效
        return False
    
    # 2. 计算评估指标
    residuals = y - gaussian2_same_wid(x, *popt)
    chi_sq = np.sum((residuals/dy)**2)
    dof = len(x) - len(popt)
    red_chi_sq = chi_sq / dof
    param_std = np.sqrt(np.diag(pcov))
    
    # 3. 校验所有条件
    valid_amp = (popt[0] > 0) and (popt[2] > 0) and (param_std[0] < popt[0]*0.5) and (param_std[2] < popt[2]*0.5)
    valid_width = (popt[1] > 1e-5) and (popt[1] < 0.001)
    valid_red_chi = red_chi_sq < 3
    
    # 计算SNR
    peak1_signal = gaussian2_same_wid(0.671644, *popt) - (popt[3]*0.671644 + popt[4])
    peak2_signal = gaussian2_same_wid(0.673081, *popt) - (popt[3]*0.673081 + popt[4])
    noise_level = np.std(residuals)
    valid_snr = (peak1_signal/noise_level >=3) and (peak2_signal/noise_level >=3)
    
    # 残差检查
    valid_residual = np.abs(np.mean(residuals)) < noise_level*0.1
    
    # 所有条件都满足则返回True
    return all([valid_amp, valid_width, valid_red_chi, valid_snr, valid_residual])

# 批量处理示例
valid_indices = []
for i in range(len(all_data_arrays)):
    x = all_x[i]
    y = all_y[i]
    dy = all_dy[i]
    if is_valid_fit(x, y, dy):
        valid_indices.append(i)

关键说明 #

• 阈值（如SNR≥3、red_chi_sq<3）可根据你的样本分布微调，建议先在已知的优质/劣质样本上测试，调整到合适的区分度。
• 若curve_fit经常不收敛，可优化初始猜测（比如用滑动窗口找峰的近似高度），或放宽maxfev值。

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