scipy
Функции фитинга с scipy.optimize curve_fit
Поиск…
Вступление
В научных приложениях часто требуется функция, которая описывает ожидаемое появление точек данных для реальных данных. Возможным оптимизатором для этой задачи является curve_fit из scipy.optimize. Ниже приводится пример применения curve_fit.
Установка функции на данные с гистограммы
Предположим, что на экспоненциально затухающем фоне имеется пик нормальных (гауссовских) распределенных данных (среднее: 3.0, стандартное отклонение: 0,3). Этот дистрибутив может быть установлен с помощью функции curve_fit в течение нескольких шагов:
1.) Импортируйте необходимые библиотеки.
2.) Определите функцию подгонки, которая должна быть привязана к данным.
3.) Получить данные из эксперимента или генерировать данные. В этом примере генерируются случайные данные для имитации фона и сигнала.
4.) Добавьте сигнал и фон.
5.) Установите функцию на данные с помощью curve_fit.
6.) (необязательно). Выделите результаты и данные.
В этом примере наблюдаемые значения y представляют собой высоты бункеров гистограммы, а наблюдаемые значения x являются центрами гистограмм ( binscenters ). Необходимо передать имя функции fit, значения x и значения y в curve_fit . Кроме того, необязательный аргумент, содержащий приблизительные оценки для параметров подгонки, может быть задан с p0 . curve_fit возвращает popt и pcov , где popt содержит результаты соответствия для параметров, а pcov - ковариационная матрица, диагональные элементы которой представляют собой дисперсию установленных параметров.
# 1.) Necessary imports.
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# 2.) Define fit function.
def fit_function(x, A, beta, B, mu, sigma):
return (A * np.exp(-x/beta) + B * np.exp(-1.0 * (x - mu)**2 / (2 * sigma**2)))
# 3.) Generate exponential and gaussian data and histograms.
data = np.random.exponential(scale=2.0, size=100000)
data2 = np.random.normal(loc=3.0, scale=0.3, size=15000)
bins = np.linspace(0, 6, 61)
data_entries_1, bins_1 = np.histogram(data, bins=bins)
data_entries_2, bins_2 = np.histogram(data2, bins=bins)
# 4.) Add histograms of exponential and gaussian data.
data_entries = data_entries_1 + data_entries_2
binscenters = np.array([0.5 * (bins[i] + bins[i+1]) for i in range(len(bins)-1)])
# 5.) Fit the function to the histogram data.
popt, pcov = curve_fit(fit_function, xdata=binscenters, ydata=data_entries, p0=[20000, 2.0, 2000, 3.0, 0.3])
print(popt)
# 6.)
# Generate enough x values to make the curves look smooth.
xspace = np.linspace(0, 6, 100000)
# Plot the histogram and the fitted function.
plt.bar(binscenters, data_entries, width=bins[1] - bins[0], color='navy', label=r'Histogram entries')
plt.plot(xspace, fit_function(xspace, *popt), color='darkorange', linewidth=2.5, label=r'Fitted function')
# Make the plot nicer.
plt.xlim(0,6)
plt.xlabel(r'x axis')
plt.ylabel(r'Number of entries')
plt.title(r'Exponential decay with gaussian peak')
plt.legend(loc='best')
plt.show()
plt.clf()
