Как вы подгоняете данные измеренных линий эмиссии с помощью функции Гаусса в Python? (Атомная спектроскопия)

Question

В физическом лабораторном проекте я измеряю различные эмиссионные линии от различных элементов. Пики высокой интенсивности возникают на определенных длинах волн. Моя цель - установить гауссовскую функцию в python, чтобы определить, на какой длине волны интенсивность достигает пика.

Я уже пытался использовать функцию нормы из библиотеки scipy.stats. Ниже приведен код и полученный график.

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

mean, std = norm.fit(he3888_1[:,0])
plt.plot(he3888_1[:,0], he3888_1[:,1], color='r')
x = np.linspace(min(he3888_1[:,0]), max(he3888_1[:,0]), 100)
y = norm.pdf(x, mean, std)
plt.plot(x, y)
plt.xlabel("Wavelength (Angstroms)")
plt.ylabel("Intensity")
plt.show()

Может ли это быть из-за низкой интенсивности в течение относительно длительного периода до нее?

Answer 1

Lmfit кажется хорошим вариантом для вашего случая. Приведенный ниже код имитирует пик Гаусса с добавленным линейным фоном и показывает, как вы можете извлечь параметры с помощью lmfit. Последний имеет ряд других встроенных моделей (Lorentzian, Voight и т. Д.), Которые можно легко комбинировать друг с другом.

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

def generate_gaussian(amp, mu, sigma_sq, slope=0, const=0):
    x = np.linspace(mu-10*sigma_sq, mu+10*sigma_sq, num=200)
    y_gauss = (amp/np.sqrt(2*np.pi*sigma_sq))*np.exp(-0.5*(x-mu)**2/sigma_sq)
    y_linear = slope*x + const
    y = y_gauss + y_linear
    return x, y

# Gaussiand peak generation
amplitude = 6
center = 3884
variance = 4
slope = 0
intercept = 0.05
x, y = generate_gaussian(amplitude, center, variance, slope, intercept)


#Create a lmfit model: Gaussian peak + linear background
gaussian = GaussianModel()
background = LinearModel()
model = gaussian + background

#Find what model parameters you need to specify
print('parameter names: {}'.format(model.param_names))
print('independent variables: {}'.format(model.independent_vars))

#Model fit
result = model.fit(y, x=x, amplitude=3, center=3880,
                   sigma=3, slope=0, intercept=0.1)
y_fit = result.best_fit #the simulated intensity

result.best_values #the extracted peak parameters

# Comparison of the fitted spectrum with the original one
plt.plot(x, y, label='model spectrum')
plt.plot(x, y_fit, label='fitted spectrum')
plt.xlabel('wavelength, (angstroms')
plt.ylabel('intensity')
plt.legend()

Вывод:

parameter names: ['amplitude', 'center', 'sigma', 'slope', 'intercept']

independent variables: ['x']

result.best_values
Out[139]: 
{'slope': 2.261379140543626e-13,
 'intercept': 0.04999999912168238,
 'amplitude': 6.000000000000174,
 'center': 3883.9999999999977,
 'sigma': 2.0000000000013993}