In diesem kurzen Abschnitt lernen wir die wesentlichen Ideen bei Numpy kennen.
Eine längere und bessere (und englischere) Fassung davon sind die [Numpy fundamentals](https://numpy.org/devdocs/user/basics.html) im Numpy User Guide, die Sie am besten im Anschluss überfliegen sollten.
Ein richtiges (und sehr gutes) Lehrbuch für alle, die mit Python vertraut sind (und als das betrachten wir uns zu diesem Punkt der Vorlesung) ist [Nicolas Rougier's From Python to Numpy](https://www.labri.fr/perso/nrougier/from-python-to-numpy/), wo für uns zunächst Kapitel 3 relevant ist.
Kurz lässt sich sagen, dass mit Numpy [Array-orientierte Programmierung](https://en.wikipedia.org/wiki/Array_programming)(auch:*Vektorisierung*) in Python möglich wird.
## Arrays
<!--
Zur Probability mit Numpy ist das hier (auf deutsch) gut geeignet:
Es gibt auch eine `arange`-Methode, sie erzeugt aber keine Range-Objekte in Numpy:
%% Cell type:code id: tags:
``` python
myRange=np.arange(1,7)
print(myRange,type(myRange),myRange.dtype)
myRange
```
%% Output
[1 2 3 4 5 6] <class 'numpy.ndarray'> int64
array([1, 2, 3, 4, 5, 6])
%% Cell type:markdown id: tags:
Der Numpy-Array `ndarray` trägt im Gegensatz zur Python-Liste einen festen Datentyp, den alle Elemente gemeinsam haben. Das können alle Numpy-datatypes (dtypes) sein, z.B. `double` (kompatibel mit dem Python-`float`) oder `float32` (auf den meisten Plattformen kompatibel mit dem C-`float`) oder `long` (kompatibel mit dem Python-`int`).
Es gibt noch mehr Wege, Arrays zu erzeugen, außer mit `arange` oder durch konvertieren einer Python-Sequenz. Z.B. lässt sich mit `ones` ein Array gefüllt mit $1$ und mit `zeros` ein Array gefüllt mit $0$ erzeugen.
%% Cell type:markdown id: tags:
Der Grund dafür, dass das Array den Namen `ndarray` trägt, ist, dass es für \`\`$n$-dimensional array'' steht.
Wenn man in Python eine Matrix speichern möchte, würde man das als Liste der Zeilenvektoren (oder der Spaltenvektoren) tun, etwa
%% Cell type:code id: tags:
``` python
matrix=[[1,0,0],[0,1,0],[0,0,1]]# Einheitsmatrix
print(matrix,type(matrix))
quarkix=matrix# eine Kopie
print(quarkix[0][0])
quarkix[0][0]=0# wir ändern den oberen linken Eintrag
print(matrix[0][0])
```
%% Output
[[1, 0, 0], [0, 1, 0], [0, 0, 1]] <class 'list'>
1
0
%% Cell type:markdown id: tags:
Daran sehen wir ein Problem: Python behandelt unsere Matrix wie eine Liste (so haben wir es ja auch hingeschrieben), also wird beim kopieren der Liste der Inhalt (die Zeilenvektoren) nicht mitkopiert (sondern nur die Pointer darauf).
%% Cell type:code id: tags:
``` python
npmatrix=np.identity(3,int)# Einheitsmatrix
print(npmatrix,type(npmatrix))
npquarkix=npmatrix[:]
print(npquarkix[0][0])
npquarkix[0][0]=2# wir ändern den oberen linken Eintrag
print(npmatrix[0][0])
# Mit einer echten Kopie wäre das nicht passiert:
real_copy=npmatrix.copy()
real_copy[0][0]=1
assertnpmatrix[0][0]!=1
```
%% Output
[[1 0 0]
[0 1 0]
[0 0 1]] <class 'numpy.ndarray'>
1
2
%% Cell type:markdown id: tags:
Es ist wichtig, festzustellen, dass der Numpy-Array das gleiche Verhalten an den Tag legt wie unsere Python-Liste-von-Listen. Wir können gleich damit indizieren und slicen, und es gibt das gleiche Problem beim Kopieren über die Slicing-Syntax.
Der `shape`-Parameter sagt uns, welche Form unser Numpy-Array hat. Dabei handelt es sich um ein $d$-Tupel, wobei $d$ die Dimension ist. Eine Matrix ist $2$-dimensional, ein Vektor $1$-dimensional und ein Skalar $0$-dimensional.
%% Cell type:code id: tags:
``` python
print(npmatrix.shape)
print(npmatrix[0],npmatrix[0].shape)
print(npmatrix[0][0],npmatrix[0][0].shape)
```
%% Output
(3, 3)
[2 0 0] (3,)
2 ()
%% Cell type:markdown id: tags:
In Numpy kann man noch etwas feiner slicen.
Die Allgemeine Syntax ist `[start:stop:step, ..]` wobei man mit dem Komma getrennt über die Achsen geht. Ein zweidimensionaler Array hat zwei Achsen, wobei Achse $0$ von oben nach unten und Achse $1$ von links nach rechts indiziert ist. Während "step" auch mit Python-Listen funktioniert, ist das indizieren mit mehreren Achsen eine Spezialität von Numpy.
In klassischem Python-Code würden wir auf einen Eintrag einer Matrix zugreifen mit `matrix[x][y]`, und das funktioniert so auch in Numpy. Allerdings wird dabei zunächst ein weiteres Listenobjekt `matrix[x]` erzeugt (beim Slicing auch zusätzlicher Speicher dafür belegt) und dann darauf `[y]` aufgerufen. Es ist daher grundsätzlich effizienter, direkt Numpy's `[x,y]` zu verwenden.
Numpy erzeugt bewusst keine Kopien beim Slicing, sondern nur eine andere Sichtweise auf den gleichen Speicherbereich (daher auch das oben beobachtete Verhalten bei `[:]`). Ob zwei Arrays auf den gleichen Speicherbereich verweisen, lässt sich mit `np.may_share_memory` prüfen. Dabei bedeutet ein positives Ergebnis keineswegs, dass die Arrays voneinander abhängig sind - so verweisen die erste und die zweite Spalte einer Matrix auch auf den gleichen Speicherbereich, nämlich die ganze Matrix. Wenn man nun einen der beiden Vektoren ändert, bleibt der andere unverändert - die ganze Matrix aber ändert sich mit.
%% Cell type:code id: tags:
``` python
npmatrix=np.array(list(range(1,10))).reshape(3,3)
candidates=(npmatrix[0::2],npmatrix[1])
print("May share memory (but actually don't):",
candidates,np.may_share_memory(*candidates))
print(type([0,0,0]))# vor der Zuweisung
npmatrix[1]=[0,0,0]
print(type(npmatrix[1]))# nach der Zuweisung
print(npmatrix)# die ganze Matrix ist wie verändert
```
%% Output
May share memory (but actually don't): (array([[1, 2, 3],
[7, 8, 9]]), array([4, 5, 6])) True
<class 'list'>
<class 'numpy.ndarray'>
[[1 2 3]
[0 0 0]
[7 8 9]]
%% Cell type:markdown id: tags:
Ausführliche Informationen zum Slicing und Indizieren liefert [die Dokumentation](https://numpy.org/doc/stable/user/basics.indexing.html).
%% Cell type:markdown id: tags:
## Broadcasting
Während für Python-Listen der Additionsoperator die Listenkonkatenation ist, und damit die Multiplikation von Listen mit Skalaren definiert ist, ist die Multiplikation von zwei Listen undefiniert.
Für Numpy-Arrays sind deutlich mehr arithmetische Operationen verfügbar:
%% Cell type:code id: tags:
``` python
E=np.identity(3,int)
print(E,"= E")
A=np.ones((3,3),int)
A[0]=[0,0,0]
print(A,"= A")
print(E+A,"= E + A")
print(E*A,"= EA")
print((E+A)**2,"= (E+A)(E+A)")
```
%% Output
[[1 0 0]
[0 1 0]
[0 0 1]] = E
[[0 0 0]
[1 1 1]
[1 1 1]] = A
[[1 0 0]
[1 2 1]
[1 1 2]] = E + A
[[0 0 0]
[0 1 0]
[0 0 1]] = EA
[[1 0 0]
[1 4 1]
[1 1 4]] = (E+A)(E+A)
%% Cell type:markdown id: tags:
Wenn man das aufmerksam nachverfolgt, stellt man fest, dass diese Rechnungen keine Matrizenmultiplikationen sind,
sondern schlicht elementweise erfolgt sind - so sind die Operationen auf Arrays definiert.
Besonders tückisch ist dies:
%% Cell type:code id: tags:
``` python
v=np.ones(3,int)
print(v,"= v")
print(A*v,"= A*v (aber nicht die Matrixmultiplikation)")
```
%% Output
[1 1 1] = v
[[0 0 0]
[1 1 1]
[1 1 1]] = A*v (aber nicht die Matrixmultiplikation)
%% Cell type:markdown id: tags:
Um explizit mit Matrizenkalkül zu rechnen, hat man früher den Numpy-Datentyp `matrix` verwendet, aber dieser ist als `deprecated` (veraltet) markiert und wird in zukünftigen Numpy-Versionen abgeschafft. Heutzutage nutzt man die Methode `np.matmul` oder den Infix-Operator `@`.
%% Cell type:code id: tags:
``` python
print(np.matmul(A,v))
print(A@v)
```
%% Output
[0 3 3]
[0 3 3]
%% Cell type:markdown id: tags:
Für viele Probleme ist es sehr hilfreich, nicht den Matrizenkalkül zu verwenden, sondern elementweise arithmetische Operationen auszuführen.
*Broadcasting* ist ein Mechanismus, der diesen elementweisen Kalkül etwas praktischer macht.
So ist die Operation ($n \times n$-array) * ($n$-Vektor) automatisch interpretiert, indem der $n$-Vektor $n$-fach kopiert wird, sodass die Multiplikation einer jeden Zeile des linken Arrays mit dem Vektor (elementweise) durchgeführt wird.
Dazu ist es wirklich hilfreich, einmal [die Dokumentation](https://numpy.org/doc/stable/user/basics.broadcasting.html) zu überfliegen.
%% Cell type:markdown id: tags:
## ufuncs
`ufunc` steht für "universal function" und bezeichnet eine Methode, die auf Numpy Arrays vektorisiert laufen kann.
Mit jeder `ufunc` lässt sich z.B. `reduce` durchführen, wo entlang einer Achse des Arrays die `ufunc` auf die resultierenden kleineren Arrays angewandt wird.
Um selbst eine `ufunc` zu schreiben, muss man [C-Code programmieren](https://numpy.org/doc/stable/user/c-info.ufunc-tutorial.html) oder aber [einen Wrapper um eine Python-Methode legen](https://numpy.org/doc/stable/reference/generated/numpy.frompyfunc.html).
Es lohnt sich, einen kurzen Blick auf alle bereits definierten `ufunc`s zu werfen:
Damit haben wir nun das Wahrscheinlichkeitsmaß. Wenn wir eine Stichprobe ziehen wollen, müssen wir noch irgendwo den "Zufall" her bekommen.
## Zufall importieren
Python bietet mit dem `random`-Modul eine Schnittstelle zu Pseudozufallszahlen. Die Methode `random.random` ist ein direkt in C implementierter [Mersenne Twister](https://de.wikipedia.org/wiki/Mersenne-Twister). Wenn man "echte" Zufallszahlen braucht, etwa für kryptografische Zwecke, gibt es dazu das `secrets`-Modul. Den Seed für den Mersenne Twister kann man angeben, und sollte man auch, um Zufallssimulationen reproduzierbar zu machen.
Damit man solche Transformationen nicht andauernd programmieren muss, kann man hier auch auf `random.randint(1,6)` oder auch auf `random.choice(range(1,7))` oder `random.randrange(1,7)` zurückgreifen. Dabei ist `randint` ein Kürzel für das entsprechende `randrange` und `choice` ist etwas allgemeiner.
Wir wollen aber festhalten: gegeben eine gleichverteilte "Zufallsvariable" $X=$`random.random` mit Werten in $[0,1)$ haben wir eine Abbildung $t =$`transform_unit_to_dice` konstruiert und implementiert, die Werte in $\Omega := \{1,2,3,4,5,6\}$ hat und $t(X)$ ist gleichverteilt. Die mathematische Abbildung $t \colon [0,1) \to \Omega$ ist eine Zufallsvariable, wir behandeln die [Verknüpfung](https://de.wikipedia.org/wiki/Komposition_(Mathematik)) $t \circ X$ als Zufallsvariable, die den Würfel modelliert.
Nun könnte man sich beschweren: `random.random()` nimmt gar keinen Parameter, ist also keine mathematische Abbildung von einem Definitionsbereich in die Menge $[0,1)$. Tatsächlich müssen wir uns vorstellen, dass es eine Abbildung $X \colon \Omega' \to [0,1)$ ist, und auf $\Omega'$ ein irgendwie geartetes Wahrscheinlichkeitsmaß definiert ist, sodass durch $X$ auf $[0,1)$ die Gleichverteilung induziert wird. Die Menge $\Omega'$ spielt für uns keine konkrete Rolle - da "kommt der Zufall her" und in der Notation `random.random()` sehen wir schon, dass wir eben kein konkretes Element von $\Omega'$ einsetzen, sondern pseudozufällig eins ziehen und das in $X$ einsetzen.
Die Methode `random.sample(population, k)` erlaubt es eine Stichprobe der Größe $k$ aus einer Population (einer Urne) zu ziehen - ziehen mit Zurücklegen. Für $k=1$ entspricht das einer Gleichverteilung auf der Population.
Mit der Methode `random.choices(population, weights=None, *, cum_weights=None, k=1)` kann man das ziehen aus der Population $k$-mal sampeln (und dabei anstelle einer geeigneten Population auch Gewichte vergeben).
