多高斯拟合振幅固定比例设置及标准差负值问题咨询

我编写了一段Python代码用于拟合MCA数据，通过read_numeric_data函数读取数据。目前有两个需求需要解决：

1. 多高斯拟合中为各峰振幅设置固定比例关系，具体规则为：

   peak_amplitudes = [1, 0.5996, 0.8841, 0.8840.51837]
   即第二个高斯的振幅为第一个的0.5996倍，第三个与第四个高斯的振幅按对应比例关联。

2. 代码中已固定高斯的均值，但拟合过程中出现部分标准差为负值的问题，需要修正。

当前拟合结果图：

原始代码 #

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

# Function to read numeric data from a file
def read_numeric_data(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

file_path = '/content/drive/MyDrive/ProjetoXRF_Manuel/April2024/20240402_In.mca'
data, _, _ = read_numeric_data(file_path)

a = -0.0188026396003431
b = 0.039549044037714
Data = data
dim = 0

# Function to convolve multiple Gaussians
def multi_gaussian(x, *params):
    eps = 1e-10
    y = np.zeros_like(x)
    for i in range(0, len(params), 3):
        amplitude, mean, sigma = params[i:i+3]
        # Initial parameters for Gaussian peaks
        peak_means = [24.21, 24.002, 27.276, 27.238]  # mean known value
        mean=peak_means[i//3] # to fix the mean value whith the known values
        y += amplitude * np.exp(-(x - mean)**2 / (2 * sigma**2 + eps))
    return y

sigma = []
area = []

# Function to plot the convolved energy spectrum
def plot_convolved_spectrum(Data, a, b, i, ax=None):
    maxim = np.max(Data)
    workingarray = np.squeeze(Data)

    # Define peak points
    peaks, _ = find_peaks(workingarray / maxim, height=0.1)
    peak_values = workingarray[peaks] / maxim
    peak_indices = peaks

    # Calculate energy values corresponding to the peaks
    energy_spectrum = a + b * peak_indices

    # Prepare data for convolution
    data = workingarray[:885] / maxim
    data_y = data / data.max()
    data_x = a + b * np.linspace(0, 885, num=len(data_y))

    # Adjust initial guesses for peak amplitudes based on desired relationships
    peak_amplitudes = [1, 0.5996, 0.884*1, 0.884*0.51837]
    peak_amplitudes[1] = 0.5996 * peak_amplitudes[0]
    peak_amplitudes[2] = 0.884 * peak_amplitudes[0]
    peak_amplitudes[3] = 0.884 * 0.51837 * peak_amplitudes[0]

    peak_sigmas = [0.1] * 4  # to fix the initial values for sigma


    peak_means = [24.21, 24.002, 27.276, 27.238]  # mean known value

    params_init = list(zip(peak_amplitudes, peak_means, peak_sigmas))
    params_init = np.concatenate(params_init)

    # Fit multiple Gaussians to the energy spectrum
    params, params_cov = curve_fit(multi_gaussian, data_x, data_y, p0=params_init)

    # Get a fine interpolation of the fit
    x_fine = np.linspace(data_x.min(), data_x.max(), num=20000)

    # Plot each individual Gaussian with a different color
    for j in range(0, len(params), 3):
        amplitude, mean, sigma = params[j:j+3]
        ax.plot(x_fine, amplitude * np.exp(-(x_fine - mean)**2 / (2 * sigma**2)),
                label=f"Gaussian {j//3 + 1}", linewidth=1.5)

    # Other plotting configurations
    ax.set_xlabel('Energy (KeV)')
    ax.set_ylabel('Normalized Data')
    ax.legend()
    ax.set_title('Convolved Energy Spectrum')

    # Print information
    sigmas_array = params[2::3]
    amplitudes_array = params[::3]
    means_array = params[1::3]
    areas = [peak_amplitudes[i] * np.sqrt(2 * np.pi) * sigmas_array[i] for i in range(len(sigmas_array))]
    area.append(areas)
    total_area = np.sum(areas)
    print("amplitudes",amplitudes_array)
    print("Standard deviations:", sigmas_array)
    print("Areas:", areas)

# Plotting
fig, ax = plt.subplots()
plot_convolved_spectrum(Data, a, b, dim, ax=ax)
ax.set_xlim([a + b * (peak_indices[0] - 30), a + b * (peak_indices[-1] + 50)])
plt.show()

一、实现振幅固定比例 #

原始代码仅设置了初始猜测的比例，但curve_fit仍会将所有振幅作为独立参数优化。需要修改高斯函数，让所有振幅依赖于一个主参数，严格按比例计算：

修改后的multi_gaussian函数 #

def multi_gaussian(x, *params):
    eps = 1e-10
    y = np.zeros_like(x)
    # 固定振幅比例
    amp_ratios = [1, 0.5996, 0.884, 0.884 * 0.51837]
    peak_means = [24.21, 24.002, 27.276, 27.238]  # 固定均值
    # 参数结构：[主振幅, sigma1, sigma2, sigma3, sigma4]
    main_amp = params[0]
    sigmas = params[1:]
    
    for i in range(4):
        amplitude = main_amp * amp_ratios[i]
        mean = peak_means[i]
        sigma = sigmas[i]
        y += amplitude * np.exp(-(x - mean)**2 / (2 * sigma**2 + eps))
    return y

调整初始参数与拟合逻辑 #

# 初始猜测：主振幅设为1，四个sigma初始为0.1
params_init = [1] + [0.1] * 4
# 拟合时仅优化主振幅和四个sigma
params, params_cov = curve_fit(multi_gaussian, data_x, data_y, p0=params_init)

二、解决标准差负值问题 #

标准差本质为正数，出现负值是因为curve_fit未设置参数约束，推荐两种解决方法：

方法1：设置参数边界（推荐） #

通过bounds参数强制sigma大于0：

# 参数下界：主振幅≥0，sigma≥极小值（避免除以0）
lower_bounds = [0] + [1e-6] * 4
# 参数上界：无限制
upper_bounds = [np.inf] + [np.inf] * 4
params, params_cov = curve_fit(multi_gaussian, data_x, data_y, p0=params_init, bounds=(lower_bounds, upper_bounds))

方法2：取sigma绝对值 #

在高斯函数中对sigma取绝对值，即使拟合出负值也能保证函数合理性：

y += amplitude * np.exp(-(x - mean)**2 / (2 * (abs(sigma))**2 + eps))

修改后的完整代码 #

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

def read_numeric_data(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

file_path = '/content/drive/MyDrive/ProjetoXRF_Manuel/April2024/20240402_In.mca'
data, _, _ = read_numeric_data(file_path)

a = -0.0188026396003431
b = 0.039549044037714
Data = data
dim = 0

def multi_gaussian(x, *params):
    eps = 1e-10
    y = np.zeros_like(x)
    amp_ratios = [1, 0.5996, 0.884, 0.884 * 0.51837]
    peak_means = [24.21, 24.002, 27.276, 27.238]
    main_amp = params[0]
    sigmas = params[1:]
    
    for i in range(4):
        amplitude = main_amp * amp_ratios[i]
        mean = peak_means[i]
        sigma = sigmas[i]
        y += amplitude * np.exp(-(x - mean)**2 / (2 * sigma**2 + eps))
    return y

sigma = []
area = []

def plot_convolved_spectrum(Data, a, b, i, ax=None):
    maxim = np.max(Data)
    workingarray = np.squeeze(Data)

    peaks, _ = find_peaks(workingarray / maxim, height=0.1)
    peak_values = workingarray[peaks] / maxim
    peak_indices = peaks

    data = workingarray[:885] / maxim
    data_y = data / data.max()
    data_x = a + b * np.linspace(0, 885, num=len(data_y))

    # 初始参数调整为[主振幅, sigma1, sigma2, sigma3, sigma4]
    params_init = [1] + [0.1]*4
    lower_bounds = [0] + [1e-6]*4
    upper_bounds = [np.inf] + [np.inf]*4
    params, params_cov = curve_fit(multi_gaussian, data_x, data_y, p0=params_init, bounds=(lower_bounds, upper_bounds))

    x_fine = np.linspace(data_x.min(), data_x.max(), num=20000)
    amp_ratios = [1, 0.5996, 0.884, 0.884 * 0.51837]
    peak_means = [24.21, 24.002, 27.276, 27.238]
    
    for j in range(4):
        amplitude = params[0] * amp_ratios[j]
        mean = peak_means[j]
        sigma = params[j+1]
        ax.plot(x_fine, amplitude * np.exp(-(x_fine - mean)**2 / (2 * sigma**2)),
                label=f"Gaussian {j+1}", linewidth=1.5)

    ax.set_xlabel('Energy (KeV)')
    ax.set_ylabel('Normalized Data')
    ax.legend()
    ax.set_title('Convolved Energy Spectrum')

    sigmas_array = params[1:]
    amplitudes_array = [params[0] * ratio for ratio in amp_ratios]
    areas = [amp * np.sqrt(2 * np.pi) * sigma for amp, sigma in zip(amplitudes_array, sigmas_array)]
    area.append(areas)
    
    print("amplitudes", amplitudes_array)
    print("Standard deviations:", sigmas_array)
    print("Areas:", areas)

fig, ax = plt.subplots()
plot_convolved_spectrum(Data, a, b, dim, ax=ax)
ax.set_xlim([a + b * (peak_indices[0] - 30), a + b * (peak_indices[-1] + 50)])
plt.show()

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