scipy
Fonctions d'ajustement avec scipy.optimize curve_fit

Introduction

Adapter une fonction aux données d'un histogramme

Supposons qu'il y ait un pic de données distribuées normalement (gaussiennes) (moyenne: 3,0, écart-type: 0,3) dans un contexte en décomposition exponentielle. Cette distribution peut être équipée de curve_fit en quelques étapes:

1.) Importez les bibliothèques requises.

2.) Définissez la fonction d'ajustement à adapter aux données.

3.) Obtenir des données d’expérience ou générer des données. Dans cet exemple, des données aléatoires sont générées afin de simuler le fond et le signal.

4.) Ajoutez le signal et l'arrière-plan.

5.) Ajuster la fonction aux données avec curve_fit.

6.) (Facultatif) Tracez les résultats et les données.

Dans cet exemple, les valeurs de y observées sont les hauteurs des tranches d'histogramme, tandis que les valeurs de x observées sont les centres des binscenters histogramme ( binscenters ). Il est nécessaire de passer le nom de la fonction fit, les valeurs x et les valeurs y à curve_fit . De plus, un argument facultatif contenant des estimations approximatives pour les paramètres d'ajustement peut être donné avec p0 . curve_fit renvoie popt et pcov , où popt contient les résultats d'ajustement pour les paramètres, tandis que pcov est la matrice de covariance, dont les éléments diagonaux représentent la variance des paramètres ajustés.

# 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()