Я пытаюсь разложить сложные сигналы газовой хроматограммы на отдельные гауссовские сигналы. Вот пример, где пунктирная линия представляет сигнал, который я пытаюсь развернуть.
Мне удалось написать код для этого с помощью scipy.optimize.curve_fit; однако, однажды примененные к реальным данным, результаты были ненадежными. Я считаю, что возможность устанавливать границы для моих параметров улучшит мои результаты, поэтому я пытаюсь использовать lmfit, который позволяет это. У меня проблема с тем, чтобы lmfit работал с переменным количеством параметров. Сигналы, с которыми я работаю, могут иметь произвольное количество базовых гауссовских компонентов, поэтому количество необходимых мне параметров будет варьироваться. Я нашел здесь подсказки, но до сих пор не могу понять ...
Создание модели lmfit Python с произвольным числом параметров а>
Вот код, с которым я сейчас работаю. Код будет запущен, но оценки параметров не изменятся, когда модель подходит. Кто-нибудь знает, как заставить работать мою модель?
import numpy as np
from collections import OrderedDict
from scipy.stats import norm
from lmfit import Parameters, Model
def add_peaks(x_range, *pars):
y = np.zeros(len(x_range))
for i in np.arange(0, len(pars), 3):
curve = norm.pdf(x_range, pars[i], pars[i+1]) * pars[i+2]
y = y + curve
return(y)
# generate some fake data
x_range = np.linspace(0, 100, 1000)
peaks = [50., 40., 60.]
a = norm.pdf(x_range, peaks[0], 5) * 2
b = norm.pdf(x_range, peaks[1], 1) * 0.1
c = norm.pdf(x_range, peaks[2], 1) * 0.1
fake = a + b + c
param_dict = OrderedDict()
for i in range(0, len(peaks)):
param_dict['pk' + str(i)] = peaks[i]
param_dict['wid' + str(i)] = 1.
param_dict['mult' + str(i)] = 1.
# In case, you'd like to see the plot of fake data
#y = add_peaks(x_range, *param_dict.values())
#plt.plot(x_range, y)
#plt.show()
# Initialize the model and fit
pmodel = Model(add_peaks)
params = pmodel.make_params()
for i in param_dict.keys():
params.add(i, value=param_dict[i])
result = pmodel.fit(fake, params=params, x_range=x_range)
print(result.fit_report())