PyData IT-Landschaft: Ein Software-Rundgang

Lässt man den Blick über die moderne Softwarelandschaft streifen, so erspäht der aufmerksame Beobachter, dass dort neben den zahlreichen Daten-Silos und Hadoop Data Lakes mehr und mehr Data Labs entstehen. Hier arbeiten Dateningenieure und Datenwissenschaftler Hand in Hand daran, aus den rohen und meist unstrukturierten Datenmassen wertvolle Datenprodukte zu formen. Besonders viele dieser Data Labs gehören dabei zur PyData Gruppe, einer offenen und gut vernetzten Community von Open Source Freunden, die eine Vielzahl an Bibliotheken und Tools einsetzt, um der universellen Programmiersprache Python an die Herausforderung von Datenanalyse und Machine Learning anzupassen. Auf diese Weise lässt sich die einsteigerfreundliche Syntax und weite Verbreitung von Python nutzen, ohne die eher langsame Ausführungsgeschwindigkeit der interpretierten Sprache in Kauf nehmen zu müssen. Zudem hat der offene Austausch und die rasante Entwicklung dazu geführt, dass das PyData Toolset numerischen Sprachen wie R oder Matlab in kaum noch etwas nachsteht und insbesondere im Machine- und Deep Learning Bereich zum verbreiteten Standard geworden ist.

NumPy - Das numerische Fundament

In diesem Beitrag wollen wir die PyData-Softwaregebilde genauer inspizieren und dabei starten wir mit dem Fundament, welches Python fit für komplexe mathematische Operationen und effiziente numerische Berechnung macht: NumPy. Inkludiert man NumPy's direktes Vorgängerprojekt Numeric in die Zeitrechnung, so ist dieses Python Fundament schon über 25 Jahre alt und entsprechend stabil sowie weit verbreitet. Neben einer zugänglichen Schreibweise für Matrixoperationen basiert NumPy's Erfolg insbesondere auf seiner Performanz.

Als interpretierte und dynamisch typisierte Sprache ist Python für Entwickler ausgesprochen komfortabel. Gleichzeitig erzeugt dieser Komfort aber zusätzlichen Aufwand beim Interpreter, etwa für Speichermanagement und Type Checking. Offensichtlich wird dies bei der mehrfachen Ausführung einer trivialen Funktion: Möchte man etwa eine Millionen Fließkommazahlen addieren, so sähe eine reine Python-Implementierung etwa wie folgt aus.

In [ ]:
import random
In [ ]:
%%time
python_floats = [random.uniform(0,1) for _ in range(1000000)]
In [ ]:
%%time
sum(python_floats)

Zunächst wird hier mit Hilfe einer List Comprehension ein Listenobjekt mit zufälligen Werten erstellt, anschließend wird die native Python Implementierung zur Summenbildung genutzt. Alternativ dazu können auch entsprechende NumPy Funktionen genutzt werden, wodurch die Ausführungszeit deutlich sinkt:

In [ ]:
import numpy as np
In [ ]:
%%time
numpy_floats = np.random.uniform(low=0, high=1., size=(1000000))
In [ ]:
%%time
np.sum(numpy_floats)

NumPy erreicht diese Geschwindigkeitsvorteile, indem es die Zahlen in einer eigenen Datenstruktur ablegt, dem sogenannten ndarray. Innerhalb dieses haben alle Elemente einen einheitlichen Datentyp (dtype) sowie eine definierte Größe (shape), welche die maximale Ausprägung einer jeden Dimension ausgibt. In folgenden Beispiel handelt es sich um ein 1-dimensionales ndarray vom type float64, jedoch können ndarrays prinzipiell beliebige Dimensionalität haben.

In [ ]:
print(type(numpy_floats)) # Element-Typ
print(numpy_floats.dtype) # Typ 
print(numpy_floats.shape) # (1000000,)

Basierend auf dieser Datenstruktur bietet NumPy nun diverse Operationen, welche mittels Cython implementiert sind, einer speziellen Python Erweiterung, die C-ähnliche Geschwindigkeit ermöglicht.

Auch wenn NumPy als Grundlage für den PyData Stack unumgänglich ist, so ähnelt die Nutzung tatsächlich der eines Kellers: Manch einer verbringt Tage und Wochen hier, um an Projekten zu tüfteln, andere gehen quasi nie hinab, sondern verbringen ihre Zeit lieber etwas weiter oben, wo es meist optisch ansprechender gestaltet ist und weitere Bequemlichkeiten zur Verfügung stehen. Dennoch hilft es immer, sich auch im Keller gut auszukennen, denn der Grundriss des gesamten Hauses wird schließlich hier bestimmt, weshalb wir uns nun noch einige weitere wichtige NumPy Funktionalitäten anschauen.

Numpy bietet diverse praktische Funktionen, um Vektorräume beliebiger Dimensionalität aufzuspannen. Mit Funktionen wie zeros oder ones lassen sich schnell und einfach ndarrays mit einem einheitlichen Wert und einfach konfigurierbarer Dimension erzeugen. Diese besitzen zudem praktische Hilfsfunktionen, etwa um alle Elemente der verschachtelten Arraystruktur in eine einfache, flache Liste zu transformieren.

In [ ]:
np.linspace(0,10,5) # [0.,2.5,5.,7.5,10.]
In [ ]:
zeros = np.zeros((2, 3, 4)) # Erzeugung eines 3-Dimensionalen Vektorraums mit insgesamt 24 Elementen (2*3*4)
zeros.shape # Gibt die Dimensionalität aus: (2,3,4)
len(zeros.flatten()) # Anzahl der Einträge: 24

Darüber hinaus bietet NumPy zahlreiche weitere Funktionen um auch komplexe Datensequenzen zu generieren. Besonders häufig wird np.arange genutzt, um Zahlenketten mit einem festen Abstand zu erzeugen. Soll nicht der Abstand sondern die Anzahl der Elemente fixiert und gleichmäßig verteilt werden, so bietet sich np.linespace an. Zudem bietet np.tile die Möglichkeit gesamte Arrays zu wiederholen, während np.repeat die Wiederholung auf jedes einzelne Element anwendet.

In [ ]:
np.arange(0,3) # array([0, 1, 2, 3])
In [ ]:
np.linspace(1,10, num=5) #array([ 1.  ,  3.25,  5.5 ,  7.75, 10.  ])
In [ ]:
np.repeat(np.arange(0,3),2)
In [ ]:
np.tile(np.arange(0,3),2)

Ein weitere typischer Anwendungsfall ist die Generierung von zufälligen Werten, wofür wir bereits zum Beginn des Artikels das Modul np.random genutzt haben. Dieses bietet darüber hinaus über 30 unterschiedliche Funktionen zur Generierung statistischer Verteilungen wie etwa Binär-, Normal- oder Geometrischer Verteilung. Die so erzeugten Werte können anschließend mit Hilfe der reshape Funktion in die gewünschte mehrdimensionale Struktur überführt werden.

In [ ]:
np.random.seed=42
three_dimensions = np.random.uniform(size=30).reshape(5,3,2) # Erzeugt eine Dreifach verschachtelte Struktur, die 5 Blöcke erzeugt, welche jeweils die Struktur 3*2 haben.

Hat man das gewünschte Array erzeugt, kann es direkt mittels Standardoperatoren oder NumPy-Funktionen verändert werden. Von einfacher Summenbildung, über Matrix-Multiplikation, elementweisen Vergleichsoperationen bis hin zu trigonometrischen und Exponentialfunktionen bietet NumPy alle üblichen mathematischen Operationen. Die Syntax ist dabei nicht nur intuitiv, sondern die Ausführung auch meist wesentlich schneller, als eine äquivalente Funktion aus den Python Standardbibliotheken.

Möchte man etwa jeden Wert einer Liste in Python um 1 erhöhen, so wäre hier erneut eine Schleife notwendig, etwa in Form einer List-Comprehension. Numpy hingegen erlaubt es, den skalaren Wert 1 einfach zu addieren, wandelt diesen dann intern in eine entsprechende Vektordarstellung um und addiert diese hinzu - mit einem Geschwindigkeitsvorteil Faktor 100.

In [ ]:
%timeit numpy_floats + 1
In [ ]:
%timeit [f+1 for f in python_floats]

Neben den Operatoren der Grundrechenarten sind Vergleichsoperatoren wichtig. Diese geben einen Boolean Wert zurück, welcher dann insbesondere zur Datenselektion genutzt werden kann. Hierfür wird eine ähnliche Syntax wie für Python Listen genutzt:

In [ ]:
numpy_floats[1:4] # Gibt das zweite, dritte und vierte Element zurück
numpy_floats[numpy_floats > 0.9] # Gibt alle Werte größer 0.9 zurück
numpy_floats[(0.99 > numpy_floats) & (numpy_floats < 0.01)] # Gibt alle Werte kleiner als 0.01 oder größer als 0.99 zurück

Ein häufiger Anwendungsfall mehrdimensionaler Zahlenstrukturen ist die Bestimmung von Extremwerten oder anderer statistischer Kennzahlen. Hier werden die Stärken NumPys nicht nur bei der Geschwindigkeit deutlich (Faktor 30), sondern es zeigt sich auch die Mächtigkeit der Syntax. Während Aggregationen standardmäßig auf allen Elementen ausgeführt werden, können sie auch mit Hilfe des Parameters axis nur für einzelne Subgruppen gebildet werden:

In [ ]:
%timeit max(python_floats)
In [ ]:
%timeit numpy_floats.max()
In [ ]:
three_dimensions.max() # Das Maximum aller Werte
three_dimensions.max(axis=0) # 6 Maximalwerte
three_dimensions.max(axis=1) # 10 Maximalwerte
three_dimensions.max(axis=2) # 15 Maximalwerte
three_dimensions.max(axis=2).max(axis=1) # 5 Maximalwerte (das Maximum für jeden äußeren Block)

Natürlich lassen sich diese Filter-, Vergleichs- und Aggregationsfunktionen auch miteinander kombinieren. Möchten wir etwa von den 5 Blöcken nur solche nutzen, welche einen Durchschnittswert von über 0.7 haben, so könnten wir das Arithmetrische Mittel über die Achsen 1 und 2 bilden, um diese mittels Vergleichsoperator in True/False Werte zu konvertieren und anschließend zur Filterung zu nutzen:

In [ ]:
three_dimensions[three_dimensions.mean(axis=2).mean(axis=1) < 0.4]

Pandas - Die eigenen n Wände

Wie bereits erwähnt, bildet NumPy zwar das essenzielle Fundament, viele Nutzer-/innen haben mit diesem jedoch kaum syntaktische Berührungspunkte. Vielmehr ist die auf NumPy aufbauende Bibliothek Pandas weit verbreitet und bietet eine erweiterte Syntax, in der sich viele zu Hause fühlen.

Pandas wurde von Wes McKinney zunächst aufgrund seines beruflichen Eigenbedarfs an einer performanten und flexiblen Bibliothek zur Arbeit mit tabellarischen Daten und Zeitreihen entwickelt und im Jahr 2008 open-sourced. Neben Import- und Exportfunktionalitäten bietet Pandas vor allem erweiterte Transformationen wie etwa Slicing, Merges oder Pivotisierungen. Darüber hinaus verfügt Pandas über erweiterte Indizierungsmöglichkeiten zur schnellen Filterung und Gruppierung von Daten und vereinfacht den Umgang mit fehlenden Werten.

In Pandas existieren zwei elementare Datenstrukturen:

  • Eine Serie (pd.Series) ist ein schlanker Wrapper um ein 1-dimensionales NumPy Array. Entsprechend besitzt eine Serie stets nur einen dtype, stellt also eine Gruppierung homogener Daten des selben Types dar.
  • Ein DataFrame (pd.DataFrame) ist im Kern eine Ansammlung vieler Serien, welche ähnlich wie bei einer Tabelle als Spalten zusammengefasst werden. In einem DataFrame können also auch unterschiedliche Datentypen erfasst werden.

Pandas unterstützt beim Import eine Vielzahl an Formaten wie Parquet, Avro, Feather oder auch Excel. Das meist genutzte Format ist aber vermutlich noch immer der CSV-Standard. Dabei kann die Datei entweder lokal gespeichert sein oder direkt von einer URL geladen werden:

In [ ]:
import pandas as pd
df = pd.read_csv('https://s3.amazonaws.com/tripdata/201801-citibike-tripdata.csv.zip', compression='zip')
df.columns = df.columns.str.replace(" ","_")

Beim Einlesen werden die Daten geparsed, automatisch interpretiert und jede Spalte nach dem Prinzip der kleinsten allgemeingültigen Gemeinsamkeit einem Datentypen (dtype) zugewiesen. Beinhaltet eine Spalte etwa ausschließlich ganzzahlige Werte, wird der int64-Datentyp abgeleitet, ist eine Kommazahl darunter, wird die gesamte Serie dem Typ float64 zugeordnet. Alle Spalten stehen anschließend gemeinsam in Form eines DataFrames im Arbeitsspeicher zur Verfügung und sind nach Zeilennummer indiziert. Mit der info-Methode können die interpretierten Datentypen kontrolliert und einige andere erste Informationen wie die Anzahl der NULL-Werte, die Zahl der Zeilen und Spalten sowie der Speicherbedarf abgelesen werden.

In [ ]:
df.info()

DataFrames und Serien nutzen zwecks Datenhaltung die ndarray-Struktur aus NumPy und erben damit auch die Überlagerung zahlreicher Standardoperatoren, so dass etwa Vergleiche zeilenweise durchgeführt werden. Dazu erweitert Pandas ndarrays um zusätzliche, intuitive Zugriffsmöglichkeiten und Hilfsfunktionalitäten. Soll etwa eine Spalte eines DataFrames (welche einer Serie entspricht) ausgewählt werden, so ist dies über die von Python Dictionaries bekannte Klammer-Notation wie auch die von Datenklassen bekannte Punkt-Notation möglich.

In [ ]:
df.tripduration == df["tripduration"] # Vergleicht die tripduration-Spalte mit sich selbst und gibt für jede Zeile entsprechend True zurück

Pandas bietet darüber hinaus unterschiedliche Wege um Daten (auch) zeilenweise zu filtern:

  • .iloc ermöglicht den Datenzugriff anhand der Position innerhalb der Serie bzw. des DataFrames
  • .loc liest die selektierten Zeilen anhand ihres Index Wertes. Darüber hinaus akzeptiert loc aber auch Serien-Objekte der gleichen Länge mit Boolean Werten.
  • .query bietet eine SQL-ähnliche Syntax zur Selektion von Daten

Die unterschiedliche Nutzung verdeutlichen die folgenden Beispiel-Abfragen:

In [ ]:
df.iloc[:3, [0,3]] # Die ersten drei Zeilen und dazu die 3. und 4. Spalte (Indizierung startet bei 0)
df.loc[(df.usertype == "Subscriber") & (df.start_station_latitude>40.76), ["starttime","stoptime"]] # Die Start- und Stopzeit von Abonnenten deren Tour östlich des Breitengrads 40.76 begann
df.query("tripduration > 3600 and birth_year>2000") # Alle Fahrten die länger als eine Stunde gingen und von Kindern nach Geburtsjahr 2000 unternommen wurden
df.loc[:100].query("start_station_id == end_station_id") # Die Kreisfahrten (gleiche End- wie Startstation) innerhalb der ersten 100 Einträge

Während die loc und iloc Funktionalität sowohl für Serien als auch für DataFrames funktionieren, ist die query Notation nur auf Tabellen verfügbar. Ähnlich verhält es sich mit zahlreichen weiteren praktischen Funktionen, welche teils für beide Datentypen, teils nur für einen zur Verfügung stehen. Bei DataFrames kann meist mit Hilfe des axis Attributes ähnlich wie bei NumPy entschieden werden, ob eine Funktion Zeilen- oder Spaltenweise durchgeführt wird.

In [ ]:
# Funktionalitäten für DataFrames und Serien
df.count() # Zählt die nicht leeren Einträge jeder Spalte
df.count(axis=1) # Zählt die nicht leeren Einträge jeder Zeile
df.end_station_id.nunique() # Gibt die Anzahl der in den Daten vorkommenden Endstationen

# Spannende Serien Funktionalitäten
df.birth_year.value_counts() # Zählt die Anzahl der Einträge nach Geburtsjahr
df.groupby("usertype").birth_year.mean() # Gibt das durchschnittliche Geburtsjahr der unterschiedlichen Nutzertypen aus
df.groupby("usertype").birth_year.describe() # Gibt statistische Werte wie Durchschnitt, Extremwerte etc. je Nutzertypen aus

Pandas bietet darüber hinaus hilfreiche Funkionen, etwa für den Umgang mit Text, Datumsangaben und Zeitreihen. Hierfür folgend einige Beispiele:

In [ ]:
df["end_station_name"].str.contains("Broadway") # Gibt für jede Zeile zurück, ob der Stationsname des Endpunkts Broadway beinhaltet
pd.to_datetime(df.starttime).dt.weekday # Der Wochentag an dem eine Fahrt statt gefunden hat
df.groupby(pd.to_datetime(df.starttime).dt.dayofyear).size().rolling(7).mean() # Auf 7 Tage geglättete Anzahl an Fahrten je Tag

Zudem können Informationen nicht nur als Zahlenwerte zurückgegeben werden, sondern auch direkt in eine visuelle Darstellung überführt werden:

In [ ]:
df.groupby(pd.to_datetime(df.starttime).dt.dayofyear).birth_year.agg(["mean","median"]).plot.line()
In [ ]:
df.groupby(pd.to_datetime(df.starttime).dt.day).birth_year.mean().rolling(7, center=True).agg([np.mean, lambda x: x.iloc[3]]).rename(columns={"<lambda>": "raw"}).plot.line()
In [ ]:
df.groupby(pd.to_datetime(df.starttime).dt.day).birth_year.mean().rolling(7).agg([np.mean, lambda x: x.iloc[0]]).rename(columns={"<lambda>": "raw"}).plot.line()

Scikit-Learn - Die Data Science Grundversorgung

Die Analyse eines Datensatzes kann zwar an sich bereits einen großen Wert bieten, jedoch lassen sich häufig durch den Einsatz von Machine Learning (ML) viele weitere, wertvolle Erkenntnisse generieren. Ähnlich wie ein stabiles Haus erst mit fließend Wasser und Strom wirklich bewohnbar wird, so erzeugt der Einsatz von Regressions-, Klassifikations- und Clusteringverfahren aus einem vorverarbeiteten, gesäuberten Datensatz einen echten wirtschaftlichen Nutzen. Generell bietet Python zahlreiche ML-Frameworks und insbesondere Deep Learning Frameworks wie Tensorflow und PyTorch haben entschieden zur Popularität der PyData Plattform beigetragen. Dennoch ist Deep Learning nur ein spezielles Teilgebiet des Maschinellen Lernens, welches zwar mächtig aber auch sehr daten- und resourcenhungrig ist, so dass klassische Machine Learning Verfahren auf tabellarischen Daten häufig schneller und einfacher wertvolle Ergebnisse erzielen.

Die wichtigsten und gängigsten Algorithmen stehen dabei im scikit-learn Framework zur Verfügung und können über ein einfaches, wohl strukturiertes Interface angesprochen werden. Darüber hinaus hält scikit-learn eine große Anzahl weiterer notwendiger Funktionalitäten wie One-Hot-Encoding, Datenskalierung, Umgang mit Fehlwerten und Train-Test-Split Erzeugung bereit. Diese Transformatoren folgen in scikit-learn genauso wie die eigentlichen ML-Alogirthmen dem fit-transform bzw. fit-predict Muster. Der fit Schritt entspricht dem Training, bei dem die Parameter so angepasst werden, dass die Ausgaben möglichst gut dem Zeilwerten (Label) entsprechen. Der transform (bei Vorverarbeitungsschritten) bzw. predict (bei den ML-Algorithmen) hingegen bezeichnet die Anwendung des gefitteten Transformators bzw. Modells auf neueen Daten, um möglichst passende Ausgaben zu generieren.

Ein klassisches Minimalbeispiel für die Nutzung von scikit-Learn sieht wie folgt aus:

In [ ]:
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import accuracy_score, f1_score

target_col = "usertype"
X_train, X_test, y_train, y_test = train_test_split(df.drop(columns=[target_col]).select_dtypes(np.number), df[target_col], test_size=0.2)

gbc = GradientBoostingClassifier(subsample=0.1, n_estimators=50)
gbc.fit(X_train, y_train)
y_pred = gbc.predict(X_test)
accuracy_score(y_true=y_test, y_pred=y_pred)

Als erstes unterteilen wir unseren Datensatz in 80% Trainings- und 20% Testdaten, wobei diese jeweils in Feature-Matrix (X) und Zielvektor (y) getrennt werden. In diesem Schritt löschen wir zudem auch alle nicht numerischen Spalten, da die meisten Modelle hiermit nicht umgehen können. Anschließend erstellen wir ein GradientBoosting-Modellobjekt, ein beliebter und häufig erfolgreicher Klassifikator, der zudem mit fehlenden und nicht skalierten Werten umgehen kann. Zwecks Demonstration reduzieren wir die Ausführungszeit durch die Beschränkung auf nur 10% der Trainingsdaten und maximal 50 Ebenen. Anschließend trainieren wir das Modell, lassen eine Prognose berechnen und bestimmen die Genauigkeit (Accuracy) des Modells, also den Anteil der richtigen Vorhersagen. Ein Wert von 98,4% klingt zwar schon sehr gut, jedoch wird hier ein typisches Problem des Accuarcy Scores offensichtlich: Bei starkem Klassen-Ungleichgewicht sind sehr hohe Werte einfach zu erreichen. Zur Bewertung sollte daher immer eine sogenannte Baseline als Referenz herangezogen werden. Im Fall einer binären Klassifikation eignet sich hierfür einfach stets den häufigsten Wert zu prognostizieren:

In [ ]:
accuracy_score(y_test, np.repeat("Subscriber", len(y_test)))

Unser Modell schlägt also unsere Baseline, allerdings wirkt die Zahl von 98,4% nun nicht mehr so überwältigend, wenn wir wissen, dass wir auch in 96,9% der Fälle richtig liegen, wenn wir einfach immer nur behaupten, der Nutzer wäre ein Subscriber.

PyViz - Pythons Gestaltungs- & Farbvielfalt

Ähnlich wie bei der Wohnungseinrichtung eine Abwägung zwischen funktionalen und ästhetischen Ansprüchen gefunden werden muss, so müssen auch Datenvisualisierungen einen Kompromiss aus Aussagemächtigkeit und Intuitivität finden, welche von relativ einfachen Diagrammen zur Ergebnisdarstellung, bis zu komplexen, interaktiven Visualisierungen zur Datenexploration reichen kann. Bei derart vielfältigen Anforderungen überrascht die breite Palette an unterschiedlichen Visualisierungsbibliotheken im Python Ökosystem nicht. Bei der Auswahl lassen sich objektive Kriterien wie etwa das Abstraktionslevel bzw. die Mächtigkeit und der Funktionsumfang genauso heranziehen, wie subjektive Präferenzen bezüglich Syntax und ästethischer Geschmack.

Eine umfassende Vorstellung der verfügbaren Bibliotheken würde diesen Artikel sprengen, viel mehr sei hierfür auf das PyViz Projekt verwiesen, welches das Ziel hat, Nutzern alle notwendigen Informationen bereitzustellen, um die für sie optimale Plotting-Bibliothek auszuwählen. Die Website listet dabei nicht nur die gängigen Tools auf, sondern bietet auch eine kategorisierte Übersicht an Bibliotheken, die sich für Dashboarding eignen und solche, die die besonders einfache Pandas High-Level-API unterstützen.

Zur groben Orientierung, folgend eine Übersicht über die gängigsten Bibliotheken, die wir im kommenden Abschnitt nutzen werden:

  • MatplotLib ist wohl das NumPy der Datenvisualisierung: Langjährig bewährt und essentielle Grundlage für zahlreiche weitere Bibliotheken, gleichzeitig aber etwas umständlich für komplexere Aufgaben. MatplotLib wurde 2002 von John Hunter entwickelt und direkt open sourced, so dass es Schnell eine große Nutzerbass aufbauen konnte, welche ihr bis heute treu geblieben ist und insbesondere die zahlreichen detailierten Einstellungsmöglichkeiten schätzt. Der direkte Umgang mit Matplotlib wird zudem auch bei einigen darauf aufbauenden Bibliotheken noch häufig gebraucht, wenn das höhere Abstraktionsniveau an seine Grenzen kommt.

  • Seaborn ist der wohl bekannteste Kandidat, für eine auf MatplotLib aufbauende Bibliothekn mit höherem Abstraktionsniveau, welches vor allem komplexe Darstellungen einfach macht. Dabei bietet es sowohl komplexe Darstellungen wie Violin- oder Distributionplots wie auch eine große Anzahl an Parameter in jeder Funktion, womit die Graphiken einfach an die eigenen Wünsche angepasst werden können.

  • AltAir ist ein eine noch eher junges Projekt (2015) und folgt dem "Grammer of Graphics" Prinzip, welches auch schon als Vorlage für die in der R Welt beliebte ggplot2 Bibliothek dient. Dies ermöglicht eine sehr mächtige und konsistente Deklaritions-Syntax. Technisch basiert Altair anders als Seaborn auf dem Visualisierungsframework Vegas.

  • Plotly und Bokeh sind zwei weitere Framework, welche bei der Erzeugung interaktiver Darstellung brillieren. Während Bokeh ein klassisches open-source Projekt ist, wird Plotly wird aktiv von einer kleinen Firma vorangetrieben, welche die eigentliche Plotting-Bibliothek ebenfalls frei zur Verfügung stellt. und nutzt einen eigenen Renderer. Darüber hinaus bieten die Plotly Macher mit Plotly Express eine Variante mit sehr intuitiver Syntax sowie mit Dash ein Framework zur Erstellung von Dashboards basierend aus zahlreichen Plotly-Darstellungen.

Praktische Implementierung

Bei der Suche, nach einem perfekten Haus oder einer perfekten Wohnung, ist es wichtig alle Fakten und Zahlen zusammenzutragen und einige Fotos zu sehen, um sich ein grobes Bild machen zu können. Dennoch ist eine Wohnungsbesichtigung unerlässlich, denn erst wenn man selbst hindurch läuft kann man sich vorstellen, wie es sein könnte hier zu leben. Aus diesem Grund soll dieser Artikel nicht mit den theoretischen Beschreibungen der Frameworks enden, sondern ihr Zusammenwirken in einem praktischen Beispiel demonstrieren. Dafür wollen wir versuchen im bereits bekannten NYC Bikesharing Datensatz das Alter der Kunden zu prognostizieren.

Data Exploration & Data Cleaning

Jedes Machine Learning Projekt startet zunächst mit einer Exploration der Daten. Aus dem Pandas Absatz haben wir bereits eine grobe Vorstellung von den verfügbaren Spalten. Dieses mal soll sich unser Modell aber nicht ausschließlich auf nativ numerische Werte beschränken, daher lohnt sich ein erneuter Blick auf die Datentypen:

In [ ]:
df.dtypes

Die Start- und Stopzeit wird von Pandas nicht automatisch als Zeitpunkt erkannt, jedoch lässt sich dies einfach ändern. Zudem sind die Attribute bezüglich Start- und Endstation schwierig. Diesen widmen wir uns später. Der Nutzertyp ist uns bereits besser bekannt, hier existieren zwei Ausprägungen, welche wir einfach in 1 und 0 transformieren können:

In [ ]:
df["starttime"] = pd.to_datetime(df.starttime)
df["stoptime"] = pd.to_datetime(df.stoptime)

df["subscriber"] = (df.usertype == "Subscriber").astype(int)

Des Weiteren fällt auf, dass das Geschlecht nicht als Boolean, sondern als Integer-Typ gespeichert ist. Hier lohnt sich ein genauerer Blick, wofür ein Histogram bezüglich der Zielvariable ein gängiges Mittel ist:

In [ ]:
import plotly.express as px
px.histogram(df.sample(frac=.9), color="gender", x="birth_year")

Die häufige Kombination des Geschlechts 0 mit dem Geburtsjahr 1969 ist äußerst unplausibel und legt nahe, dass dies eventuell die Default-Werte im Registrierungsprozess darstellen könnten. In jedem Fall dürften diese Angaben schlicht nicht stimmen. Es gibt nun verschiedene Wege mit diesen falschen Werten umzugehen, die einfachste ist jedoch, die entsprechenden Einträge einfach zu entfernen:

In [ ]:
df = df.query("not (gender == 0 and birth_year==1969)")

Ein besonderes Augenmerk verdient zudem stets die Zielvariable:

In [ ]:
df.birth_year.describe()

Bei deren Betrachtung fällt auf, dass einerseits manche Geburtsjahre sehr unplausibel scheinen, denn das Geburtsjahr 1885 entspräche zum Ausleihzeitpunkt einem Alter von 132 Jahren. Des Weiteren ist der Proxy Geburtsjahr für die Zielvariable Alter etwas umständlich, so dass sich eine direkte Transformation lohnt. In diesem Zuge ist es wichtig, die Ursprungsvariable Geburtsjahr nicht während es Modelltrainings zu nutzen, daher löschen wir sie am besten direkt:

In [ ]:
df["age"] = df.starttime.dt.year - df.birth_year
df.drop(columns="birth_year", inplace=True)
df = df.query("age < 100")

Neben dem Historgram eignet sich auch ein sogenannter (Gaussian) Kernel Density Estimation Plot (KDE), welcher eine möglichst passende Dichtefunktion schätzt und darstellt:

In [ ]:
df.age.plot.kde(bw_method=0.1, xlim=(0, 100))

Diese Darstellung zeigt recht deutlich die Schieflage der Altersverteilung: In den Kinderjahren ist die Nutzerzahl nahe Null, steigt dann rapide auf den Höchstwert um etwa 30 Jahre an und fällt dann langsam nach hinten ab. Diese Schieflage (skewness) wird auch deutlich, wenn man Mittelwert und Arithmetisches Mittel vergleicht:

In [ ]:
print("Häufigster Wert", df.age.mode()[0])
print("Mittelwert", df.age.median())
print("Durchschnitt", df.age.mean())

Ein erstes Modell

Für den Anfang sollen uns diese Informationen für ein erstes Modell genügen. Data Science ist bekanntlich ein iterativer Prozess, bei dem konstant zwischen Analyse, Modelltraining und Evaluation gewechselt wird. Der Quellcode für das Modell zur Altersprognose ähnelt stark dem zur Usertype-Klassifikation, jedoch handelt es sich mit dem Alter nun um eine kontinuierliche Zielvariable, entsprechend müssen wir ein Regressionsverfahren wählen und die Bewertungsfunktion anpassen:

In [ ]:
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_absolute_error

target_col = "age"
X_train, X_test, y_train, y_test = train_test_split(
    df.drop(columns=[target_col]).select_dtypes(np.number),
    df[target_col],
    test_size=0.2,
    random_state=0,
)

gbr = GradientBoostingRegressor(subsample=0.1, n_estimators=50, random_state=0)
gbr.fit(X_train, y_train)
y_pred = gbr.predict(X_test)
mean_absolute_error(y_true=y_test, y_pred=y_pred)

Dieses mal erhalten wir als Bewertungsmetrik 9,71, was die durchschnittliche Abweichung unserer Prognose vom echten Alter des Nutzers bezeichnet. Auch hier lohnt sich das Heranziehen einer Baseline, wofür sich gerade bei einer schiefen Verteilungsfunktion der Median besser eignet als der Durchschnitt:

In [ ]:
train_median = y_train.median()
mean_absolute_error(y_test, np.repeat(train_median, len(y_test)))

Die Prognose des mittleren Alters erzielt einen absoluten Fehler von 9,84, wir sind also mit unserem Modell etwas besser, der Unterschied von 0,13 bedeutet aber, dass wir uns im Mittel nur bei jeder achten Person um ein Jahr weniger stark verschätzen.

Nun stehen uns diverse nächste Schritte zur Auswahl, um unsere Prognose zu verbessern: Wir können uns auf die Suche nach besseren Hyperparametern für unser Modell machen oder andere Modelle ausprobieren. Wir könnten uns auch erneut der Datengrundlage widmen, Daten bereinigen oder weitere Jahre hinzunehmen. Beides wird unser Modell vermutlich verbessern, jedoch ist das Hyperparameter-Tuning inzwischen ein weitgehend automatisierter, rechenaufwändiger Prozess und die Datenbereinigung meist Fleißarbeit. Wir wollen im Folgenden daher lieber versuchen weitere Features zu generieren, welche unserem Modell wertvolle Informationen darreicht.

Dafür müssen wir zunächst verstehen, welche Informationen unser Modell derzeit als wichtig erachtet. Bei Baum-Algoirthmen ist die feature_importance, welche jedem Eingangsfeature ein bestimmtes Gewicht zurechnet, der erste Anlaufpunkt:

In [ ]:
pd.Series(gbr.feature_importances_, X_train.columns).sort_values().plot.bar()

Es zeigt sich deutlich, dass der Ort scheinbar wichtig ist, denn das Modell rechnet den Koordinaten das höchste Gewicht zu, obwohl sie aufgrund ihrer unabhängigen Darstellung der Achsen nicht optimal aufbereitet sind. Auch fällt auf, dass diverse wichtige Informationen, wie etwa die Uhrzeit, derzeit noch garnicht genutzt werden. Wir berechnen daher aus unseren Timestamp-Spalten den Wochentag sowie die Stunde zu der das Fahrrad ausgeliehen bzw. zurückgegeben wurde:

In [ ]:
df["weekday"] = df.starttime.dt.weekday

df["hour_start"] = df.starttime.dt.hour
df["hour_end"] = df.stoptime.dt.hour

Mit diesen einfachen Features verbessern wir uns bereits auf 9,64. Als nächstes widmen wir uns der Kombination aus Ort und Zeit und berechnen mit Hilfe des geopy-Pakets für jede Fahrt anhand der Start- und Zielkoordinaten die zurückgelegte Strecke. Da diese Berechnung relativ aufwändig ist und nicht vektorisiert zur Verfügung steht bestimmen wir zunächst die Anzahl aller existierenden Streckenkombinationen und berechnen jeweils die Distanz. Anschließend joinen wir diese Distanz-Information anhand der IDs der Start- und Zielstation. Auf diese Weise sparen wir uns nicht nur die Distanzberechnung für mehrfach existierende ID-Kombinationen, sondern können die Distanz-Informationen etwa auch für zukünftige Läufe cachen und skalieren besser, falls wir später weitere Fahrten hinzunehmen:

In [ ]:
import geopy.distance

get_distance = lambda r: geopy.distance.distance(
    [r.start_station_latitude, r.start_station_longitude],
    [r.end_station_latitude, r.end_station_longitude],
).m

distances = df[
    [
        "start_station_id",
        "end_station_id",
        "start_station_latitude",
        "start_station_longitude",
        "end_station_latitude",
        "end_station_longitude",
    ]
].drop_duplicates()
distances["distance"] = distances.apply(get_distance, axis=1)
distances = (
    distances[["start_station_id", "end_station_id", "distance"]]
    .drop_duplicates(subset=["start_station_id", "end_station_id"])
    .set_index(["start_station_id", "end_station_id"])
)

df = df.join(distances, on=["start_station_id", "end_station_id"], how="left")

Mit den Informationen über Dauer und Strecke einer Fahrt können wir nun auch sehr einfach eine Geschwindigkeit in km/h berechnen:

In [ ]:
df["speed"] = df.distance/1000/(df.tripduration/3600)
df.speed.hist(bins=100)

Die Verteilung der Geschwindigkeit entspricht einer Normalverteilung mit Erwartungswert knapp unter 10. Die niedrige Durchschnittsgeschwindigkeit ist zum einen vermutlich dem New Yorker Verkehr geschuldet, des Weiteren berechnen wir nur die Luftliniendistanz zwischen Start- und Zielstation und ignorieren jegliche Umwege. Besonders deutlich wird dies bei den Fahrten mit einer Geschwindigkeit von Null, welche fast ausnahmslos darauf zurückzuführen ist, dass das Fahrrad an der gleichen Station abgeholt und zurückgegeben wurde. Entsprechend sollten diese Fahrten mit einem zusätzlichen Feature gekennzeichnet werden:

In [ ]:
df["roundtrip"] = (df.start_station_id == df.end_station_id).astype(int)
In [ ]:
px.histogram(df.sample(frac=.1), color="roundtrip", x="speed")

Führen wir erneut unser Training aus, reduziert sich unser durchschnittlicher Fehler mit den neuen Features auf 9.55 und die Geschwindigkeit wird vom Modell nun als wichtigstes Feature bewertet:

In [ ]:
X_train, X_test, y_train, y_test = train_test_split(
    df.drop(columns=[target_col]).select_dtypes(np.number),
    df[target_col],
    test_size=0.2,
    random_state=0,
)

gbr.fit(X_train, y_train)
y_pred = gbr.predict(X_test)
mean_absolute_error(y_true=y_test, y_pred=y_pred)
In [ ]:
pd.Series(gbr.feature_importances_, X_train.columns).sort_values().plot.bar()

Trotz der Bedeutung der Geschwindigkeit bewertet unser Modell die Koordinaten noch immer als sehr wichtig. Entsprechend sollten wir versuchen eine bessere Repräsentation zu entwickeln. Eine visuelle Darstellung ist hierfür meist hilfreich: Wir summieren die Fahrten nach ihren Startpunkten auf und berechnen für jeden das mittlere Alter sowie die Anzahl der von hier getätigten Fahrten. Das Ergebnis lässt sich am besten mit Hilfe von Plotly-Express durch eine Punktwolke vor einer entsprechenden Landkarte visualisieren.

In [ ]:
plot_df = (
    df.query("start_station_latitude <41")
    .groupby(["start_station_latitude", "start_station_longitude"])
    .age.agg(["median", "size"])
    .reset_index()
    .dropna()
)

plot_df = (
    df.query("start_station_latitude <41")
    .groupby(["start_station_latitude", "start_station_longitude"])
    .age.agg(["median", "size"])
    .reset_index()
    .dropna()
)

try:
    px.set_mapbox_access_token(open(".mapbox_token").read())
    fig = px.scatter_mapbox(
        plot_df,
        lat="start_station_latitude",
        lon="start_station_longitude",
        color="median",
        size="size",
        zoom=10,
        width=10
    )
    fig.show()
except:
    print("this chart requires a mapbox token")
    

Es zeigt sich, dass das mittlere Alter zwischen den Gebieten deutlich variiert: Rund um die Universität ist das mittlere Alter etwa deutlich jünger als im Finanzdistrikt. Diese Zusammenhänge sind zwar mittels Clustering detektierbar und könnten dann in Form eines Model-Esembles integriert werden, wir nutzen aber einen einfacheren Ansatz und ignorieren den tatsächlichen örtlichen Zusammenhang der Stationen. Vielmehr nutzen wir die Stations-IDs und berechnen für jede Station den Altersmedian und joinen ihn anschließend an unsere Fahrten. Wichtig ist hierbei kein "information leakage" zu erzeugen, also strikt zwischen Trainings- und Testdaten zu trennen, wobei uns das Konstrukt eines Scikit-Learn Transformers unterstützt. Dieses Interface folgt auch dem fit-transform Design: Im fit-Schritt berechnen wir auf den Trainingsdaten für jede Station einen Altersmittelwert, im transform-Schritt werden diese Durchschnittswerte an die Test-Daten angefügt. Taucht eine Station nur im Testdatensatz auf und es liegt kein spezifischer Mittelwert vor, fügen wir einen globalen default-Wert ein:

In [ ]:
from sklearn.base import TransformerMixin, BaseEstimator

class StatsTransformer(BaseEstimator, TransformerMixin):
    def __init__(self, column=None, func=np.median, default_val=None):
        self.column = column
        self.func = func
        self.groups = None
        self.default_val = None
        super(StatsTransformer, self).__init__()

    def fit(self, X, y=None, **kwargs):
        data = X.copy()
        data["target"] = y
        self.groups = (
            data.groupby(self.column)
            .target.apply(self.func)
            .rename(f"{y.name}_{self.func.__name__}_on_{self.column}")
        )
        if self.default_val is None:
            self.default_val = self.func(data.target)
        return self

    def transform(self, X, y=None, **kwargs):
        X = (
            X.join(self.groups, on=self.column)
            .fillna(self.default_val)
            .drop(columns=self.column)
        )
        return X

    def fit_transform(self, X, y=None, **kwargs):
        ab = self.fit(X, y)
        return ab.transform(X, y)

Diesen Transformer müssen wir genauso wie ein Modell zuerst instanziieren und anschließend fitten. Um dies nicht für jeden Schritt einzeln ausführen zu müssen, bietet Scitkit-Learn Pipelines, mit welchen verschiedene aufeinander aufbauende Schritte gemeinsam ausgeführt werden können.

In [ ]:
from sklearn.pipeline import make_pipeline

pipeline = make_pipeline(
    StatsTransformer("start_station_id"), 
    StatsTransformer("end_station_id"), 
    gbr
)

pipeline.fit(X_train, y_train)

y_pred = pipeline.predict(X_test)
mean_absolute_error(y_true=y_test, y_pred=y_pred)

Mit den neuen Features schafft es unser Modell nun auf 9,33. Ein Blick auf die Feature-Importance zeigt deutlich die Bedeutung der neuen Features und wie entsprechend den reinen Koordinaten nur noch wenig Gewicht beigemessen wird:

In [ ]:
pd.Series(
    pipeline.steps[2][1].feature_importances_,
    list(X_train.columns.drop(["start_station_id", "end_station_id"]))
    + ["age_median_on_start", "age_median_on_end"],
).sort_values().plot.bar()

Mit den erstellten Features konnten wir die Qualität unseres Modells stetig verbessern, jedoch ist der Fehler des Modells im Vergleich zur Baseline nur um 0,5 verbessert. Natürlich besteht die Möglichkeit uns durch anderen Modelle, bessere Hyperparameter oder zusätzliche Features weiter zu verbessern, jedoch wird bereits offensichtlich, dass dies nicht automatisch passiert. Vielmehr ist strukturierte Arbeit erforderlich, bei der die unterschiedlichen Tools zur Datentransformation, Visualisierung und zum Modelltraining eng ineinander greifen. Die möglichen nächsten Schritte sind dabei vielfältig und ihr Erfolg im Vorhinein schwer abzuschätzen. Umso wichtiger ist ein strukturiertes Vorgehen und ein pflichtbewusstest geführtes Protokol der durchgeführten Experimente, denn nicht ohne Grund heißt dieses Aufgabengebiet Data Science.