Установка двух пиков с гауссом в python - PullRequest
Новая
       72

Установка двух пиков с гауссом в python

0 голосов
/ 03 апреля 2020

Curve_fit не подходит должным образом. Я пытаюсь согласовать экспериментальные данные с curve_fit. Данные импортируются из файла .txt в массив: 

d = np.loadtxt("data.txt")
data_x          = np.array(d[:, 0])
data_y          = np.array(d[:, 2])
data_y_err      = np.array(d[:, 3])

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

def model_dGauss(x, xc, A, y0, w, dx):
    P = A/(w*np.sqrt(2*np.pi))
    mu1 = (x - (xc - dx/3))/(2*w**2)
    mu2 = (x - (xc + 2*dx/3))/(2*w**2)    
    return y0 + P * ( np.exp(-mu1**2) + 0.5 * np.exp(-mu2**2))

Использование значений для догадки очень чувствительно к моим значениям угадывания. Где находится точка подбора данных, если только почти идеальный параметр подгонки даст результат? Или я делаю что-то совершенно не так? 

t = np.linspace(8.4, 10, 300)
guess_dG = [32, 1, 10, 0.1, 0.2]
popt, pcov = curve_fit(model_dGauss, data_x, data_y, p0=guess_dG, sigma=data_y_err, absolute_sigma=True)
A, xc, y0, w, dx = popt

Печать данных 

plt.scatter(data_x, data_y)
plt.plot(t, model_dGauss(t1,*popt))
plt.errorbar(data_x, data_y, yerr=data_y_err)

дает:

Результат печати

Результат - прямая линия внизу моего графика, в то время как оцениваемые параметры не так уж плохи. Как это может быть?

1 Ответ

0 голосов
/ 08 апреля 2020

Всегда приветствуется полный пример кода (и, как правило, здесь, на SO). Чтобы избежать путаницы в использовании curve_fit здесь, позвольте мне предположить, что вам будет проще использовать lmfit (https://lmfit.github.io/lmfit-py) и особенно его встроенные функции модели и использование именованных параметры. С помощью lmfit ваш код для двух гауссианов плюс постоянное смещение может выглядеть следующим образом: 

from lmfit.models import GaussianModel, ConstantModel

# start with 1 Gaussian + Constant offset:
model = GaussianModel(prefix='p1_') + ConstantModel()

# this model will have parameters named:  
# p1_amplitude, p1_center, p1_sigma, and c.  
# here we give initial values to these parameters
params = model.make_params(p1_amplitude=10, p1_center=32, p1_sigma=0.5, c=10)

# optionally place bounds on parameters (probably not needed here):
params['p1_amplitude'].min = 0.
## params['p1_center'].vary = False # fix a parameter from varying in fit

# now do the fit (including weighting residual by 1/y_err):
result = model.fit(data_y, params, x=data_x, weights=1.0/data_y_err)

# print out param values, uncertainties, and fit statistics, or get best-fit
# parameters from `result.params`
print(result.fit_report())

# plot results
plt.errorbar(data_x, data_y, yerr=data_y_err, label='data')
plt.plot(data_x, result.best_fit, label='best fit')
plt.legend()
plt.show()

Чтобы добавить второй гауссиан, вы можете просто сделать 

model = GaussianModel(prefix='p1_') + GaussianModel(prefix='p2_') + ConstantModel()

# and then:  
params = model.make_params(p1_amplitude=10, p1_center=32, p1_sigma=0.5, c=10, 
                           p2_amplitude=2, p2_center=31.75, p1_sigma=0.5)

и так далее.

Ваша модель имеет два гауссовых разделения или, по крайней мере, имеет "связанные" значения - значения sigma должны быть одинаковыми для двух пиков, а амплитуда второго равна половине амплитуды первого. Как определено до сих пор, 2-гауссова модель имеет все параметры, являющиеся независимыми. Но в lmfit есть механизм для установки ограничений на любой параметр, давая алгебраическое выражение c в терминах других параметров. Так, например, вы могли бы сказать, 

params['p2_sigma'].expr = 'p1_sigma'
params['p2_amplitude'].expr = 'p1_amplitude / 2.0'

Теперь, p2_amplitude и p2_sigma не будут независимо изменяться в подгонке, но будут ограничены, чтобы иметь эти значения.

...