Как регрессировать несколько гауссовских пиков из спектрограммы с помощью scipy?

Вот MCVE, показывающий, как регрессировать переменное количество пиков из значительной части вашего сигнала.

import itertools
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import signal, stats, optimize

Мы загружаем вам данные и подвергаем их постобработке, чтобы получить физические ряды:

data = pd.read_csv("20240402_In.mca", encoding = "latin-1", names=["y"]).loc[12:2058]
data["y"] = pd.to_numeric(data["y"])
data["y"] /= data["y"].max()
data["x"] = -0.0188026396003431 + 0.039549044037714 * data.index

data = data.loc[(20. <= data["x"]) & (data["x"] <= 30.), :]

Мы создаем модель, которая принимает переменное количество пиков (определяется решением find_peaks):

def peak(x, A, s, x0):
    law = stats.norm(loc=x0, scale=s)
    return A * law.pdf(x) / law.pdf(x0)

def model(x, *parameters):
    n = len(parameters) // 3
    y = np.zeros_like(x)
    for i in range(n):
        y += peak(x, *parameters[(i * 3):((i + 1) * 3)])
    return y

Выявляем вершины:

indices, bases = signal.find_peaks(data.y, prominence=0.0125)
#(array([ 67, 118, 195, 210]),
# {'prominences': array([0.01987281, 0.99920509, 0.18918919, 0.03338633]),
#  'left_bases': array([  1,   1, 179, 205]),
#  'right_bases': array([ 76, 160, 219, 219])})

Это ключевой момент: на основе идентификации мы создаем обоснованное предположение:

x0s = data.x.values[indices]
As = bases["prominences"]
ss = (data.x.values[bases["right_bases"]] - data.x.values[bases["left_bases"]]) / 8.
p0 = list(itertools.chain(*zip(As, ss, x0s)))
#[0.01987281399046105,
# 0.37077228785356864,
# 22.682348638047493,
# 0.9992050874403816,
# 0.7860372502495658,
# 24.699349883970907,
# 0.1891891891891892,
# 0.19774522018856988,
# 27.744626274874886,
# 0.033386327503974564,
# 0.06921082706599968,
# 28.337861935440593]

Теперь мы регрессируем модель относительно вашего набора данных:

popt, pcov = optimize.curve_fit(model, data.x, data.y, p0=p0)
#array([1.03735804e-02, 9.61270732e-01, 2.29030214e+01, 9.92651381e-01,
#   1.59694755e-01, 2.46561911e+01, 1.85332645e-01, 1.21882422e-01,
#   2.77807838e+01, 3.67805911e-02, 1.21890416e-01, 2.83583849e+01])

Теперь мы можем оценить модель:

yhat = model(data.x, *popt)

Графически это приводит к:

Это вполне удовлетворительно для неограниченной подгонки, но определенно не соответствует ограничениям, которые вы хотите обеспечить.

Не могли бы вы опубликовать в своем посте свои данные, готовые для кодирования, в виде массивов?

На данный момент это самый эффективный метод фиксации параметров в Python, который я обнаружил. В этом случае я использую тот же набор данных, что и раньше. Однако я подгоняю параметры только для первого гауссиана, не учитывая тот факт, что на самом деле это комбинация двух гауссиан.

Я разделил оценку на две части: приблизительную и уточненную. Уточненная аппроксимация ограничивает данные в пределах ±3 сигмы, поскольку она имеет тенденцию давать лучшую производительность.

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes


# Define the Gaussian function
def gaussian(x, amplitude, mean, stddev):
    return amplitude * np.exp(-((x - mean) / stddev)**2 / 2)

# Data that is going to be analysed

x_data= a + b * np.linspace(0, 885, 885)
y_data =data[:885]/np.max(data[:885])


#### Bounds for the first gaussian

# Define bounds for the parameters
bounds = ([0, 24.21, 0], [1, 24.22, 2])
y_data_error =np.sqrt(data[:885])/np.max(data[:885])
# Fit the data to the Gaussian function with bounds
params, covariance = curve_fit(gaussian, x_data, y_data, bounds=bounds)

# Create the model
gaussian_model = Model(gaussian)

# Create parameters with initial values
A_value = 1.0  # Initial value for amplitude
mu_value = 24.21  # Initial value for mean
sigma_value = 1.0  # Initial value for standard deviation
params = gaussian_model.make_params(amplitude=A_value, mean=mu_value, stddev=sigma_value)

# Fix amplitude and mean
#params['amplitude'].vary = False
params['mean'].vary = False

# Fit the model to the data
result = gaussian_model.fit(y_data, params, x=x_data)

# Extract the parameters
amplitude = result.params['amplitude'].value
mean = result.params['mean'].value
stddev = result.params['stddev'].value

# Data interpolation
x_fine = np.linspace(x_data.min(), x_data.max(), num=1000)
y_fit= result.eval(x=x_fine)

ylower_limit = ( (24.21  - 3 * stddev)-a )/b
yupper_limit = ( (24.21  + 3 * stddev)-a )/b

# Selecciona los valores dentro del intervalo

y_limitado = y_data[int(ylower_limit):int(yupper_limit)]

x_limitado = a + b*np.linspace(int(ylower_limit), int(yupper_limit), 24)

ylimit_error=np.sqrt(data[int(ylower_limit):int(yupper_limit)])/np.max(data[int(ylower_limit):int(yupper_limit)])


params1, covariance1 = curve_fit(gaussian, x_limitado, y_limitado, bounds=bounds)

# Create the model
gaussian_model = Model(gaussian)

# Create parameters with initial values
A_value1 = 1.0  # Initial value for amplitude
mu_value1 = 24.21  # Initial value for mean
sigma_value1 = 1.0  # Initial value for standard deviation
params1 = gaussian_model.make_params(amplitude=A_value1, mean=mu_value1, stddev=sigma_value1,sigma=ylimit_error)

# Fix amplitude and mean
#params['amplitude'].vary = False
params1['mean'].vary = False

# Fit the model to the data
result1 = gaussian_model.fit(y_limitado, params1, x=x_limitado)

# Extract the parameters
amplitude1 = result1.params['amplitude'].value
mean1 = result1.params['mean'].value
stddev1 = result1.params['stddev'].value


# Data interpolation
x_fine1 = np.linspace(x_limitado.min(), x_limitado.max(), num=1000)
y_fit1 = result1.eval(x=x_fine1)


fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(x_limitado,y_limitado, color='black', marker='o')
ax.plot(x_fine1,y_fit1,color='blue')
ax.plot(x_fine,y_fit,color='green')

ax.errorbar(x_limitado,y_limitado,yerr=ylimit_error, fmt='o', color='red', label='errorbar')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Gaussian Fit to Data')
ax.set_xlim((24.21  - 3 * stddev, 24.21  + 3 * stddev))
ax.legend()
ax.grid(True)
plt.show()


print("Amplitude:", amplitude)
print("Mean:", mean)
print("Standard Deviation:", stddev)
print("Amplitude:", amplitude1)
print("Mean:", mean1)
print("Standard Deviation:", stddev1)

Редактировать

Теперь у меня есть возможность зафиксировать больше параметров, что позволяет мне точно настроить подгонку суммы первых двух гауссиан. Однако меня озадачивает наличие дополнительного пика в посадке. Я не уверен в его происхождении.

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

# Define Gaussian function
def gaussian(x, amplitude, mean, stddev):
    return amplitude * np.exp(-((x - mean) / stddev)**2 / 2)

# Generate x values for the data
x_data = a + b * np.linspace(0, 885, 885)

# Normalize y data
y_data = data[:885] / np.max(data[:885])
y_data_error = np.sqrt(data[:885]) / np.max(data[:885])

# Fit first Gaussian
bounds = ([1 - 0.0001, 24.21, 0], [1, 24.22, 2])  # Define bounds for the parameters
params, covariance = curve_fit(gaussian, x_data, y_data, bounds=bounds)  # Fit Gaussian to data

# Create the model
gaussian_model = Model(gaussian)

# Initial values for parameters
A_value = 1.0
mu_value = 24.21
sigma_value = 1.0
params = gaussian_model.make_params(amplitude=A_value, mean=mu_value, stddev=sigma_value)
params['mean'].vary = False  # Fix mean parameter

# Fit the model to the data
result = gaussian_model.fit(y_data, params, x=x_data)

# Extract parameters
amplitude = result.params['amplitude'].value
mean = result.params['mean'].value
stddev = result.params['stddev'].value

# Data interpolation
x_fine = np.linspace(x_data.min(), x_data.max(), num=1000)
y_fit = result.eval(x=x_fine)

# Determine y limits for data subset
ylower_limit = ((24.21 - 3 * stddev) - a) / b
yupper_limit = ((24.21 + 3 * stddev) - a) / b

# Select subset of data within the limits
y_limitado = y_data[int(ylower_limit):int(yupper_limit)]
x_limitado = a + b * np.linspace(int(ylower_limit), int(yupper_limit), 24)
ylimit_error = np.sqrt(data[int(ylower_limit):int(yupper_limit)]) / np.max(data[int(ylower_limit):int(yupper_limit)])

# Constrained fit
bounds_sum = ([0, 24.21, 0, 0, 24.0, 0], [1, 24.22, 1, 1, 24.1, 1])  # Bounds for parameters

# Define function for sum of Gaussians
def sum_of_gaussians(x, amplitude1, mean1, stddev1, amplitude2, mean2, stddev2):
    gaussian1 = amplitude1 * np.exp(-((x - mean1) / stddev1)**2 / 2)
    gaussian2 = amplitude2 * np.exp(-((x - mean2) / stddev2)**2 / 2)
    return gaussian1 + gaussian2

# Fit the sum of Gaussians to the subset of data
params1, covariance1 = curve_fit(sum_of_gaussians, x_limitado, y_limitado, bounds=bounds_sum, sigma=ylimit_error)

# Create the model
multigaussian_model = Model(sum_of_gaussians)

# Initial values for parameters
A_value1 = amplitude
mu_value1 = 24.22
sigma_value1 = 0.1

A_value2 = amplitude * 0.5388
mu_value2 = 24.002
sigma_value2 = 0.1

# Set parameters for the model
params1 = multigaussian_model.make_params(amplitude1=A_value1, mean1=mu_value1, stddev1=sigma_value1,
                                          amplitude2=A_value2, mean2=mu_value2, stddev2=sigma_value2, sigma=ylimit_error)

# Fix certain parameters
params1['amplitude1'].vary = False
params1['mean1'].vary = False
params1['amplitude2'].vary = False
params1['mean2'].vary = False

# Fit the model to the data
result1 = multigaussian_model.fit(y_limitado, params1, x=x_limitado)

# Extract parameters
amplitude1 = result1.params['amplitude1'].value
mean1 = result1.params['mean1'].value
stddev1 = result1.params['stddev1'].value

amplitude2 = result1.params['amplitude2'].value
mean2 = result1.params['mean2'].value
stddev2 = result1.params['stddev2'].value

# Data interpolation
x_fine1 = np.linspace(x_limitado.min(), x_limitado.max(), num=1000)
y_fit1 = result1.eval(x=x_fine1)

# Individual Gaussians
x = np.linspace(min(mean1 - 3 * stddev1, mean2 - 3 * stddev2),
                max(mean1 + 3 * stddev1, mean2 + 3 * stddev2), 1000)
y1 = gaussian(x, amplitude1, mean1, stddev1)
y2 = gaussian(x, amplitude2, mean2, stddev2)

# Plotting
fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(x_limitado, y_limitado, color='black', marker='o')
ax.plot(x_fine1, y_fit1, color='blue')
ax.errorbar(x_limitado, y_limitado, yerr=ylimit_error, fmt='o', color='red', label='errorbar')
ax.plot(x, y1, label='Gaussian 1')
ax.plot(x, y2, label='Gaussian 2')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Gaussian Fit to Data')
ax.set_xlim((24.21 - 3 * stddev, 24.21 + 3 * stddev))
ax.legend()
ax.grid(True)
plt.show()

# Print parameters
print("Amplitude1:", amplitude1)
print("Mean1:", mean1)
print("Standard Deviation1:", stddev1)
print("")
print("Amplitude2:", amplitude2)
print("Mean2:", mean2)
print("Standard Deviation2:", stddev2)