scipy
Passingsfuncties met scipy.optimize curve_fit
Zoeken…
Invoering
Het aanpassen van een functie die het verwachte voorkomen van gegevenspunten aan echte gegevens beschrijft, is vaak vereist in wetenschappelijke toepassingen. Een mogelijke optimizer voor deze taak is curve_fit van scipy.optimize. In het volgende wordt een voorbeeld van toepassing van curve_fit gegeven.
Een functie aanpassen aan gegevens uit een histogram
Stel dat er een piek is van normaal (Gaussiaanse) gedistribueerde gegevens (gemiddelde: 3,0, standaardafwijking: 0,3) op een exponentieel vervallende achtergrond. Deze verdeling kan in enkele stappen worden uitgerust met curve_fit:
1.) Importeer de vereiste bibliotheken.
2.) Definieer de aanpassingsfunctie die moet worden aangepast aan de gegevens.
3.) Gegevens verkrijgen uit experiment of gegevens genereren. In dit voorbeeld worden willekeurige gegevens gegenereerd om de achtergrond en het signaal te simuleren.
4.) Voeg het signaal en de achtergrond toe.
5.) Pas de functie aan op de gegevens met curve_fit.
6.) (Optioneel) Plot de resultaten en de gegevens.
In dit voorbeeld zijn de waargenomen y-waarden de hoogten van de histogrambins, terwijl de waargenomen x-waarden de centra van de histogrambins ( binscenters ) zijn. Het is noodzakelijk om de naam van de aanpassingsfunctie, de x-waarden en de y-waarden door te curve_fit aan curve_fit . Bovendien kan een optioneel argument met ruwe schattingen voor de fitparameters worden gegeven met p0 . curve_fit retourneert popt en pcov , waarbij popt de passende resultaten voor de parameters bevat, terwijl pcov de covariantiematrix is, waarvan de diagonale elementen de variantie van de gepaste parameters vertegenwoordigen.
# 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()
