Commit cf55c205 authored by Konrad Völkel's avatar Konrad Völkel
Browse files

+expectation maximization (& micro fixes)

parent 29c9b909
Pipeline #92948 passed with stages
in 5 minutes and 28 seconds
......@@ -63,4 +63,5 @@ parts:
- file: lineare-regression
- file: kernel-methoden
- file: svm
- file: model-selection
\ No newline at end of file
- file: model-selection
- file: expectation-maximization
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
%% Cell type:markdown id:66aa7da7-bfa4-44d9-8d1f-b289931187f1 tags:
## Iris
# Iris
Mit dem Iris-Datensatz werden wir uns noch häufiger beschäftigen, auch in den Übungsaufgaben.
Hier sehen wir uns erst einmal eine Plot-Technik "jittering" an, die es ermöglicht Daten mit einem Scatterplot darzustellen, die häufig den selben Wert annehmen. Das ist eine Alternative zu Histogramm-Techniken sowie einer Anpassung der Transparenz (des Alpha-Werts).
%% Cell type:code id:db48fff4-f9c6-4a19-95a6-7dc177351a15 tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
iris = datasets.load_iris()
x = iris.data
classes = iris.target
print(x[:5])
print(classes[:5])
```
%% Output
[[5.1 3.5 1.4 0.2]
[4.9 3. 1.4 0.2]
[4.7 3.2 1.3 0.2]
[4.6 3.1 1.5 0.2]
[5. 3.6 1.4 0.2]]
[0 0 0 0 0]
%% Cell type:code id:34796b74-d57a-499f-ab9b-c561cb1916cd tags:
``` python
def jitter(xs):
"""jitter your values if they are too often the same"""
b = np.diff(np.sort(xs))
distance = b[b>0].min()
return xs + distance * (2*np.random.random(len(xs)) - 1)
fig, (ax1,ax2) = plt.subplots(2, figsize=(20,20))
petlen, petwid = jitter(x[:, 2]), jitter(x[:, 3]) # jitter!
#petlen, petwid = x[:, 2], x[:, 3] # no jitter here, it looks like a grid <---- Probieren Sie den Unterschied!
plmin, plmax = np.min(petlen), np.max(petlen)
cmap = plt.cm.cividis_r
ax1.scatter(petlen, petwid, c=classes, cmap=cmap, alpha=0.7, s=200)
for c in sorted(list(set(classes))):
ax2.hist(petlen[classes == c], range=(plmin,plmax), bins=50, color=cmap(c/2-0.1))
ax2.set_xlabel("Petal length")
ax2.set_ylabel("Petal width")
ax1.set_ylabel("Petal width")
plt.show()
```
%% Output
......
......@@ -1515,7 +1515,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
......@@ -1529,7 +1529,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.2"
"version": "3.8.10"
}
},
"nbformat": 4,
......
%% Cell type:markdown id:5fa2e75d-6697-4c5a-84f4-437f76847c86 tags:
### Konfidenzintervalle mit Python
## Konfidenzintervalle mit Python
%% Cell type:code id:9616c62b-f52b-464b-9e9a-da8ee864646a tags:
``` python
import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
```
%% Cell type:markdown id:c4807851-d434-4534-a7e5-faf35b3b1629 tags:
Um Intervalle zu plotten, kann man mit `fill_between` arbeiten:
%% Cell type:code id:852cf532-a5e0-42ae-811f-21718215b18c tags:
``` python
x = np.arange(0,20,0.1)
y = x**2
ci = 40
fig, ax = plt.subplots()
ax.plot(x, y)
ax.fill_between(x, (y-ci), (y+ci), alpha=0.2, color='black')
fig.suptitle("Intervall um die Werte von x²")
plt.savefig("images/fiktives-konfidenzintervall.png")
plt.show()
```
%% Output
%% Cell type:code id:f9e00974-6d0d-4355-9a76-d63dffafb36e tags:
``` python
x = st.norm.rvs(5, 4, 1000)
fig, axs = plt.subplots(1,2)
fig.suptitle("Darstellungen einer Stichprobe einer Normalverteilung")
axs[0].boxplot(x)
axs[0].set_title("Boxplot")
axs[1].hist(x, orientation='horizontal')
axs[1].set_title("Histogramm")
plt.savefig("images/boxplot-histogramm-vergleich.png")
plt.show()
```
%% Output
%% Cell type:markdown id:ce5211d6-25cb-4443-a01b-55e7002ace99 tags:
Achtung: diese Quantile sind eben *keine* Konfidenzintervalle.
Um für eine Stichprobe einer Normalverteilung ein Konfidenzintervall berechnen zu lassen, können wir Scipy benutzen:
%% Cell type:code id:b7d0ec37-0892-4df2-ae7e-deff140e500c tags:
``` python
np.random.seed(1234321)
# normalverteilte Stichprobe generieren:
real_mu, real_sigma = 5, 4
n = 25
x = st.norm.rvs(real_mu, real_sigma, n)
# Erwartungstreue Punktschätzer:
mean, sigma = x.mean(), x.std(ddof=1)
# Intervall aus der t-Verteilung:
ci = st.t.interval(0.95, n-1, loc=mean,
scale=st.sem(x))
from matplotlib.patches import Rectangle
fig, ax = plt.subplots(figsize=(10,5))
fig.suptitle("Darstellungen eines Konfidenzintervalls des Erwartungswerts einer normalverteilten Stichprobe der Größe "+str(n))
ax.scatter(x, np.full(n, 1), s=2, label="Daten")
ax.hist(x, bins=20, alpha=0.5, label="Häufigkeit")
ax.plot(mean, 1, "+",
color='black', label="Geschätzter Mittelwert")
ax.plot(real_mu, 1, "x", color='black', label="Tatsächlicher Erwartungswert")
ax.add_patch(Rectangle((ci[0], 0), ci[1]-ci[0], 5,
alpha=0.5, color='red',
label="Konfidenzintervall zu α=0.05"))
ax.legend()
plt.savefig("images/konfidenzintervall-normalverteilung.png")
plt.show()
```
%% Output
%% Cell type:markdown id:6e88ac32-4c15-4d65-bd0d-5d42cc0dd73d tags:
Alternativ kann man auch `errorbar` benutzen:
%% Cell type:code id:f43f67f8-85b8-479c-a9a2-c2e6584c4571 tags:
``` python
fig, ax = plt.subplots(figsize=(10,2))
ax.scatter(x, np.full(n, 1), s=2, label="Daten")
ax.set_yticklabels([])
ax.plot(mean, 1, "+",
color='black', label="Geschätzter Mittelwert")
ax.plot(real_mu, 1, "x", color='black', label="Tatsächlicher Erwartungswert")
ax.errorbar(mean, 1, xerr=(ci[1]-ci[0])/2, yerr=None,
label="Konfidenzintervall zu α=0.05",
color="red", alpha=0.5)
ax.legend()
plt.show()
```
%% Output
%% Cell type:markdown id:c4807f35-f019-44e1-889f-50c79e23d499 tags:
Das Python-Modul `statsmodels` liefert noch mehr Möglichkeiten, Konfidenzintervalle zu berechnen, und mit `seaborn` lassen sich diese komfortabel plotten.
Anstatt Konfidenzintervalle zu visualisieren, werden diese meist implizit genutzt, bei der Konstruktion von Tests.
......
## Konfidenzintervalle
# Konfidenzintervalle
Wir betrachten zur Motivation die Grafik {ref}`warming-without-confidence` und setzen Sie danach in Bezug zu {ref}`warming-with-confidence`
......
%% Cell type:markdown id:e1e42003-c0b9-4e9f-8ac0-d350876ba047 tags:
## LDA (linear discriminant analysis)
# LDA (linear discriminant analysis)
Im Gegensatz zur PCA ist die LDA ein supervised learning Verfahren, d.h. es muss eine Klassifikation gegeben sein.
Man kann sich vorstellen, dass LDA mit einem gelabelten (d.h. klassifizierten) Trainings-Datensatz gefüttert wird,
dann eine Klassifikation ausspuckt, die man auch auf ungelabelte (d.h. noch nicht klassifizierte) Test-Daten anwenden kann.
Die Art der Klassifikation sind lineare Grenzen zwischen den Klassen.
Bevor wir LDA erklären, stellen wir uns etwas einfacheres vor:
Von bereits klassifizierten Daten berechnen wir Zentroide (jedes Zentroid ist der Mittelwert aller Punkte innerhalb einer Klasse),
und sehen die Menge der Zentroide nun als eigenen Datensatz an,
auf den wir PCA anwenden können.
Die erste Hauptkomponente ist dann die Richtung der größten Varianz innerhalb der Zentroide.
Diese Koordinaten erleichtern es, die Klassen linear voneinander zu trennen.
Die Koordinatentransformation, die sich bei der Zentroide-PCA ergibt, lässt sich auch auf die ursprünglichen Daten anwenden.
Die erste Hauptkomponente ist dann eine Richtung mit großer Varianz zwischen den Klassen, während die Varianz innerhalb der Klassen ignoriert wird.
Dieses Verfahren ist der LDA sehr ähnlich; der wesentliche Unterschied ist nun, dass bei der LDA direkt die Richtung größter Varianz zwischen den Klassen berechnet wird.
Grundvoraussetzung bei der LDA ist, dass die Klassen näherungsweise normalverteilt sind (mit gleicher Kovarianzmatrix $\Sigma$ für jede Klasse, man spricht von *Homoskedastizität*). Wenn die Voraussetzung der Homoskedastizität nicht gegeben ist, kann ein verallgemeinertes Verfahren, die quadratische Diskriminanzanalyse, QDA, verwendet werden - die besprechen wir hier nicht.
### LDA für zwei Klassen
## LDA für zwei Klassen
Wenn es genau zwei Klassen gibt, ist das Ergebnis der LDA eine affine Hyperebene (d.h. ein Unterraum eine Dimension kleiner als der ursprüngliche Vektorraum, der den Ursprung nicht enthalten muss), die die beiden Klassen optimal separiert. Äquivalent (dual) dazu ist das Ergebnis der LDA ein Vektor mit Aufpunkt (ein Normalenvektor der Hyperebene an einer Stelle) sodass die Klasseneinteilung durch Projektion auf die so beschriebene Gerade entschieden werden kann (links vom Aufpunkt = Klasse 1, rechts vom Aufpunkt = Klasse 2). Zur Erinnerung: ein Normalenvektor zu einem Untervektorraum ist ein Vektor, der auf allen Vektoren des Untervektorraums orthogonal ist (= einen rechten Winkel hat), der also gewissermaßen in eine andere Dimension zeigt.
Im konkreten Fall, wo die Daten dreidimensional sind, ist die Hyperebene eine gewöhnliche Ebene im dreidimensionalen Raum.
Im konkreten Fall, wo die Daten zweidimensional sind, ist die Hyperebene eine gewöhnliche Gerade in der Ebene.
In Abbildung {ref}`lda-example` (aus {cite}`esl`) ist ein Beispiel zweidimensionaler Daten, die in zwei normalverteilte Klassen fallen. Die Projektion auf die y-Achse im linken Bild erlaubt nicht so eine gute Unterscheidung wie die Projektion auf den Normalenvektor zur Gerade, die am besten die beiden Klassen separiert.
```{figure} images/sensible_rule.png
:width: 480px
:align: center
:name: lda-example
Lineare Diskriminanzanalyse (rechts)
```
Ein Algorithmus zur Berechnung der LDA geht so:
1. Wir bestimmen eine Kovarianzmatrix, die für jeden Cluster die Kovarianz beschreiben soll (wenn wir für jeden Cluster einzeln die Kovarianzmatrix berechnen, werden wir unterschiedliche bekommen - wir wollen aber nur genau eine, da wir von Homoskedastizität ausgehren). Ein Weg ist, jeden Cluster einzeln zu zentrieren und für die zentrierten Daten gemeinsam die Kovarianzmatrix zu berechnen.
2. Führe eine Hauptkomponententransformation durch mit der zuvor bestimmten Kovarianzmatrix - danach sollten die Daten innerhalb jedes Clusters standardnormalverteilt aussehen.
3. Bestimme die Mittelpunkte der Cluster im so transformierten Koordinatensystem
4. Verbinde die Mittelpunkte mit einer Geraden. Abhängig davon, welche Klasse mehr Punkte enthält, liegt der Entscheidungspunkt auf dieser Geraden genau in der Mitte oder näher an einem der beiden Clustermittelpunkte. Der Entscheidungspunkt ist der Punkt, an dem man eine Senkrechte zu der Mittelpunkte-verbindenden Geraden einzeichnen kann, sodass die beiden Seiten der Senkrechten möglichst gut den Klassen entsprechen.
Manchmal wird dies auch anders formuliert, und anstelle einer Hauptachsentransformation werden die Abstände in Schritt 4 nicht dem euklidischen Abstand gemessen, sondern mit dem **Mahalanobis-Abstand**, das ist genau der euklidische Abstand im transformierten Raum. Im nicht-transformierten Raum sind die Cluster nicht sphärisch, daher ist die separierende Hyperebene dann auch nicht orthogonal zur Mittelpunkt-verbindenden Geraden.
Um daraus analog zur PCA eine Methode zu machen, die das Koordinatensystem so transformiert, dass man die hinteren, weniger wichtigen Dimensionen wegwerfen kann, führt man die LDA auf der Hyperebene erneut durch, iterativ, bis man genug Richtungen hat. Die verbleibende Hyperbene wirft man dann weg.
### LDA für mehr als zwei Klassen
## LDA für mehr als zwei Klassen
Um LDA auf mehr als zwei Klassen zu verallgemeinern, müssen wir die Gerade, die die beiden Clustermittelpunkte verbunden hat, auf die Situation mit mehreren Mittelpunkten verallgemeinern. Dabei ergibt sich im Allgemeinen ein Unterraum der Dimension $k-1$, wenn es $k$ Klassen gibt (machen Sie sich das für $k=2,3,4$ klar!).
Wenn $\mu_i$ die Mittelpunkte der Cluster sind, und $\mu$ der Mittelpunkt der $\mu_i$, so berechne die Kovarianzmatrix der Clustermittelpunkte:
$$
\Sigma_b = \frac{1}{k} \sum_{i=1}^k \left(\mu_i - \mu \right) {\left(\mu_i - \mu \right)}^T
$$
Wenn $\Sigma$ die Kovarianzmatrix innterhalb eines Clusters ist (die wir ja alle gleich annehmen), dann bestimmen wir für $\Sigma^{-1}\Sigma_b$ die $k-1$ größten Eigenwerte und zugehörige Eigenvektoren. Der davon aufgespannte Untervektorraum entspricht der Verbindungsgraden zweier Mittelpunkte im Fall von zwei Klassen.
Wenn $\Sigma^{-1}\Sigma_b w = \alpha w$ eine entsprechende Eigenwertgleichung ist,
so gilt also $\Sigma_b w = \alpha \Sigma w$, d.h. der Vektor $w$ wird von $\Sigma_b$ und von $\Sigma$ gleich transformiert.
### LDA mit Scikit-learn
## LDA mit Scikit-learn
Wir basteln erstmal synthetische Daten, ähnlich wie zuvor bei dem Experiment des Standardskalierens, PCA und k-Means:
%% Cell type:code id:e99721c5-12a9-4fd2-aa98-0528450bf84b tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(1234)
# craft some dummy data:
dummycenters = np.asarray([[10, -1, -1],
[12, 3, 0.3],
[15, 2, 0.2],
[20, 4, 0.1]])
dummycov = np.asarray([[3, 2, 1.5],
[2, 3, 1],
[1.5, 1, 1]])
dummydata = np.concatenate([np.random.multivariate_normal(center,
dummycov,
size=200)
for center in dummycenters])
original_labels = np.concatenate([np.full(200, k)
for k, _ in enumerate(dummycenters)])
dummymu = np.mean(dummydata, axis=0)
```
%% Cell type:code id:8acc8544-c2ae-4e59-87ca-5433d5b47e8b tags:
``` python
# let's take a look:
fig, ax = plt.subplots(figsize=(8,8))
ax.scatter(dummydata[:,0], dummydata[:,1], c=original_labels)
ax.plot(dummymu[0], dummymu[1],
"X", ms=15, alpha=0.9, c='red')
ax.plot(dummycenters[:,0], dummycenters[:,1],
"X", ms=15, alpha=0.9, c='black')
plt.show()
```
%% Output
%% Cell type:code id:b751a776-ea54-4eec-a5db-e085610708b3 tags:
``` python
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
LDA = LinearDiscriminantAnalysis(solver="svd", store_covariance=True)
lda_fit = LDA.fit(dummydata, original_labels)
new_labels = lda_fit.predict(dummydata)
lda_data = lda_fit.transform(dummydata)
lda_centers = lda_fit.transform(dummycenters)
fig, ax = plt.subplots()
ax.scatter(lda_data[:,0], lda_data[:,1], c=original_labels)
ax.plot(lda_centers[:,0], lda_centers[:,1],
"X", ms=15, alpha=0.9, c='black')
plt.show()
```
%% Output
%% Cell type:markdown id:17a2c48d-f0cd-4cac-9885-37b80b1d78a4 tags:
```{admonition} Bemerkung
Wenn man auf die Daten eine Standardskalierung anwendet, erhält man das gleiche Ergebnis - lediglich mit skalierten Achsen. Daher ist eine Skalierung bei Anwendung der LDA sinnlos.
```
Wir wollen nun noch sehen, wie gut LDA zur Klassifikation taugt. Dazu teilen wir die Daten in Trainings- und Testdaten auf und lassen die Testdaten von der auf den Trainingsdaten gefitteten LDA klassifizieren.
%% Cell type:code id:438f8032-ef67-46a4-9619-6439ece8c22b tags:
``` python
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(dummydata, original_labels, test_size=0.33)
lda_fit = LDA.fit(X_train, y_train)
new_test_labels = lda_fit.predict(X_test)
wrongly_classified = X_test[new_test_labels != y_test]
correct_label_of_wrongs = y_test[new_test_labels != y_test]
fig, ax = plt.subplots(figsize=(6,6))
ax.scatter(X_test[:,0], X_test[:,1],
c=new_test_labels,
alpha=0.6)
ax.scatter(wrongly_classified[:,0], wrongly_classified[:,1], marker='x',
c=correct_label_of_wrongs,
alpha=0.6)
plt.show()
```
%% Output
%% Cell type:markdown id:92a5656d-d60f-4a2a-a826-ee5017e3560d tags:
Zuletzt schauen wir uns wieder die Wein-Daten an:
%% Cell type:code id:f57985d6-3f83-43db-8222-503f1d0caa9a tags:
``` python
from sklearn.datasets import load_wine
a, target = load_wine(return_X_y=True)
fig, ax = plt.subplots()
ax.scatter(a[:,0], a[:,1], c=target)
plt.show()
```
%% Output
%% Cell type:code id:55112d22-413f-470c-91bf-b6574a91b16e tags:
``` python
lda_fit = LDA.fit(a, target)
lda_data = lda_fit.transform(a)
fig, ax = plt.subplots()
ax.scatter(lda_data[:,0], lda_data[:,1], c=target)
plt.show()
```
%% Output
%% Cell type:markdown id:b747286a-ce23-4976-ae7c-d663d9f838f7 tags:
Und das sieht tatächlich so aus, als könnte man die drei Klassen mit zwei Geraden voneinander separieren.
%% Cell type:markdown id:9c89d7af-60cb-4b74-82e1-06f31a90dbc7 tags:
### Überblick Scikit-Learn
## Überblick Scikit-Learn
[vgl. Scikit-learn User Guide](https://scikit-learn.org/stable/user_guide.html)
Wir haben nun schon viele Bestandteile von Scikit-Learn kennen gelernt:
1\. Supervised learning: discriminant_analysis.LinearDiscriminantAnalysis
2\. Unsupervised learning: clustering.KMeans
3\. Model selection and evaluation: model_selection.train_test_split
6\. Dataset transformations: preprocessing.StandardScaler
6\. Dataset transformations: decomposition.PCA
7\. Dataset loading utilities: datasets
Wir haben uns nicht mit den Teilen von Scikit-Learn beschäftigt, die im User Guide unter Inspection, Visualizations, Computing, Persistence aufgeführt sind.
Ein Aspekt, den wir auch nicht untersucht haben, der aber bei komplexeren Projekten mit Scikit-Learn hilfreich wird, ist das Basteln von Pipelines.
......
......@@ -1066,7 +1066,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
......@@ -1080,7 +1080,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.2"
"version": "3.8.10"
}
},
"nbformat": 4,
......
......@@ -641,7 +641,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
......@@ -655,7 +655,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.2"
"version": "3.8.10"
}
},
"nbformat": 4,
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment