scipy
Funktionen mit scipy.optimieren curve_fit

Einführung

Anpassen einer Funktion an Daten aus einem Histogramm

Angenommen, es gibt einen Peak von normal (Gauß'schen) verteilten Daten (Mittelwert: 3,0, Standardabweichung: 0,3) in einem exponentiell abklingenden Hintergrund. Diese Verteilung kann in wenigen Schritten mit curve_fit versehen werden:

1.) Importieren Sie die erforderlichen Bibliotheken.

2.) Definieren Sie die Fit-Funktion, die an die Daten angepasst werden soll.

3.) Holen Sie Daten aus dem Experiment oder generieren Sie Daten. In diesem Beispiel werden Zufallsdaten generiert, um den Hintergrund und das Signal zu simulieren.

4.) Fügen Sie das Signal und den Hintergrund hinzu.

5.) Passen Sie die Funktion mit curve_fit an die Daten an.

6.) (Optional) Zeichnen Sie die Ergebnisse und die Daten auf.

In diesem Beispiel sind die beobachteten y-Werte die Höhen der Histogramm-Bins, während die beobachteten x-Werte die Mitten der Histogramm-Bins ( binscenters ) sind. Es ist notwendig, den Namen der Fit-Funktion, die x-Werte und die y-Werte an curve_fit zu curve_fit . Außerdem kann mit p0 ein optionales Argument angegeben werden, das grobe Schätzungen für die Anpassungsparameter p0 . curve_fit gibt popt und pcov , wobei popt die Anpassungsergebnisse für die Parameter enthält, während pcov die Kovarianzmatrix ist, deren diagonale Elemente die Varianz der angepassten Parameter darstellen.

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