Die lineare Regressionsanalyse ist ein häufiger Einstieg ins maschinelle Lernen um stetige Werte vorherzusagen (Prediction bzw. Prädiktion). Hinter der Regression steht oftmals die Methode der kleinsten Fehlerquadrate und die hat mehr als eine mathematische Methode zur Lösungsfindung (Gradientenverfahren und Normalengleichung). Alternativ kann auch die Maximum Likelihood-Methode zur Regression verwendet werden. Wir wollen uns in diesem Artikel nicht auf die Mathematik konzentrieren, sondern uns direkt an die Anwendung mit Python Scikit-Learn machen:
Haupt-Lernziele:
- Einführung in Machine Learning mit Scikit-Learn
- Lineare Regression mit Scikit-Learn
Neben-Lernziele:
- Datenvorbereitung (Data Preparation) mit Pandas und Scikit-Learn
- Datenvisualisierung mit der Matplotlib direkt und indirekt (über Pandas)
Was wir inhaltlich tun:
Der Versuch einer Vorhersage eines Fahrzeugpreises auf Basis einer quantitativ-messbaren Eigenschaft eines Fahrzeuges.
Die Daten als Download
Für dieses Beispiel verwende ich die Datei “Automobil_data.txt” von Kaggle.com. Die Daten lassen sich über folgenden Link downloaden, nur leider wird ein (kostenloser) Account benötigt:
https://www.kaggle.com/toramky/automobile-dataset/downloads/automobile-dataset.zip
Sollte der Download-Link unerwartet mal nicht mehr funktionieren, freue ich mich über einen Hinweis als Kommentar 🙂
Die Entwicklungsumgebung
Ich verwende hier die Python-Distribution Anaconda 3 und als Entwicklungs-Umgebung Spyder (in Anaconda enthalten). Genauso gut funktionieren jedoch auch Jupyter Notebook, Eclipse mit PyDev oder direkt die IPython QT-Console.
Zuerst einmal müssen wir die Daten in unsere Python-Session laden und werden einige Transformationen durchführen müssen. Wir starten zunächst mit dem Importieren von drei Bibliotheken NumPy und Pandas, deren Bedeutung ich nicht weiter erläutern werde, somit voraussetze.
import matplotlib.pyplot as plt # Die Nr.1 der Bibliotheken zur Datenvisualisierung
import numpy as np # Bibliothek "Nummerisches Python"
import pandas as pd # Bibliothek "Panel Data"
Wir nutzen die Pandas-Bibliothek, um die “Automobile_data.txt” in ein pd.DataFrame zu laden.
dataSet = pd.read_csv("Automobile_data.txt", # Hier liegt die Datei im selben Verzeichnis wie das Python-Skript!
delimiter = ',',
thousands = None,
decimal = '.')
Schauen wir uns dann die ersten fünf Zeilen in IPython via dataSet.head().
In : dataSet.head()
Out:
symboling normalized-losses make fuel-type aspiration num-of-doors \
0 3 ? alfa-romero gas std two
1 3 ? alfa-romero gas std two
2 1 ? alfa-romero gas std two
3 2 164 audi gas std four
4 2 164 audi gas std four
body-style drive-wheels engine-location wheel-base ... engine-size \
0 convertible rwd front 88.6 ... 130
1 convertible rwd front 88.6 ... 130
2 hatchback rwd front 94.5 ... 152
3 sedan fwd front 99.8 ... 109
4 sedan 4wd front 99.4 ... 136
fuel-system bore stroke compression-ratio horsepower peak-rpm city-mpg \
0 mpfi 3.47 2.68 9.0 111 5000 21
1 mpfi 3.47 2.68 9.0 111 5000 21
2 mpfi 2.68 3.47 9.0 154 5000 19
3 mpfi 3.19 3.4 10.0 102 5500 24
4 mpfi 3.19 3.4 8.0 115 5500 18
highway-mpg price
0 27 13495
1 27 16500
2 26 16500
3 30 13950
4 22 17450
[5 rows x 26 columns]
Hinweis: Der Datensatz hat viele Spalten, so dass diese in der Darstellung mit einem Backslash \ umgebrochen werden.
Gleich noch eine weitere Ausgabe dataSet.info(), die uns etwas über die Beschaffenheit der importierten Daten verrät:
In : dataSet.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 205 entries, 0 to 204
Data columns (total 26 columns):
symboling 205 non-null int64
normalized-losses 205 non-null object
make 205 non-null object
fuel-type 205 non-null object
aspiration 205 non-null object
num-of-doors 205 non-null object
body-style 205 non-null object
drive-wheels 205 non-null object
engine-location 205 non-null object
wheel-base 205 non-null float64
length 205 non-null float64
width 205 non-null float64
height 205 non-null float64
curb-weight 205 non-null int64
engine-type 205 non-null object
num-of-cylinders 205 non-null object
engine-size 205 non-null int64
fuel-system 205 non-null object
bore 205 non-null object
stroke 205 non-null object
compression-ratio 205 non-null float64
horsepower 205 non-null object
peak-rpm 205 non-null object
city-mpg 205 non-null int64
highway-mpg 205 non-null int64
price 205 non-null object
dtypes: float64(5), int64(5), object(16)
memory usage: 41.7+ KB
Einige Spalten entsprechen hinsichtlich des Datentypes nicht der Erwartung. Für die Spalten ‘horsepower’ und ‘peak-rpm’ würde ich eine Ganzzahl (Integer) erwarten, für ‘price’ hingegen eine Fließkommazahl (Float), allerdings sind die drei Spalten als Object deklariert. Mit Trick 17 im Data Science, der Anzeige der Minimum- und Maximum-Werte einer zu untersuchenden Datenreihe, kommen wir dem Übeltäter schnell auf die Schliche:
dataSet['horsepower'].min()
Out: '100'
dataSet['horsepower'].max()
Out: '?'
Datenbereinigung
Für eine Regressionsanalyse benötigen wir nummerische Werte (intervall- oder ratioskaliert), diese möchten wir auch durch richtige Datentypen-Deklaration herstellen. Nun wird eine Konvertierung in den gewünschten Datentyp jedoch an den (mit ‘?’ aufgefüllten) Datenlücken scheitern.
Schauen wir uns doch einmal die Datenreihen an, in denen in der Spalte ‘peak-rpm’ Fragezeichen stehen:
dataSet[dataSet['peak-rpm'] == '?'][['engine-type', 'num-of-cylinders']]
Out:
engine-type num-of-cylinders
130 ohc four
131 ohc four
Zwei Datenreihen sind vorhanden, bei denen ‘peak-rpm’ mit einem ‘?’ aufgefüllt wurde. Nun könnten wir diese Datenreihen einfach rauslöschen. Oder mit sinnvollen (im Sinne von wahrscheinlichen) Werten auffüllen. Vermutlichen haben beide Einträge – beide sind OHC-Motoren mit 4 Zylindern – eine ähnliche Drehzahl-Angabe wie vergleichbare Motoren. Mit folgendem Quellcode, gruppieren wir die Spalten ‘engine-type’ und ‘num-of-cylinders’ und bilden für diese Klassen den arithmetischen Mittelwert (.mean()) für die ‘peak-rpm’.
dataSet_rpm = dataSet[dataSet['peak-rpm'] != '?'][['engine-type', 'num-of-cylinders','peak-rpm']]
dataSet_rpm['peak-rpm'] = dataSet_rpm['peak-rpm'].astype(float)
dataSet_rpm_grouped = dataSet_rpm.groupby(['engine-type', 'num-of-cylinders'])
dataSet_rpm_grouped['peak-rpm'].mean()
Und schauen wir uns das Ergebnis an:
dataSet_rpm_grouped['peak-rpm'].mean()
Out:
engine-type num-of-cylinders
dohc four 5700.000000 --
six 5050.000000
dohcv eight 5750.000000
l four 4668.181818
three 5100.000000
ohc five 5081.818182
four 5155.468750
six 4821.428571
ohcf four 4775.000000
six 5900.000000
ohcv eight 4625.000000
six 5212.500000
twelve 5000.000000
rotor two 6000.000000
Name: peak-rpm, dtype: float64
Ein Vier-Zylinder-OHC-Motor hat demnach durchschnittlich einen Drehzahl-Peak von 5155 Umdrehungen pro Minute. Ohne nun (fahrlässigerweise) auf die Verteilung in dieser Klasse zu achten, nehmen wir einfach diesen Schätzwert, um die zwei fehlende Datenpunkte zu ersetzen.
Wir möchten jedoch die Original-Daten erhalten und legen ein neues DataSet (dataSet_c) an, in welches wir die Korrekturen vornehmen:
dataSet_c = dataSet.copy() # das "c"-Anhängsel steht für "corrected"
Nun können wir die fehlenden Peak-RPM-Einträge mit unserem Schätzwert ersetzen:
dataSet_c.loc[dataSet_c['peak-rpm'] == '?', 'peak-rpm'] = 5155
Was bei einer Drehzahl-Angabe noch funktionieren mag, ist für anderen Spalten bereits etwas schwieriger: Die beiden Spalten ‘price’ und ‘horsepower’ sind ebenfalls vom Typ Object, da sie ‘?’ enthalten. Verzichten wir einfach auf die betroffenen Zeilen:
dataSet_c = dataSet_c[dataSet_c['price'] != '?'] # entsprechende Zeilen herausfiltern
dataSet_c['price'] = dataSet_c['price'].astype(float) # Typ-Konvertierung zu Float
dataSet_c = dataSet_c[dataSet_c.horsepower != '?'] # entsprechende Zeilen herausfiltern
dataSet_c['horsepower'] = dataSet_c['horsepower'].astype(float) # Typ-Konvertierung in Int
Datenvisualisierung mit Pandas
Wir wollen uns nicht lange vom eigentlichen Ziel ablenken, dennoch nutzen wir die Visualisierungsfähigkeiten der Pandas-Library (welche die Matplotlib inkludiert), um uns dann die Anzahlen an Einträgen nach Hersteller der Fahrzeuge (Spalte ‘make’) anzeigen zu lassen:
dataSet_grouped_make = dataSet_c.groupby('make')
dataSet_grouped_make['make'].count().plot(kind = 'bar', figsize = (10, 10))
plt.show() # Besser jedes Plot abschließen! Auch wenn es in Pandas entstanden ist.
Oder die durchschnittliche PS-Zahl nach Hersteller:
(dataSet_c.groupby('make'))['horsepower'].mean().plot(kind = 'barh',
title = 'Mean Horsepower',
figsize = (10, 10))
plt.show()
Vorbereitung der Regressionsanalyse
Nun kommen wir endlich zur Regressionsanalyse, die wir mit Scikit-Learn umsetzen möchten. Die Regressionsanalyse können wir nur mit intervall- oder ratioskalierten Datenspalten betreiben, daher beschränken wir uns auf diese. Die “price”-Spalte nehmen wir jedoch heraus und setzen sie als unsere Zielgröße fest.
""" ----- Vorbereitung für die Regressionsanalyse ----- """
cols_ratio = ['horsepower', 'wheel-base', 'length', 'width', 'height', 'curb-weight', 'engine-size', 'compression-ratio', 'city-mpg', 'highway-mpg']
cols_target = ['price']
dataSet_ratio = dataSet_c.loc[:, cols_ratio]
dataSet_target = dataSet_c[cols_target]
Interessant ist zudem die Betrachtung vorab, wie die einzelnen nummerischen Attribute untereinander korrelieren. Dafür nehmen wir auch die ‘price’-Spalte wieder in die Betrachtung hinein und hinterlegen auch eine Farbskala mit dem Preis (höhere Preise, hellere Farben).
grr = pd.plotting.scatter_matrix(dataSet_c[cols_target + cols_ratio]
,c = dataSet_target
,figsize=(15, 15)
,marker = 'o'
,hist_kwds={'bins' : 20}
,s = 60
,alpha = 0.8)
plt.show()
Die lineare Korrelation ist hier sehr interessant, da wir auch nur eine lineare Regression beabsichtigen.
Wie man in dieser Scatter-Matrix recht gut erkennen kann, scheinen einige Größen-Paare nahezu perfekt zu korrelieren, andere nicht.
Korrelation…
- …nahezu perfekt linear: highway-mpg vs city-mpg (mpg = Miles per Gallon)
- … eher nicht gegeben: highway-mpg vs height
- … nicht linear, dafür aber nicht-linear: highway-mpg vs price
Nun, wir wollen den Preis eines Fahrzeuges vorhersagen, wenn wir eine andere quantitative Größe gegeben haben. Auf den Preis bezogen, erscheint mir die Motorleistung (Horsepower) einigermaßen linear zu korrelieren. Versuchen wir hier die lineare Regression und setzen somit die Spalte ‘horsepower’ als X und ‘price’ als y fest.
X = dataSet_ratio[['horsepower']] # doppelte [], da eine Liste von Spalten zu übergeben ist
y = dataSet_c[cols_target]
Die gängige Konvention ist übrigens, X groß zu schreiben, weil hier auch mehrere x-Dimensionen enthalten sein dürfen (multivariate Regression). y hingegen, ist stets nur eine Zielgröße (eine Dimension).
Die lineare Regression ist ein überwachtes Verfahren des maschinellen Lernens, somit müssen wir unsere Prädiktionsergebnisse mit Test-Daten testen, die nicht für das Training verwendet werden dürfen. Scitkit-Learn (oder kurz: sklearn) bietet hierfür eine Funktion an, die uns das Aufteilen der Daten abnimmt:
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y,
test_size = 0.3, # 70% der Daten für das Training
random_state = None) # bei Bedarf kann hier "dem Zufall auf die Sprünge geholfen" werden
Zu beachten ist dabei, dass die Daten vor dem Aufteilen in Trainings- und Testdaten gut zu durchmischen sind. Auch dies übernimmt die train_test_split-Funktion für uns, nur sollte man im Hinterkopf behalten, dass die Ergebnisse (auf Grund der Zufallsauswahl) nach jedem Durchlauf immer wieder etwas anders aussehen.
Lineare Regression mit Scikit-Learn
Nun kommen wir zur Durchführung der linearen Regression mit Scitkit-Learn, die sich in drei Zeilen trainieren lässt:
""" ----- Lineare Regressionsanalyse ------- """
from sklearn.linear_model import LinearRegression # importieren der Klasse
lr = LinearRegression() # instanziieren der Klasse
lr.fit(X_train, y_train) # trainieren
Aber Vorsicht! Bevor wir eine Prädiktion durchführen, wollen wir festlegen, wie wir die Güte der Prädiktion bewerten wollen. Die gängigsten Messungen für eine lineare Regression sind der MSE und R².
Ein großer MSE ist schlecht, ein kleiner gut.
Ein kleines R² ist schlecht, ein großes R² gut. Ein R² = 1.0 wäre theoretisch perfekt (da der Fehler = 0.00 wäre), jedoch in der Praxis unmöglich, da dieser nur bei absolut perfekter Korrelation auftreten würde. Die Klasse LinearRegression hat eine R²-Messmethode implementiert (score(x, y)).
print('------ Lineare Regression -----')
print('Funktion via sklearn: y = %.3f * x + %.3f' % (lr.coef_[0], lr.intercept_))
print("Alpha: {}".format(lr.intercept_))
print("Beta: {}".format(lr.coef_[0]))
print("Training Set R² Score: {:.2f}".format(lr.score(X_train, y_train)))
print("Test Set R² Score: {:.2f}".format(lr.score(X_test, y_test)))
print("\n")
Die Ausgabe (ein Beispiel!):
------ Lineare Regression -----
Funktion via sklearn: y = 170.919 * x + -4254.701 # Die Funktion ist als y = 171 * x - 4254.7
Alpha: [-4254.70114803] # y-Achsenschnitt bei x = 0
Beta: [ 170.91919086] # Steigung der Gerade
Training Set R² Score: 0.62
Test Set R² Score: 0.73
Nach jedem Durchlauf ändert sich mit der Datenaufteilung (train_test_split()) das Modell etwas und auch R² schwankt um eine gewisse Bandbreite. Berauschend sind die Ergebnisse dabei nicht, und wenn wir uns die Regressionsgerade einmal ansehen, wird auch klar, warum:
plt.figure(figsize=(10,10))
plt.scatter(X_train, y_train, color = 'blue') # Blaue Punkte sind Trainingsdaten
plt.scatter(X_test, y_test, color = 'green') # Grüne Punkte sind Testdaten
plt.plot(X_train, lr.predict(X_train), color = 'red') # Hier ensteht die Gerade (x, y) = (x, lr.predict(x)
plt.xlabel(X_train.columns[0])
plt.ylabel(cols_target[0])
plt.show()
Bei kleineren Leistungsbereichen, etwa bis 100 PS, ist die Preis-Varianz noch annehmbar gering, doch bei höheren Leistungsbereichen ist die Spannweite deutlich größer. (Nachträgliche Anmerkung vom 06.05.2018: relativ betrachtet, bleibt der Fehler über alle Wertebereiche ungefähr gleich [relativer Fehler]. Die absoluten Fehlerwerte haben jedoch bei größeren x-Werten so eine Varianz der möglichen y-Werte, dass keine befriedigenden Prädiktionen zu erwarten sind.)
Egal wie wir eine Gerade in diese Punktwolke legen, wir werden keine befriedigende Fehlergröße erhalten.
Nehmen wir einmal eine andere Spalte für X, bei der wir vor allem eine nicht-lineare Korrelation erkannt haben: “highway-mpg”
X = dataSet_ratio[['highway-mpg']]
y = dataSet_c[cols_target]
Wenn wir dann das Training wiederholen:
------ Lineare Regression -----
Funktion via sklearn: y = -868.787 * x + 40575.036
Alpha: [ 40575.03556055]
Beta: [-868.7869183]
Training Set R² Score: 0.49
Test Set R² Score: 0.40
Die R²-Werte sind nicht gerade berauschend, und das erklärt sich auch leicht, wenn wir die Trainings- und Testdaten sowie die gelernte Funktionsgerade visualisieren:
Die Gerade lässt sich nicht wirklich gut durch diese Punktwolke legen, da letztere eher eine Kurve als eine Gerade bildet. Im Grunde könnte eine Gerade noch einigermaßen gut in den Bereich von 22 bis 43 mpg passen und vermutlich annehmbare Ergebnisse liefern. Die Wertebereiche darunter und darüber jedoch verzerren zu sehr und sorgen zudem dafür, dass die Gerade auch innerhalb des mittleren Bereiches zu weit nach oben verschoben ist (ggf. könnte hier eine Ridge-/Lasso-Regression helfen).
Richtig gute Vorhersagen über nicht-lineare Verhältnisse können jedoch nur mit einer nicht-linearen Regression erreicht werden.
Nicht-lineare Regression mit Scikit-Learn
Nicht-lineare Regressionsanalysen erlauben es uns, nicht-lineare korrelierende Werte-Paare als Funktion zu erlernen. Im folgenden Scatter-Plot sehen wir zum einen die gewohnte lineare Regressionsgerade (y = a * x + b) in rot, eine polinominale Regressionskurve dritten Grades (y = a * x³ + b * x² + c * x + d) in violet sowie einen Entscheidungsweg einer Entscheidungsbaum-Regression in gelb.
Nicht-lineare Regressionsanalysen passen sich dem Verlauf der Punktwolke sehr viel besser an und können somit in der Regel auch sehr gute Vorhersageergebnisse liefern. Ich ziehe hier nun jedoch einen Gedankenstrich, liefere aber den Quellcode für die lineare Regression als auch für die beiden nicht-linearen Regressionen mit:
Python Script Regression via Scikit-Learn
Weitere Anmerkungen
- Bibliotheken wie Scitkit-Learn erlauben es, machinelle Lernverfahren schnell und unkompliziert anwenden zu können. Allerdings sollte man auch verstehen, wei diese Verfahren im Hintergrund mathematisch arbeiten. Diese Bibliotheken befreien uns also nicht gänzlich von der grauen Theorie.
- Statt der “reinen” lineare Regression (LinearRegression()) können auch eine Ridge-Regression (Ridge()), Lasso-Regression (Lasso()) oder eine Kombination aus beiden als sogenannte ElasticNet-Regression (ElasticNet()). Bei diesen kann über Parametern gesteuert werden, wie stark Ausreißer in den Daten berücksichtigt werden sollen.
- Vor einer Regression sollten die Werte skaliert werden, idealerweise durch Standardisierung der Werte (sklearn.preprocessing.StandardScaler()) oder durch Normierung (sklearn.preprocessing.Normalizer()).
- Wir haben hier nur zwei-dimensional betrachtet. In der Praxis ist das jedoch selten ausreichend, auch der Fahrzeug-Preis ist weder von der Motor-Leistung, noch von dem Kraftstoffverbrauch alleine abhängig – Es nehmen viele Größen auf den Preis Einfluss, somit benötigen wir multivariate Regressionsanalysen.
(Random Initialization). Wir berechnen die Ausgabe (Forwardpropogation) und vergleichen sie über eine Verlustfunktion (z. B. über die Funktion Mean Squared Error) mit dem tatsächlich korrekten Wert. Auf Grund der zufälligen Initialisierung haben wir eine nahe zu garantierte Falschheit der Ergebnisse und somit einen Verlust. Für die Verlustfunktion berechnen wir den Gradienten für gegebene Eingabewerte. Voraussetzung dafür ist, dass die Funktion ableitbar ist. Wir bewegen uns entgegen des Gradienten in Richtung Minimum der Verlustfunktion. Ist dieses Minimum (fast) gefunden, spricht man auch davon, dass der Lernalgorithmus konvergiert.
) der stets mit einer Eingangskonstante multipliziert und somit als Wert erhalten bleibt. Der Bias ist das Alpha 
in einer Schulmathe-tauglichen Formel wie 
.
ist die Steigung, der Gradient, der Funktion.
und den darauffolgenden Abgleich mit dem tatsächlichen 
herauszufinden. Anfangs behaupten wir beispielsweise einfach, sowohl 
. Folglich wird 
bzw. partiell nach jedem vorhandenen 
:
am Ende her haben? Das ergibt sie aus der Kettenregel: Die äußere Funktion wurde abgeleitet, so wurde aus 
dann 
. Jedoch muss im Sinne eben dieser Kettenregel auch die innere Funktion abgeleitet werden. Da wir nach 
erhalten.
ist der Gradient der Funktion!)
(eta) bezeichnet:
enthalten sind. Der Fehler ist dann aber sehr groß (sollte maximal sein, im Vergleich zu zukünftigen Epochen). Die Gewichte werden also angepasst, die Gerade somit besser in die Punktwolke platziert. Mit jeder Epoche wird die Gerade erneut in die Punktwolke gelegt, der Gesamtfehler (über alle 
konvergiert:
