多高斯曲线拟合结果异常、标准差近乎一致问题求助

问题描述 #

编写代码拟合MCA数据，使用自定义read_file函数提取数据，但高斯拟合结果不符合预期，所有拟合的高斯函数标准差几乎一致（约1e-09），拟合曲线与实际数据偏差较大。

原始代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import os

file_path='/content/drive/MyDrive/ProjetoXRF_Manuel/April2024/Amostra4.mca' # You need to change this to your directory

def read_file(file_path):
    # Abrir el archivo
    with open(file_path, 'r', encoding='latin1') as f:
        data = []
        live_time = None  # Inicializar variable live_time
        real_time = None  # Inicializar variable real_time
        for line in f:
            if 'LIVE_TIME' in line:
                live_time = line.strip()  # Almacenar la línea que contiene LIVE_TIME
            elif 'REAL_TIME' in line:
                real_time = line.strip()  # Almacenar la línea que contiene REAL_TIME
            elif all(char.isdigit() or char in ".-+e" for char in line.strip()):
                # Analizar y almacenar datos numéricos como flotante único
                row = [float(x) for x in line.split()]
                data.extend(row)

    # Convertir 'data' a un arreglo de numpy
    return np.array(data), live_time, real_time

data, live_time, real_time = read_numeric_data(file_path)
a = -0.0076702251954372525
b = 0.03952691936704189
Data=data

peaks, _ = find_peaks(Data[:, 0] / maxim, height=0.1)

# Get peak values and indices
peak_values = Data[peaks] / maxim
peak_indices = peaks

# Definición de la función gaussiana
def gaus(x, a, x0, sigma):
    return a * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))
limite=(10-a)/b
print(limite)
# Datos
x_values = a + b * np.linspace(0, 270, 270)
y_values = np.squeeze(Data[:270]) / np.max(Data[:270])

# Parámetros iniciales
sigma = 2 # initial value
standar=[]
# Ajuste de curvas para diferentes conjuntos de datos
for i in range(len(peak_values)):
    mean = a + b * peak_indices[i] 
    amplitude = peak_values[i][0]
    
    # Ajustar la curva para cada conjunto de datos
    popt, pcov = curve_fit(gaus, x_values, y_values, p0=[amplitude, mean, 1])
    
    # Visualización del ajuste
    m,n,p= popt
    plt.plot(x_values, gaus(x_values,amplitude,mean,p), label='fit')

    print(mean, amplitude,"Desviación standar",p)

    standar.append(p)

plt.scatter(x_values, y_values, label='Data')
plt.xlim(4, 10)
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.grid(True)
plt.show()

问题排查与修复 #

1. 函数调用错误 #

定义了read_file函数，但调用时写成了read_numeric_data，直接导致数据无法读取，需修正为：

data, live_time, real_time = read_file(file_path)

2. 缺失依赖与未定义变量 #

• 未导入find_peaks，需从scipy.signal导入：

  from scipy.signal import find_peaks
  

• maxim变量未定义，需先计算数据最大值：

  maxim = np.max(Data[:270])
  

3. 数据维度错误 #

Data是一维数组，Data[:, 0]索引方式无效，需改为直接使用Data：

peaks, _ = find_peaks(Data[:270] / maxim, height=0.1)
peak_values = Data[peaks] / maxim

4. 全局拟合导致优化失败 #

每次拟合都使用全部数据，噪声干扰大，优化器无法收敛到正确的sigma值。应限制拟合范围为峰值附近的局部区域：

# 对每个峰值，只拟合附近±5个通道的区域
window = 5
for i in range(len(peak_values)):
    peak_idx = peak_indices[i]
    # 确定拟合的x、y范围
    fit_x = x_values[peak_idx - window : peak_idx + window]
    fit_y = y_values[peak_idx - window : peak_idx + window]
    
    mean = a + b * peak_idx 
    amplitude = peak_values[i]
    
    # 仅对局部区域拟合
    popt, pcov = curve_fit(gaus, fit_x, fit_y, p0=[amplitude, mean, 0.5])
    m, n, p = popt
    
    # 绘制拟合曲线时，用拟合参数而非初始值
    plt.plot(x_values, gaus(x_values, m, n, p), label=f'fit {i+1}')

5. 绘图参数错误 #

原代码绘图时使用初始参数（amplitude, mean）而非拟合得到的popt参数，导致拟合曲线无法反映优化结果，需修正为使用m, n, p（即popt的三个值）。

修复后完整代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.signal import find_peaks
import os

file_path='/content/drive/MyDrive/ProjetoXRF_Manuel/April2024/Amostra4.mca' # 修改为你的文件路径

def read_file(file_path):
    with open(file_path, 'r', encoding='latin1') as f:
        data = []
        live_time = None
        real_time = None
        for line in f:
            if 'LIVE_TIME' in line:
                live_time = line.strip()
            elif 'REAL_TIME' in line:
                real_time = line.strip()
            elif all(char.isdigit() or char in ".-+e" for char in line.strip()):
                row = [float(x) for x in line.split()]
                data.extend(row)
    return np.array(data), live_time, real_time

# 读取数据
data, live_time, real_time = read_file(file_path)
a = -0.0076702251954372525
b = 0.03952691936704189
Data = data

# 计算最大值，用于归一化
maxim = np.max(Data[:270])
# 寻找峰值
peaks, _ = find_peaks(Data[:270] / maxim, height=0.1)
peak_values = Data[peaks] / maxim
peak_indices = peaks

# 高斯函数定义
def gaus(x, amp, x0, sigma):
    return amp * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))

# 生成x轴值
x_values = a + b * np.linspace(0, 270, 270)
y_values = np.squeeze(Data[:270]) / maxim

standar = []
# 局部拟合每个峰值
window = 5  # 峰值附近的拟合窗口大小
plt.scatter(x_values, y_values, label='Data', s=5)

for i in range(len(peak_values)):
    peak_idx = peak_indices[i]
    # 确保窗口不超出数据范围
    start_idx = max(0, peak_idx - window)
    end_idx = min(len(x_values), peak_idx + window)
    fit_x = x_values[start_idx:end_idx]
    fit_y = y_values[start_idx:end_idx]
    
    mean = a + b * peak_idx 
    amplitude = peak_values[i]
    
    # 拟合局部数据
    popt, pcov = curve_fit(gaus, fit_x, fit_y, p0=[amplitude, mean, 0.5])
    amp_fit, mean_fit, sigma_fit = popt
    
    # 绘制拟合曲线
    plt.plot(x_values, gaus(x_values, amp_fit, mean_fit, sigma_fit), label=f'Peak {i+1} fit')
    
    print(f"峰值位置: {mean:.2f}, 峰值幅度: {amplitude:.2f}, 标准差: {sigma_fit:.4f}")
    standar.append(sigma_fit)

plt.xlim(4, 10)
plt.xlabel('Energy')
plt.ylabel('Normalized Intensity')
plt.legend()
plt.grid(True)
plt.show()

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