NumPy Tutorial

Code & Zahlen

NumPy ist eine Bibliothek für Python. Ich empfehle Python 3 zu benutzen. Es werden unter Ubuntu folgende Pakete benötigt:

  • python3-matplotlib
  • python3-numpy
  • python3-scipy

Offizielle Dokumentation:

Beispiel: Anpassen einer Gaußfunktion

Gegeben sei die Datei gauss.csv in der $x$ und $y$ Werte stehen. Dies sind die ersten paar Zeilen der Datei:

-5.000000000000000000e+00 6.500006597214588178e-01
-4.989989989989989994e+00 5.909619296200073757e-01
-4.979979979979979987e+00 5.267851553980463786e-01
-4.969969969969969981e+00 4.393687797793060512e-01
-4.959959959959959974e+00 8.199626613472296155e-01

Dazu beginne ich eine neue Textdatei, auswertung.py:

import matplotlib.pyplot as pl
import numpy as np
import scipy.optimize as op

def normal_distribution(x, integral, mean, sigma):
    return integral / np.sqrt(2 * np.pi) / sigma * np.exp(- x**2 / sigma**2)

Mit import numpy as np importierte ich die NumPy Bibliothek unter dem Kurznamen np. Das gleiche mache ich ebenfalls mit der Plotbibliothek matplotlib.

Danach definiere ich die Funktion normal_distribution die einer Gaußfunktion entspricht:

$$N(x; I, \mu, \sigma) := \frac{I}{\sqrt{2\pi} \sigma} \exp\left(-\frac{x^2}{\sigma^2}\right)$$

Dieses Programm kann immer mit python3 auswertung.py ausgeführt werden.

Die Daten können nun wie folgt eingelesen werden:

data = np.loadtxt('gauss.csv')

Dies erstellt ein numpy.array, ein Objekt in dem Arrays beliebiger Größe gespeichert werden können, sowohl 1D als auch 2D. Daraus müssen nun $x$ und $y$ extrahiert werden:

x = data[:, 0]
y = data[:, 1]

Die Syntax mit den […] ist ähnlich wie in C auch. In NumPy ist die Syntax allerdings erweitert, so dass direkt Zeilen und Spalten ausgewählt werden kann. Der erste Index ist die Zeile, der zweite die Spalte. Für alle Elemente gibt es das Symbol :. Die oberen zwei Zeilen Quelltext nehmen also die erste (Index 0) und zweite (Index 1) Spalte aus der Datenmatrix heraus.

Die Funktion curve_fit passt eine beliebige Funktion an die Daten an:

popt, pconv = op.curve_fit(normal_distribution, x, y)

popt ist ein Vektor (Array) mit den Fitparametern, pconv ist die Kovarianzmatrix, aus der die Fehler für den Fit berechnet werden können.

Hier ist popt[0] also die Amplitude der Gaußfunktion, popt[1] der Mittelwert und popt[2] die Standardabweichung.

Für die Fehler der Parameter braucht man noch eine Zeile mehr:

fehler = np.sqrt(pconv.diagonal())

Dies nimmt die Diagonalen Elemente und zieht die Wurzel aus diesen. Das Array fehler enthält in der gleichen Reihenfolge die Fehler.

Das ganze soll jetzt noch in einer Grafik dargestellt werden. Dazu wird matplotlib benutzt:

pl.plot(x, y)

Die Grafik kann jetzt interaktiv angezeigt werden:

pl.show()

Oder gespeichert werden:

pl.savefig('gauss.pdf')

Es können alle üblichen Dateiformate wie .png oder .pdf benutzt werden.

Dies hat jetzt nur die Originaldaten geplottet. Es soll noch die angepasste Funktion geplottet werden. NumPy kann keine Funktion direkt plotten. Man muss erst $x$-Punkte erzeugen, dann die Funktion darauf anwenden um $y = f(x)$ Werte zu erzeugen und dies dann plotten.

fitted_x = np.linspace(np.min(x), np.max(x), 1000)
fitted_y = normal_distribution(fitted_x, *popt)

Diese Zeilen bedürfen einiger Erklärung: linspace erzeugt linear gleichverteilte Punkte zwischen zwei Punkten. Die ersten beiden Parameter sind Start- und Endpunkt, der dritte Parameter die Anzahl der Stützstellen. Ich wähle immer das Minimum und Maximum meiner Eingangsdaten als Endpunkte, auf diese Weise ist der Definitionsbereich genau gleich. 1000 Stützstellen geben immer schön glatte Kurven, 100 ist oft zu wenig. Dies wird als Vektor in fitted_x gespeichert.

Dann wird die Funktion normal_distribution auf den ganzen Vektor fitted_x elementweise angewandt. Es bedarf hier keiner for-Schleife, da NumPy dies im Hintergrund automatisch macht. Das *popt ist kein Pointer, sondern bedeutet in Python, dass die Liste popt „ausgeschüttet" werden soll und die Einträge als einzelne Parameter benutzt werden sollen. Die Funktion normal_distribution erwartet ja vier Parameter, popt enthält die letzten drei. Somit passt das dann.

Dies soll in einen schönen Plot:

# Plotte die Originaldaten. Da es Messdaten sind, werden sie nicht mit
# einer Linie verbunden.
pl.plot(x, y, linestyle='none', marker='.', label='Daten')

# Plotte die Anpassungsfunktion. Diesmal ohne Punkte, aber mit einer Linie.
pl.plot(fitted_x, fitted_y, linestyle='-', marker='none', label='Anpassung')

# Aktiviere das Gitter.
pl.grid(True)

# Setze eine Legende an den optimalen Platz.
pl.legend(loc='best')

# Speichere den Plot als PDF.
pl.savefig('anpassung.pdf')

# Setze die Plot-Umgebung zurück. Dies ist wichtig, weil mehrere ``plot()``
# Befehle in einen Plot gehen, wie oben schon benutzt.
pl.clf()

Wichtige Funktionen

NumPy

Diese Funktionen können mit np. benutzt werden.

array([1, 2, 3])

Erzeugt aus einer Liste ein NumPy Array.

linspace(min, max, anzahl)

Erzeugt eine Liste mit Stützstellen.

sqrt(x), exp(x), pi, sin(x)

Normale Mathematikfunktionen, die direkt auf NumPy Arrays angewandt werden können und so elementweise angewandt.

Matplotlib

Diese Funktionen können mit pl. benutzt werden.

plot(x, y)

Plottet einfache Punkte. Es gibt noch sehr viele Extraoptionen, die man noch angeben kann, zum Beispiel linestyle, marker, label. Siehe plot Dokumentation.

errorbar(x, y, yerr, xerr)

Plottet Punkte mit Fehlerbalken. errorbar Dokumentation

set_xscale('log')

Setzt die $x$-Achse auf logarithmisch.

set_yscale('log')

Setzt die $y$-Achse auf logarithmisch.

legend(loc='best')

Aktiviert die Legende.