Benjamin Aunkofer

Lineare Regression in Python mit Scitkit-Learn

/in , , , , , /by

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².

MSE = \frac{\sum_{i=1}^n (y_i - \hat{y_i})^2}{n}

Ein großer MSE ist schlecht, ein kleiner gut.

R^2 = 1 - \frac{MSE}{Var(y)}= \frac{\frac{1}{n} \cdot \sum_{i=1}^n (y_i - \hat{y_i})^2}{\frac{1}{n} \cdot \sum_{i=1}^n (y_i - \hat{\mu_y})^2}

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.
Benjamin Aunkofer

Entropie – Und andere Maße für Unreinheit in Daten

/in , , , , , , /by

Dieser Artikel ist Teil 1 von 4 der Artikelserie Maschinelles Lernen mit Entscheidungsbaumverfahren.

Hierarchische Klassifikationsmodelle, zu denen das Entscheidungsbaumverfahren (Decision Tree) zählt, zerlegen eine Datenmenge iterativ oder rekursiv mit dem Ziel, die Zielwerte (Klassen) im Rahmen des Lernens (Trainingsphase des überwachten Lernens) möglichst gut zu bereiningen, also eindeutige Klassenzuordnungen für bestimmte Eigenschaften in den Features zu erhalten. Die Zerlegung der Daten erfolgt über einen Informationsgewinn, der für die Klassifikation mit einem Maß der Unreinheit berechnet wird (im nächsten Artikel der Serie werden wir die Entropie berechnen!) Read more

Christoph Gresch

Der Blick für das Wesentliche: Die Merkmalsselektion

/in , , , , , , , , , /by

In vielen Wissensbasen werden Datensätze durch sehr große Merkmalsräume beschrieben. Während der Generierung einer Wissensbasis wird versucht jedes mögliche Merkmal zu erfassen, um einen Datensatz möglichst genau zu beschreiben. Dabei muss aber nicht jedes Merkmal einen nachhaltigen Wert für das Predictive Modelling darstellen. Ein Klassifikator arbeitet mit reduziertem Merkmalsraum nicht nur schneller, sondern in der Regel auch weitaus effizienter. Oftmals erweist sich ein automatischer Ansatz der Merkmalsselektion besser, als ein manueller, da durchaus Zusammenhänge existieren können, die wir selbst so nicht identifizieren können.

Die Theorie: Merkmalsselektion

Automatische Merkmalsselektionsverfahren unterscheiden 3 verschiedene Arten: Filter, Wrapper und Embedded Methods. Einen guten Überblick über Filter- und Wrapper-Verfahren bieten Kumari et al. in ihrer Arbeit “Filter versus wrapper feature subset selection in large dimensionality micro array: A review” (Download als PDF).

Der Filter-Ansatz bewertet die Merkmale unabhängig des Klassifikators. Dabei werden univariate und multivariate Methoden unterschieden. Univariate Methoden bewerten die Merkmale separat, während der multivariate Ansatz mehrere Merkmale kombiniert. Für jedes Merkmal bzw. jedes Merkmalspaar wird ein statistischer Wert berechnet, der die Eignung der Merkmale für die Klassifikation angibt. Mithilfe eines Schwellwertes werden dann geeignete Merkmale herausgefiltert. Der Filter-Ansatz bietet eine schnelle und, aufgrund der geringen Komplexität, leicht skalierbare Lösung für die Merkmalsselektion. Der Nachteil von Filter-Selektoren besteht in der Missachtung der Abhängigkeiten zwischen den Merkmalen. So werden redundante Merkmale ähnlich bewertet und verzerren später die Erfolgsrate des Klassifikators. Bekannte Beispiele für Filter-Selektoren sind unter anderem die Euklidische Distanz und der Chi-2-Test.

Der Wrapper-Ansatz verbindet die Merkmalsbewertung mit einem Klassifikator. Innerhalb des Merkmalsraumes werden verschiedene Teilmengen von Merkmalen generiert und mithilfe eines trainierten Klassifikators getestet. Um alle möglichen Teilmengen des Merkmalsraumes zu identifizieren, wird der Klassifikator mit einem Suchalgorithmus kombiniert. Da der Merkmalsraum mit Zunahme der Anzahl der Merkmale exponentiell steigt, werden heuristische Suchmethoden für die Suche nach optimalen Teilmengen genutzt. Im Gegensatz zu den Filtern können hier redundante Merkmale abgefangen werden. Die Nutzung eines Klassifikators zur Bewertung der Teilmengen ist zugleich Vor- und Nachteil. Da die generierte Teilmenge auf einen speziellen Klassifikator zugeschnitten wird, ist nicht gewährleistet, dass die Menge auch für andere Klassifikatoren optimal ist. Somit ist dieser Ansatz zumeist abhängig vom gewählten Klassifikator. Zudem benötigt der Wrapper-Ansatz eine viel höhere Rechenzeit. Wrapper-Selektoren werden beispielsweise durch Genetische Algorithmen und Sequentielle Forward/Backward-Selektoren vertreten.

Embedded-Ansätze stellen eine Sonderform der Wrapper-Methode da. Allerdings werden Merkmalssuche und Klassifikatoren-Training nicht getrennt. Die Suche der optimalen Teilmenge ist hier im Modelltraining eingebettet. Dadurch liefern Embedded-Ansätze die gleichen Vorteile wie die Wrapper-Methoden, während die Rechenzeit dabei erheblich gesenkt werden kann. Der reduzierte Merkmalsraum ist aber auch hier vom jeweiligen Klassifikator abhängig. Klassifikatoren, die den Embedded-Ansatz ermöglichen sind beispielsweise der Random-Forest oder die Support-Vector-Maschine.

Entwicklungsgrundlage

Analog zum letzten Tutorial wird hier Python(x,y) und die Datenbasis „Human Activity Recognition Using Smartphones“ genutzt. Die Datenbasis beruht auf erfassten Sensordaten eines Smartphones während speziellen menschlichen Aktivitäten: Laufen, Treppen hinaufsteigen, Treppen herabsteigen, Sitzen, Stehen und Liegen. Auf den Aufzeichnungen von Gyroskop und Accelerometer wurden mehrere Merkmale erhoben. Die Datenmenge, alle zugehörigen Daten und die Beschreibung der Daten sind frei verfügbar.

(https://archive.ics.uci.edu/ml/datasets/Human+Activity+Recognition+Using+Smartphones)

Alle Daten liegen im Textformat vor. Für ein effizienteres Arbeiten mit der Datenbasis wurden diese im Vorfeld in das csv-Dateiformat überführt.

Python-Bibliotheken

Alle für das Data Mining relevanten Bibliotheken sind in Python(x,y) bereits enthalten. Für die Umsetzung werden folgende Bibliotheken genutzt:

import numpy as np
import pandas as pd

from sklearn.cross_validation import StratifiedKFold
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import f_classif, RFECV, SelectKBest
from sklearn.svm import SVC

Die Bibliotheken NumPy und Pandas unterstützen die Arbeit mit verschiedenen Datenstrukturen und scikit-learn umfasst alle Funktionen des maschinellen Lernens.

Daten vorbereiten

Vor der Anwendung der einzelnen Verfahren werden die Daten vorbereitet. Das Data Frame wird eingelesen, die Klassen in numerische Labels überführt und das Datenfeld in Merkmale (X) und Klassenspalte (y) separiert. Weiterhin wird die informationslose Spalte subject entfernt.

index = 0
selected_features = []

# reading database
data = pd.read_csv("data/measures.csv", sep = ';', decimal = ',')

# converting textual class labels to numeric classes like description
data = data.replace({'WALKING': 1, 'WALKING_UPSTAIRS': 2, 'WALKING_DOWNSTAIRS': 3,
              'SITTING': 4, 'STANDING': 5, 'LAYING': 6})

# drop subject column
data = data.drop('subject', 1)

# remove class column from data set
print "removing class column from training set.."
X = data.drop('activity', 1)
y = data['activity']

columns = X.columns.values.tolist()

1. Verfahren: RFECV

Der RFECV (Recursive Feature Elimination with Cross Validation) ist ein Vertreter des Wrapper-Ansatzes. In diesem Beispiel wird die Merkmalsselektion mit einem Support Vector Klassifikator kombiniert. Der RFECV berechnet ein Ranking über die einzelnen Merkmale. Dabei bestimmt der Selektor selbst die optimale Menge der Merkmale. Alle Merkmale mit Platz 1 im Ranking bilden den optimalen Merkmalsraum.

''' ########## METHOD 1: RFE with cross validation and SVC ########## '''
print "create classifier for feature selection.."
svc = SVC(kernel = 'linear')

# fit the feature selector
print "create the feature selector.."
rfecv = RFECV(estimator = svc, step = 1, cv = StratifiedKFold(y, 3), scoring = 'accuracy')
print "fit the selector for data set.."
rfecv.fit(X, y)

print "The estimated number of optimal features is: " + str(rfecv.n_features_)

# get the most importent features
feat_importence = zip(rfecv.ranking_, columns)

# prepare list of selected features for new DataFrame
for i in range(len(feat_importence)):
    if(feat_importence[i][0] == 1):#>= np.nanmean(rfecv.ranking_)):
        selected_features.append(feat_importence[i][1])
        print "added feature: " + str(feat_importence[i][1]) + ".."

2. Verfahren: Random Forest-Klassifikator

Der Random-Forest-Klassifikator gehört zu den Modellen, die einen Embedded-Ansatz ermöglichen. Während des Klassifikatoren-Trainings wird jedem Merkmal ein Wert zugeordnet. Je höher der Wert, desto bedeutsamer das Merkmal. Allerdings ist hier eine manuelle Filterung notwendig, da anders als beim RFECV kein internes Optimum ermittelt wird. Mithilfe eines geeigneten Schwellwertes können die zu wählenden Merkmale bestimmt werden. In diesem Beispiel werden alle Merkmale selektiert, die eine Wichtung größer dem Mittelwert erhalten.

''' ########## METHOD 2: Random Forrest Classifier Feat Importance ########## '''
print "create classifier for feature selection.."
rfc = RandomForestClassifier(n_estimators = 500, criterion = 'entropy', max_depth = 4)
rfc = rfc.fit(X, y)

# get the most importent features
feat_importence = zip(rfc.feature_importances_, columns)

# prepare list of selected features for new DataFrame
for i in range(len(feat_importence)):
    if(feat_importence[i][0] >= np.mean(rfc.feature_importances_)):
        selected_features.append(feat_importence[i][1])
        print "added feature: " + str(feat_importence[i][1]) + ".."

3. Verfahren: Select K Best

Das Select K Best-Verfahren gehört den Filter-Ansätzen an. Daher kommt hier anders als bei den anderen beiden Verfahren kein Klassifikator zum Einsatz. Auch in diesem Verfahren wird für jedes Merkmal ein Wert berechnet, der die Wichtigkeit des Merkmals beziffert. Für die Berechnung der Werte können verschiedene Methoden verwendet werden. In diesem Beispiel wird eine Varianzanalyse genutzt (Parameter f_classif). Auch hier wird mithilfe eines manuellen Schwellwertes der reduzierte Merkmalsraum bestimmt.

''' ########## METHOD 3: Select K Best Features ########## '''
print "create classifier for feature selection.."
skb = SelectKBest(f_classif)
skb = skb.fit(X, y)

# get the most importent features
feat_importence = zip(skb.scores_, columns)

# prepare list of selected features for new DataFrame
for i in range(len(feat_importence)):
    if(feat_importence[i][0] >= np.nanmean(skb.scores_)):
        selected_features.append(feat_importence[i][1])
        print "added feature: " + str(feat_importence[i][1]) + ".."

Ergebnisse

Für die Bewertung der einzelnen Selektionsverfahren werden die einzelnen Verfahren in den Data-Mining-Prozess (siehe vorheriges Tutorial: Einstieg in das maschinelle Lernen mit Python(x,y)) integriert. Die nachfolgende Tabelle veranschaulicht die Ergebnisse der Klassifikation der einzelnen Verfahren.

 

Selektionsverfahren

Anzahl der Merkmale

Erfolgsrate Klassifikation
Ohne

561

93,96%
RFECV

314

94,03%
Random Forest

118

90,43%
Select K Best

186

92,30%

 

Durch den RFECV konnte das Ergebnis der Klassifikation leicht verbessert werden. Die anderen Selektionsverfahren, die auch deutlich weniger Merkmale nutzen, verschlechtern das Ergebnis sogar. Dies liegt vor allem an der manuellen Regulierung des Schwellwertes.

Christoph Gresch

Einstieg in das Maschinelle Lernen mit Python(x,y)

/in , , , , , , , , , , /by

Python(x,y) ist eine Python-Distribution, die speziell für wissenschaftliche Arbeiten entwickelt wurde. Es umfasst neben der Programmiersprache auch die Entwicklungsumgebung Spyder und eine Reihe integrierter Python-Bibliotheken. Mithilfe von Python(x,y) kann eine Vielzahl von Interessensbereichen bearbeitet werden. Dazu zählen unter anderem Bildverarbeitung oder auch das maschinelle Lernen. Das All-in-One-Setup für Python(x,y) ist für alle gängigen Betriebssysteme online erhältlich. Read more

Benjamin Aunkofer

Numerical Python – Einführung in wissenschaftliches Rechnen mit NumPy

/in , , , /by

NumPy steht für Numerical Python und ist eines der bekanntesten Pakete für alle Python-Programmierer mit wissenschaftlichen Hintergrund. Von persönlichen Kontakten erfuhr ich, dass NumPy heute in der Astrophysik fast genauso verwendet wird wie auch von sogenannten Quants im Investment-Banking. Das NumPy-Paket ist sicherlich ein Grundstein des Erfolges für Python in der Wissenschaft und für den häufigen Einsatz für die Implementierung von Algorihtmen des maschinellen Lernens in Python.

Die zentrale Datenstruktur in NumPy ist das mehrdimensionale Array. Dieses n-dimensionale Array (ndarray) ist eine sehr mächtige Datenstruktur und verwende ich beispielsweise in meinem Artikel über den k-Nächste-Nachbarn-Algorithmus. Die Besonderheit des NumPy-Arrays ist, dass es ein mehrdimensionaler Container für homogene Daten ist. Ein Datentyp gilt also für das gesamte Array, nicht nur für bestimmte Zeilen oder Spalten!

import numpy as np

Read more

Benjamin Aunkofer

Wahrscheinlichkeitsverteilungen – Zentralen Grenzwertsatz verstehen mit Pyhton

/in , , , , , , , /by

Wahrscheinlichkeitsverteilung sind im Data Science ein wichtiges Handwerkszeug. Während in der Mathevorlesung die Dynamik dieser Verteilungen nur durch wildes Tafelgekritzel schwierig erlebbar zu machen ist, können wir mit Programmierkenntnissen (in diesem Fall wieder mit Python) eine kleine Testumgebung für solche Verteilungen erstellen, um ein Gefühl dafür zu entwickeln, wie unterschiedlich diese auf verschiedene Wahrscheinlichkeitswerte, Varianz und Mengen an Datenpunkten reagieren und wann sie untereinander annäherungsweise ersetzbar sind – der zentrale Grenzwertsatz. Den Schwerpunkt lege ich in diesem Artikel auf die Binominal- und Normalverteilung.

Für die folgenden Beispiele werden folgende Python-Bibliotheken benötigt:

import matplotlib.pyplot as pyplot
import random as random
import math as math

Read more

Benjamin Aunkofer

ABC-XYZ-Analyse

/in , , , , , , , /by

Die ABC-XYZ-Analyse ist eine aussagekräftige Analyse für die Strategiefindung in der Warenwirtschaft und Logistik bzw. im Supply Chain Management. Die Analyse basiert auf der Vorstellung einer Pareto-Verteilung, die darauf hindeutet, dass oftmals eine kleine Menge eines großen Ganzen einen unverhältnismäßig großen Einfluss auf eben dieses große Ganze hat.

Die ABC-XYZ-Analyse beinhaltet im ersten Schritt eine ABC- und im zweiten Schritt eine XYZ-Analyse. Im dritten Schritt werden die Ergebnisse in einer Matrix zusammengeführt. In diesem Artikel erläutere ich nicht, wofür eine ABC-XYZ-Analyse dient und wie die Ergebnisse zu interpretieren sind, hier kann ich jedoch auf einen älteren Artikel “ABC-XYZ-Analyse” – www.der-wirtschaftsingenieur.de vom 3. Mai 2011 von mir verweisen, der vorher lesenswert ist, wenn kein Vorwissen zur ABC-XYZ-Analyse vorhanden ist.

Die Vorarbeit

Für die ABC- und XYZ-Analyse benötigen wir folgende Python-Bibliotheken:

import pandas as pd
import numpy as np
import random as random
import matplotlib.pyplot as pyplot

Wir laden die EKPO-Tabelle in ein DataFrame (Datenstruktur der Pandas-Bibliothek):

EKPO = pd.read_csv("[PFAD]EKPO.csv", delimiter=';', thousands='.', decimal=',')

Die Datei stammt aus einem SAP-Testsystem und steht hier zum Download bereit:

csv-icon

SAP.EKPO

Wir benötigen daraus nur folgende Zeilen:

EKPO_X = EKPO[['MATNR', 'MATKL', 'MENGE', 'PEINH', 'NETPR', 'NETWR']].copy()

Jetzt kommt der erste Kniff: Das Feld “MENGE” im SAP beschreibt die Menge in der jeweiligen Mengeneinheit (z. B. Stück, Meter oder Liter). Da wir hier jedoch nicht den genauen Verbrauch vorliegen haben, sondern nur die Einkaufsmenge (indirekt gemessener Verbrauch), sollten wir die Menge pro Preiseinheit “PEINH” berücksichtigen, denn nach dieser Preiseinheitsmenge erfolgt der Einkauf.

EKPO_X['Preiseinheitsmenge'] = EKPO_X['MENGE'] / EKPO_X['PEINH']

Für die Preiseinheitsmenge ein Beispiel:
Sie kaufen sicherlich pro Einkauf keine 3 Rollen Toilettenpapier, sondern eine oder mehrere Packungen Toilettenpapier. Wenn Sie zwei Packung Toilettenpapier für jeweils 2 Euro kaufen, die jeweils 10 Rollen beinhalten, ist die Preiseinheit = 10 und die Preiseinheitsmenge => 20 gekaufte Toilettenrollen / 10 Rollen pro Packung = 2 Packungen Toilettenpapier.

Nun haben wir also unsere für den Einkauf relevante Mengeneinheit. Jetzt sortieren wir diese Materialeinkäufe primär nach dem Umsatzvolumen “NETWR” absteigend (und sekundär nach der Preiseinheitsmenge aufsteigend, allerdings spielt das keine große Rolle):

EKPO_X = EKPO_X.sort_values(by = ['NETWR', 'Preiseinheitsmenge'], ascending=[False, True]) # Sortierung nach Umsatzvolumen pro Bestellung absteigend

Einige Störfaktoren müssen noch bereinigt werden. Erstens sollen Einträge mit Preisen oder Umsätzen in Höhe von 0,00 Euro nicht mehr auftauchen:

EKPO_X = EKPO_X[(EKPO_X.NETPR != 0) & (EKPO_X.NETWR != 0)]

Zweitens gibt es Einkäufe, die ein Material ohne Materialnummer und/oder ohne Materialklasse haben. Bei einer Zusammenfassung (Aggregation) über die Materialnummer oder die Materialklasse würden sich diese “leeren” Einträge als NULL-Eintrag bündeln. Das wollen wir vermeiden, indem wir alle NULL-Einträge mit jeweils unterschiedlichen Zufallszahlen auffüllen.

EKPO_X.MATNR[EKPO_X.MATNR.isnull() == True] = EKPO_X.MATNR[EKPO_X.MATNR.isnull() == True].apply(lambda x: random.random()) # Manche MATNR fehlen (NULL), diese füllen wir mit zufälligen Werten auf. Dabei ist es natürlich wichtig, dass die Zufallszahl für jede Zeile neu generiert wird! EKPO_X.MATNR.fillna(random.random()) funktioniert nicht, denn hier würde ein gleicher Wert alle NaN-Werte ersetzen

ABC – Analyse:

Nun geht es an die eigentliche ABC-Analyse, dafür müssen wir die Gruppierung der Materialien vornehmen. Gleich vorweg: Dies sollte man eigentlich über die einzelnen Materialnummern machen, da dies jedoch in der Visualisierung (auf Grund der hohen Anzahl und Vielfältigkeit) etwas aufwändiger ist, machen wir es über die Materialklassen. Wir gehen dabei einfach davon aus, dass die Materialklassen relativ homogene Materialien zusammenfassen und somit auch das Verbrauchs-/Einkaufverhalten innerhalb einer Gruppe nicht sonderlich viel Abweichung aufweist.

# Aggregation über die Materialklasse, Aufsummierung der Umsätze, Mengen und Volumen 
MATKL_MENGEN = (EKPO_X.MENGE.groupby(EKPO_X.MATKL).sum()).to_frame()
MATKL_PREISEINHEIT_MENGE = (EKPO_X.Preiseinheitsmenge.groupby(EKPO_X.MATKL).sum()).to_frame()
MATKL_VOLUMEN = (EKPO_X.NETWR.groupby(EKPO_X.MATKL).sum()).to_frame()

# Aggregation über die Materialklasse, Berechnung des Durchschnittpreises (ist bei einer Materialklasse, allerdings wenig sinnvoll!)
MATKL_Preise = (EKPO_X.NETPR.groupby(EKPO_X.MATKL).mean()).to_frame()EKPO_G = MATKL_MENGEN.join(MATKL_PREISEINHEIT_MENGE, how='left')

# Zusammenfügen der Ergebnisse (Left-Join)
EKPO_G = EKPO_G.join(MATKL_Preise, how='left')
EKPO_G = EKPO_G.join(MATKL_VOLUMEN, how='left')
EKPO_G = EKPO_G.sort_values(['NETWR'], ascending=False)

# Berechnung der kumulierten Umsätze und Mengen (Beachte: Vorher muss nach Umsätzen absteigend sortiert worden sein! (siehe oben)
EKPO_G['Volumen_kumuliert'] = EKPO_G.NETWR.cumsum()
EKPO_G['Menge_kumuliert'] = EKPO_G.MENGE.cumsum()

Nun können wir uns ganz im Sinne der ABC-Analyse die typische Pareto-Verteilung der kumulierten Umsätze (Umsatzgrößen absteigend sortiert) ansehen:

EKPO_G[['Menge_kumuliert','Volumen_kumuliert']].plot([EKPO_G.Menge_kumuliert, EKPO_G.Volumen_kumuliert], color=['red','pink'], figsize=[20,10], fontsize=8, title='Kumulierte Werte - Sortierung nach Materialklassen-Volumen')

abc_analyse_sap_netwr_menge_kumulierte_kurve_pareto

Die X-Achse zeigt die Materialklassen von links nach rechts in der Sortierung nach dem Umsatzvolumen (größester Umsatz links, kleinster Umsatz rechts). Die Y-Achse zeigt den Betrag der Umsatzhöhe (Euro) bzw. der Menge (Preiseinheitsmenge). Die Kurve der Menge ist mit Vorsicht zu bewerten, da primär nach dem Umsatz und nicht nach der Menge sortiert wurde.

Klassifikation:

Nun kommen wir zur Klassifikation. Hier machen wir es uns sehr einfach: Wir gehen einfach davon aus, dass 80% des Wertbeitrages aller Umsätze von etwa 20% der Materialien (hier: Materialklassen) umfassen und klassifizieren daher über feste relative Größen:

EKPO_G['ABC_Gruppe'] = "C" # Erstmal sind alle Materialien der C-Gruppe zugeordnet
EKPO_G['ABC_Gruppe'][EKPO_G.Volumen_kumuliert <= EKPO_G.NETWR.sum() / 100 * 95] = 'B' # Materialien, deren kumuliertes Volumen maximal 95% des Gesamtvolumens umfassen, sind Gruppe B
EKPO_G['ABC_Gruppe'][EKPO_G.Volumen_kumuliert <= EKPO_G.NETWR.sum() / 100 * 80] = 'A' # Materialien, deren kumuliertes Volumen maximal 80% des Gesamtvolumens umfassen, sind Gruppe A

Hinweis:
Intelligenter wird so eine Klassifikation, wenn wir den steilsten Anstieg innerhalb der kumulierten Volumen (die zuvor gezeigte Kurve) ermitteln und danach die Grenzen für die A-, B-, C-Klassen festlegen.

Optional: Farben für die Klassen festlegen (für die nachfolgende Visualisierung)

EKPO_G['Color'] = 'red'
EKPO_G['Color'][EKPO_G['ABC_Gruppe'] == 'B'] = 'orange'
EKPO_G['Color'][EKPO_G['ABC_Gruppe'] == 'C'] = 'green'

Jetzt Aggregieren wir über die ABC-Gruppe:

GruppenWerte = EKPO_G.groupby(['ABC_Gruppe'])
GruppenVolumen = (GruppenWerte.NETWR.sum()).to_frame()
GruppenMengen = (GruppenWerte.Preiseinheitsmenge.sum()).to_frame()

# Wieder zusammenfügen
GruppenVolumenMengen = GruppenVolumen.join(GruppenMengen)

Das Ergebnis:

GruppenVolumenMengen

Out:
NETWR Preiseinheitsmenge
ABC_Gruppe
A 6190725.01 175748.29
B 1231070.86 199599.24
C 408128.45 99745.63

Schauen wir uns nun die Verteilung der Werte und Mengen zwischen den Klassen A, B und C an:

GruppenVolumenMengen.plot(kind='bar', width=0.90, xlim=[0,1000], figsize=[10,5], yticks=GruppenVolumenMengen.NETWR)

 

abc_analyse_gruppen_vergleich

Es ist recht gut erkennbar, dass die Gruppe A deutlich mehr Umsatzvolumen (also Wertbeitrag) als die Gruppen B und C hat. Allerdings hat sie auch eine höhere Bestellmenge, wie jedoch nicht proportional von C über B zu A ansteigt wie das Umsatzvolumen.

Nachfolgend sehen wir die Klassifikation nochmal nicht kumuliert über die Umsatzvolumen der Materialien (Materialklassen):

EKPO_G[['NETWR']].plot(kind='bar', figsize=[20,10], legend = True, color=EKPO_G.Color, alpha=0.65, title='ABC - Analyse')

abc_analyse_sap_netwr

XYZ – Analyse

Für die XYZ-Analyse berechnen wir den arithmetischen Mittelwert, die Standardabweichung und die Summe aller Mengen pro Materialklasse [‘MATKL’] (oder alternativ, der einzelnen Materialnummern [‘MATNR’]) über eine Aggregation: 

Material_Menge = EKPO_X.Preiseinheitsmenge.groupby(EKPO_X.MATKL).agg({'mean', 'std', 'sum'})
#Oder mit dem Material: Material_Menge = EKPO_X.Preiseinheitsmenge.groupby(EKPO_X.MATNR) .agg({'mean', 'std', 'sum'})

#Leider ergeben sich einige NaNs bei der Standardabweichung, da ein Material oder eine Materialklasse nur eine einzige Buchung haben kann, diese müssen wir bereinigen (hier: mit Nullen auffüllen):
Material_Menge = Material_Menge.fillna(0)

Die XYZ-Analyse soll aufzeigen, welche Materialien (hier: Materialklassen) in stabilen Mengen verbraucht (hier: eingekauft) werden und welche größere Schwankungen hinsichtlich der Verbrauchsmenge (hier: Einkaufsmenge) aufweisen. Dazu berechnen wir den Variationskoeffizienten:

Variationskoeffizient = frac{Standardabweichung}{Mittelwert}

Wir berechnen diesen Variationskoeffizienten und sortieren das DataFrame nach diesem aufsteigend:

Material_Menge['Variationskoeffizient'] = Material_Menge['std'] / Material_Menge['mean']
Material_Menge = Material_Menge.sort_values(['Variationskoeffizient'], ascending = True)

Klassifikation:

Nun klassifizieren wir die Materialien (Materialklassen) über den Variationskoeffizienten in XYZ-Klassen. Dabei gehen wir davon aus, dass Materialien/Materialklassen, die einen Variationskoeffizienten von bis zu 70% des Maximalwertes aufweisen, in die Y-Klasse fallen. Solche, die nur maximal 20% des Maximalwertes aufweisen, fallen in die X-Klasse:

Material_Menge['XYZ_Gruppe'] = 'Z'
Material_Menge['XYZ_Gruppe'][Material_Menge.Variationskoeffizient <= Material_Menge.Variationskoeffizient.max() / 100 * 70] = 'Y'
Material_Menge['XYZ_Gruppe'][Material_Menge.Variationskoeffizient <= Material_Menge.Variationskoeffizient.max() / 100 * 20] = 'X'

Auch hier gilt analog zur ABC-Analyse: Intelligente Klassifikation erfolgt über die Analyse der Kurve der kumulierten Variationskoeffizienten. Die Grenzen der Klassen sollten idealerweise zwischen den steilsten Anstiegen (bzw. die größten Wertedifferenzen) zwischen den Werten der kumulierten Variationskoeffizienten-Liste gezogen werden.

Optional: Farben fürs Plotten setzen.

Material_Menge['Color'] = 'red'
Material_Menge['Color'][Material_Menge.XYZ_Gruppe == 'Y'] = 'orange'
Material_Menge['Color'][Material_Menge.XYZ_Gruppe == 'X'] = 'green'

Jetzt schauen wir uns mal die Verteilung der Materialien hinsichtlich des Variationskoeffizienten an:

Material_Menge.Variationskoeffizient.plot(kind='bar', width=0.90, xlim=[0,1000], figsize=[20,5], rot=90, color=Material_Menge.Color, title='XYZ - Analyse')

xyz_analyse_sap_matkl_menge

Die meisten Materialklassen haben einen recht niedrigen Variationskoeffizienten, sind im Einkauf (und daher vermutlich auch im Verbrauch) recht stabil. Die Materialklasse 0004 hingegen ist einigen Mengenschwankungen unterworfen. In der ABC-Analyse ist diese Materialklasse 0004 als B-Gruppe klassifiziert.

ABC-XYZ-Analyse

Nun möchten wir also die zuvor erstellte ABC-Klassifikation mit der XYZ-Klassifikation zusammen bringen.

Dafür fügen wir die beiden Pandas.DataFrame über den Index (hier die Materialklasse ‘MATKL’, im anderen Fall das Material ‘MATNR’) zusammen:

XYZ_ABC = pd.merge(EKPO_G, Material_Menge, left_index = True, right_index = True, how='left')

Die Zusammenfassung als Kreuztabelle:

pd.crosstab(XYZ_ABC.ABC_Gruppe, XYZ_ABC.XYZ_Gruppe, margins=True)

Out:

  X Y Z All

A 17 1 0 18

B 19 1 1 21

C 69 2 0 71

All 105 4 1 110

Für die Interpretation dieser Ergebnisse verweise ich erneut auf den Artikel bei der-wirtschaftsingenieur.de.

Benjamin Aunkofer

Benford-Analyse

/in , , , , /by

Das Benfordsche Gesetz beschreibt eine Verteilung der Ziffernstrukturen von Zahlen in empirischen Datensätzen. Dieses Gesetz, welches kein striktes Naturgesetz ist, sondern eher ein Erklärungsversuch in der Natur und in der Gesellschaft vorkommende Zahlenmuster vorherzusagen.

Das Benfordsche Gesetz beruht auf der Tatsache, dass die Ziffern in einem Zahlensystem hierarchisch aufeinander aufbauen: Es beginnt mit der 1, dann folgt die 2, dann die 3 usw. In Kombination mit bestimmten Gesetzen der Natur (der natürliche Wachstumsprozess, dabei möglichst energiesparend wachsen/überleben) oder Ökonomie (so günstig wie möglich einkaufen) ist gemäß des Benfordschen Gesetz zu erwarten, dass die Ziffer 1 häufiger vorkommt als die 2, die wiederum häufiger vorkommt als die 3. Die Ziffer 9 braucht demnach den längsten Weg und kommt entsprechend verhältnismäßig seltener vor.

Dieses Phönomen hilft uns bei echten Zufallszahlen nicht weiter, denn dann sind alle Ziffern nicht aufeinander aufbauend, sondern mehr oder weniger gleichberechtigt in ihrem Auftrauen. Bei der klassischen und axiomatischen Wahrscheinlichkeit kommen wir damit also nicht ans Ziel.

Die Benford-Analyse ist im Grunde eine Ausreißeranalyse: Wir vergleichen Ziffernmuster in Datenbeständen mit der Erwartungshaltung des Benfordschen Gesetzes. Weicht das Muster von dieser Erwartung ab, haben wir Diskussionsbedarf.

Moderne Zahlensysteme sind Stellenwert-Zahlensysteme. Neben den Dual-, Oktal- und Hexadezimalzahlensystemen, mit denen sich eigentlich nur Informatiker befassen, wird unser Alltag vom Dezimalzahlensystem geprägt. In diesem Zahlensystem hat jede Stelle die Basis 10 (“dezi”) und einen Exponenten entsprechend des Stellenwertes, multipliziert mit der Ziffer d. Es ist eine Exponentialfunktion, die den Wert der Ziffern in bestimmter Reihenfolge ermittelt:

Z =sum_{i=-n}^{m}d_{i}cdot10^{i}

Die Benford-Analyse wird meistens nur für die erste Ziffer (also höchster Stellenwert!) durchgeführt. Dies werden wir gleich einmal beispielhaft umsetzen.

Die Wahrscheinlichkeit des Auftretens der ersten “anführenden” Ziffer d ist ein Logarithmus zur Basis B. Da wir im alltäglichen Leben – wie gesagt – nur im Dezimalzahlensystem arbeiten, ist für uns B = 10.

p(d)=log_{B}left( 1 +frac{1}{d} right)

Im Standard-Python lässt sich diese Formel leicht mit einer Schleife umsetzen:

import math

[round(math.log10(1+1/float(i))*100.0, 2) for i in range(1,10)]
Out: [30.1, 17.61, 12.49, 9.69, 7.92, 6.69, 5.8, 5.12, 4.58]

Benford-Algorithmus in Python mit NumPy und Pandas

Nachfolgend setzen wir eine Benford-Analyse als Minimalbeispiel in NumPy und Pandas um.

import numpy as np

x = np.arange(1,10) # NumPy-Array erstellen
x
Out: array([ 1,  2,  3,  4,  5,  6,  7,  8,  9])

benford = np.round(np.log10(1+1/x) * 100.0, decimals=2) # Den Logarithmus auf das NumPy-Array anwenden und runden
Out: array([ 30.1 ,  17.61,  12.49,   9.69,   7.92,   6.69,   5.8 ,   5.12,
         4.58])

Nun möchten wir eine Tabelle erstellen, mit zwei Spalten, eine für die Ziffer (Digit), die andere für die relative Häufigkeit der Ziffer (Benford Law). Dazu nutzen wir das DataFrame aus dem Pandas-Paket. Das DataFrame erstellen wir aus den zwei zuvor erstellten NumPy-Arrays.

import pandas as pd

benfordFrame = pd.DataFrame({'Digit': x, 'Benford Law': benford})

benfordFrame
Out: 
   Benford Law  Digit
0        30.10      1
1        17.61      2
2        12.49      3
3         9.69      4
4         7.92      5
5         6.69      6
6         5.80      7
7         5.12      8
8         4.58      9

Man könnte sicherlich auch den natürlichen Index des DataFrames nutzen, indem wir diesen nur um jeweils 1 erhöhen, aber das verwirrt später nur und tun wir uns jetzt daher lieber nicht an…

benfordFrame.plot('Digit', 'Benford Law', kind='bar', title='Benford', legend=False)

Das Dataframe-Objekt kann direkt plotten (das läuft über die matplotlib, die wir aber nicht direkt einbinden müssen):benfordsches_gesetz_dezimalzahlen_ziffern

Diese neun Balken zeigen die Verteilung der Ziffernhäufigkeit nach dem Benfordschen Gesetz, diese Verteilung ist unsere Erwartungshaltung an andere nummerische Datenbestände, wenn diese einem natürlichen Wachstum unterliegen.

Analyse der Verteilung der ersten Ziffer in Zahlungsdaten

Jetzt brauchen wir Daten mit nummerischen Beträgen, die wir nach Benford testen möchten. Für dieses Beispiel nehme ich aus einem ein SAP-Testdatensatz die Spalte ‘DMBTR’ der Tabelle ‘BSEG’ (SAP FI). Die Spalte ‘DMBTR’ steht für “Betrag in Hauswährung’, die ‘BSEG’ ist die Tabelle für die buchhalterischen Belegsegmente.
Die Datei mit den Testdaten ist über diesen Link zum Download verfügbar (Klick) und enthält 40.000 Beträge.

Wir laden den Inhalt der Datei via NumPy.LoadTxt() und machen aus dem resultierenden NumPy-Array wieder ein Pandas.DataFrame und holen uns die jeweils erste Ziffer für alle Einträge als Liste zurück.

financialTransactions = np.loadtxt("[DEIN LOKALER PFAD]\BSEG_DMBTR.csv", skiprows=1)

financialTransactionsFrame = pd.DataFrame({'Zahlungen':financialTransactions})

firstDigits = [str(value)[0:1] for value in financialTransactionsFrame['Zahlungen']]

Die Einträge der ersten Ziffer in firstDigits nehmen wir uns dann vor und gruppieren diese über die Ziffer und ihrer Anzahl relativ zur Gesamtanzahl an Einträgen.

percentDigits = np.asarray([[i, firstDigits.count(str(i))/float(len(financialTransactionsFrame['Zahlungen']))*100] for i in range(1, 10)])

percentDigits.T[1].sum()
Out: 94.0

Wenn wir die Werte der relativen Anzahl aufsummieren, landen wir bei 94% statt 100%. Dies liegt daran, dass wir die Ziffer 0 ausgelassen haben, diese jedoch tatsächlich vorkommt, jedoch nur bei Beträgen kleiner 1.00. Daher lassen wir die Ziffer 0 außenvor. Wer jedoch mehr als nur die erste Ziffer prüfen möchte, wird die Ziffer 0 wieder mit in die Betrachtung nehmen wollen. Nur zur Probe nocheinmal mit der Ziffer 0, so kommen wir auf die 100% der aufsummierten relativen Häufigkeiten:

percentDigits = np.asarray([[i, firstDigits.count(str(i))/float(len(financialTransactionsFrame['Zahlungen']))*100] for i in range(0, 10)])
percentDigits.T[1].sum()
Out: 100.0

Nun erstellen wir ein weiteres Pandas.DataFrame, mit zwei Spalten: Die Ziffern (Digit) und die tatsächliche Häufigkeit in der Gesamtpopulation (Real Distribution):

percentDigitsFrame = pd.DataFrame({'Digit':percentDigits[:,0], 'Real Distribution':percentDigits[:,1]})

Abgleich der Ziffernhäufigkeit mit der Erwartung

Nun bringen wir die theoretische Verteilung der Ziffern, also nach dem zuvor genannten Logarithmus, und die tatsächliche Verteilung der ersten Ziffern in unseren Zahlungsdaten in einem Plot zusammen:

import matplotlib.pyplot as pyplot 

fig = pyplot.figure()

ax = fig.add_subplot(111)
ax2 = ax.twinx()

percentDigitsFrame.plot('Digit', 'Real Distribution', kind='bar', ax=ax2, width = 0.4, color="green", position=0, legend=False)
benfordFrame.plot('Digit', 'Benford Law', kind='bar', ax=ax, width = 0.4, color="blue", position=1, legend=False)

lines, labels = ax.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2, loc=0)

ax.set_ylim(1,35)
ax2.set_ylim(1,35)

pyplot.show()

In dem Plot wird deutlich, dass die Verteilung der führenden Ziffer in unseren Zahlungsdaten in ziemlich genau unserer Erwartung nach dem Benfordschen Gesetz entspricht. Es sind keine außerordentlichen Ausreißer erkennbar. Das wäre auch absolut nicht zu erwarten gewesen, denn der Datensatz ist mit 40.000 Einträgen umfassend genug, um dieses Muster gut abbilden zu können und von einer Manipulation dieser Beträge im SAP ist ebenfalls nicht auszugehen.

benford-analyse-python-fuer-sap-bseg-dmbtr

Gegenüberstellung: Computer-generierte Zufallszahlen

Jetzt wollen wir nochmal kurz darauf zurück kommen, dass das Benfordsche Gesetz für Zufallszahlen nicht unwendbar ist. Bei echten Zufallszahlen bin ich mir da auch sehr sicher. Echte Zufallszahlen ergeben sich beispielweise beim Lotto, wenn die Bälle mit Einzel-Ziffern durch eine Drehkugel hüpfen. Die Lottozahl-Ermittlung erfolgt durch die Zusammenstellung von jeweils gleichberechtigt erzeugten Ziffern.
Doch wie ist dies bei vom Computer generierten Zufallszahlen? Immerhin heißt es in der Informatik, dass ein Computer im Grunde keine Zufallszahlen erzeugen kann, sondern diese via Takt und Zeit erzeugt und dann durchmischt. Wir “faken” unsere Zahlungsdaten nun einfach mal via Zufallszahlen. Hierzu erstellen wir in NumPy ein Array mit 2000 Einträgen einer Zufallszahl (NumPy.Random.rand(), erzeugt floats 0.xxxxxxxx) und multiplizieren diese mit einem zufälligen Integer (Random.randint()) zwischen 0 und 1000.

from random import randint

financialTransactions = np.round(np.random.rand(2000) * randint(0,1000),decimals=2) # Zufallszahlen erzeugen, die auf dem ersten Blick als Zahlungsdaten durchgehen :-)

Erzeugen wir die obigen Datenstrukturen erneut, zeigt sich, dass die Verteilung der Zufallszahlen ganz anders aussieht: (vier unterschiedliche Durchläufe)

benford-analyse-python-numpy-pseudo-zufallszahlen

Anwendung in der Praxis

Data Scientists machen sich das Benfordsche Gesetz zu Nutze, um Auffälligkeiten in Zahlen aufzuspüren. In der Wirtschaftsprüfung und Forensik ist diese Analyse-Methode recht beliebt, um sich einen Eindruck von nummerischen Daten zu verschaffen, insbesondere von Finanztransaktionen. Die Auffälligkeit durch Abweichung vom Benfordchen Gesetz entsteht z. B. dadurch, dass Menschen eine unbewusste Vorliebe für bestimmte Ziffern oder Zahlen haben. Greifen Menschen also in “natürliche” Daten massenhaft (z. B. durch Copy&Paste) ein, ist es wahrscheinlich, dass sie damit auch vom Muster des Benfordschen Gesetzes abweichen. Weicht das Muster in Zahlungsströmen vom Bendfordsche Erwartungsmuster für bestimmte Ziffern signifikant ab, könnte dies auf Fälle von unnatürlichen Eingriffen hindeuten.

Die Benford-Analyse wird auch gerne eingesetzt, um Datenfälschungen in wissenschaftlichen Arbeiten oder Bilanzfälschungen aufzudecken. Die Benford-Analyse ist dabei jedoch kein Beweis, sondern liefert nur die Indizien, die Detailanalysen nach sich ziehen können/müssen.

Benjamin Aunkofer

Die Rastrigin-Funktion

/in , , , , /by

Jeder Data Scientist kommt hin und wieder mal in die Situation, einen Algorithmus trainieren bzw. optimieren zu wollen oder zu müssen, ohne jedoch, dass passende Trainingsdaten unmittelbar verfügbar wären. Zum einen kann man in solchen Fällen auf Beispieldaten zugreifen, die mit vielen Analysetools mitgeliefert werden, oder aber man generiert sich seine Daten via mathematischer Modelle selbst, die für bestimmte Eigenschaften bekannt sind, die gute Bedingungen für das Optimierungstraining liefern. 

Ein solches Modell, das man als Machine Learning Entwickler kennen sollte, ist die Rastrigin-Funktion, die laut Wikipedia von Leonard A. Rastrigin erstmalig beschrieben wurde. Dabei handelt es sich um eine Häufigkeites-/Wahrscheinlichkeitsverteilung, deren Dichte mehrere lokale Modi (Gipfel) aufweist. Ein Modus (oder Modalwert) ist in einer Häufigkeitsverteilung der häufigste Wert (“Bergspitze”) bzw. der Wert mit der höchsten Wahrscheinlichkeit.

Anmerkung des Autors: Dieser Artikel stellt zum einen die Rastrigin-Funktion und ihre Bedeutung für die Optimierungsrechnung vor, ist zum anderen aber auch eine Einführung in den Umgang mit NumPy-Matrizen (die eine Menge For-Schleifen ersparen können).

Die Rastrigin-Funktion

Mathematisch beschrieben wird die Rastrigin-Funktion wie folgt:

f(x_1 cdots x_n) = An + sum_{i=1}^n (x_i^2 -10cos(2pi x_i))

-5.12 leq x_i leq 5.12

Wobei für das globale Minimum gilt: f(0) = 0
Außerdem ist zu beachten, dass A=10 eine Konstante ist.

Die Rastrigin-Funktion im Standard-Python umsetzen und visualisieren

Die Formel lässt sich in Python (wie natürlich in jeder anderen Programmiersprache auch) einfach umsetzen:

value = 10 + x**2 - 10 * math.cos(2 * math.pi * x)

Nun können wir über den klassischen Weg der Programmierung einfach eine For-Schleife verwenden, um die Rastrigin-Funktionswerte in eine Liste zu packen und mit einem Plot zu visualsieren, dabei bin ich leider doch nicht ganz um die Verwendung des NumPy-Pakets nicht herumgekommen:

import matplotlib.pyplot as pyplot
import numpy as np # NumPy hat die Matrizen-Datenstruktur, die wir benötigen 
import math as math # Grundlegende mathematische Funktionen (hier benötigt: Kreiszahl Pi und Cosinus-Funktion)

rastriginValues = []
i = 0


for x in np.arange(-5.12, 5.12, 0.01): # Die Python-eigene range()-Funktion kann leider keine Floats, sondern nur Integer erzeugen :-/
    value = 10 + x**2 - 10 * math.cos(2 * math.pi * x)
    i += 1
    print(i, x, value)
    rastriginValues.append(value)
    
pyplot.plot(rastriginValues)
pyplot.ylim(0,50)
pyplot.xlim(0,1024)
pyplot.show()

rastrigin-line-chart

Die grafische Darstellung zeigt, dass es sich tatsächlich um eine symmetrische multimodalen Verteilung handelt.

Die Rastrigin-Funktion mehrdimensional umsetzen, mit NumPy-Matrizen-Funktionen

Die obige Umsetzung der Rastrigin-Funktion ist eindimensional (eine Variable), braucht für die Darstellung allerdings zwei Dimensionen (f(x) und die Durchlaufanzahl bzw. Zeitachse). Nun könnten wir die Zahl der Variablen von 1 (x) auf 2 (x und y) erhöhen und eine dreidimensionale Darstellung erzeugen. Eine ähnliche dreidimensionale Darstellung gab es bereits in meiner Vorstellung des k-nearest-Neighbour-Algorithmus nachzuvollziehen. Dabei müssten wir die Konstante A=10 auf A=20 verdoppeln:

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as pyplot
import numpy as np

figure = pyplot.figure()
axe = figure.add_subplot(111, projection='3d')

x = np.linspace(-5.12, 5.12, 100) # unterteilt den Bereich in 100 Schnitte, ähnlich: np.arange(-5.12, 5.12, 0.1)
y = np.linspace(-5.12, 5.12, 100)
x, y = np.meshgrid(x, y) # erzeugt ein Koordinatensystem

# Nun ohne Schleifen: Wir wenden die NumPy-Funktionen (np.cos statt math.cos und np.pi statt math.pi) 
# auf die NumPy-Arrays an (x und y) und erhalten ein NumPy-Array z zurück
z = 20 + x**2 - 10 * np.cos(2 * np.pi * x) + y**2 - 10 * np.cos(2* np.pi * y)

# Plotte die drei Variablen (x, y, z) im dreidimensionalen Raum
axe.plot_surface(x, y, z, rstride=1, cstride=1, cmap="jet", linewidth=0, antialiased=False)

pyplot.title('Rastrigin-Map')
pyplot.grid(True)
axes = pyplot.gca()
axes.set_xlim([-5.12,5.12])
axes.set_ylim([-5.12,5.12])
pyplot.show()

rastrigin-funktion-python-3d-plot

Die Rastrigin-Funktion wird gerne für Optimierungsalgorithmen eingesetzt, wofür sie wegen des großen Suchraums und der hohen Anzahl lokaler Modi ein herausforderndes Umfeld bietet. Beispielsweise wird – meines Erachtens nach – das wohl beliebteste Optimierungsverfahren im maschinellen Lernen, das Gradientenverfahren, hier keine guten Ergebnisse liefern, denn es gibt einfach zu viele lokale Minima.

 

 

Benjamin Aunkofer

Python vs R Statistics

/in , , , /by

Immer wieder wird mir von Einsteigern die Frage gestellt, ob sich der Einstieg und die Einarbeitung in die Programmiersprache Python eher lohnen würde als in R Statistics. Nun gibt es in den englischsprachigen Portalen bereits viele Diskussionen und Glaubenskriege zu diesem Vergleich – diese habe ich mir mit Absicht nicht weiter durchgelesen, sondern ich versuche hier meine Erfahrung aufs Blog zu bringen und bin auf Eure Meinungen/Erfahrungen gespannt!

Mit weniger R-Code schneller zum Ziel, und mit Python darüber hinaus

Was mir beim Einstieg in R gleich auffiel: Nach der Installation kann man sofort loslegen! Ein Plot oder eine Regressionsanalyse ist binnen weniger Code-Zeilen erledigt, denn die Sprache bringt diese Funktionen von Haus aus mit. In Python ist das Ziel auch nicht weit weg, allerdings müssen für die Plots erst die MatplotLib installiert werden, für Matrizenberechnung die Numpy-Bibliothek und um eine, mit der R-Datenstruktur Data.Frame vergleichbare Datenstruktur in Python zu erhalten, die Pandas-Bibliothek. Diese Python-Bibliotheken kann man zwar mit Fug und Recht als Bestandteil des Python-Universums ansehen, standardmäßig ausgeliefert werden sie aber nicht und auch sollten sie streng vom Standardpython in der Anwendung getrennt werden, im Klartext: Die Bibliotheken erfordern extra Einarbeitung und machen die Handhabung komplizierter, das einfache Python verliert ein Stück weit seine Einfachheit.

Auch die beliebte Entwicklungsumgebung R-Studio sucht seinesgleichen und ist IPython meiner Meinung nach hinsichtlich der Usability absolut überlegen. R ist einfach darauf ausgerichtet, Daten zu analysieren und zu visualisieren, aber beschränkt sich eben auch darauf.

“R is more about sketching, and not building,” says Michael Driscoll, CEO of Metamarkets. “You won’t find R at the core of Google’s page rank or Facebook’s friend suggestion algorithms. Engineers will prototype in R, then hand off the model to be written in Java or Python.”

Im Gegenzug ist Python eine Programmiersprache, die nicht nur an den einen Zweck gebunden ist. Mit Python können ebenfalls (Web-)Server- oder Desktop-Anwendungen und somit ohne Technologiebruch analytische Anwendungen komplett in Python entwickelt werden. Und auch wenn R ebenfalls unüberschaubar viele Packages mitbringt, bietet Python noch einiges mehr, beispielsweise zur dreidimensionalen Darstellung von Graphen.

Software-Entwickler lieben Python, Mathematiker eher R

Data Science ist ein äußerst interdisziplinäres Fachgebiet und Data Scientists können Mathematiker, Physiker, Informatiker, Ingenieure oder (wenn auch etwas seltener) Wirtschafts- oder auch Geisteswissenschaftler sein. Ein Großteil kommt aus der Mathematik oder äußerst mathematischer Fachgebiete wie der Physiker oder der Elektroingenieurwissenschaft. In diesen Studiengängen wird überwiegend mit Programmiersprachen gearbeitet, die von Mathematikern für Mathematiker entwickelt wurden, also R Statistics, MATLAB oder Octave. Beispielsweise ist meine Frau studierte Elektotechnikingenieurin und setzte alle ihre Prototypen des maschinellen Lernens in MATLAB um, sie findet sich aber auch in R gut zurecht.

Wer aus der Software-Entwicklung kommt, findet sich in Python vermutlich sehr viel schneller zurecht als in R. In meiner subjektiven Wahrnehmung stelle ich tatsächlich fest, dass diejenigen Data Scientists, die aus der Mathematik zum Data Science gekommen sind, meistens R präferieren und diejenigen, die aus der Anwendungsentwicklung kommen, eher mit Python arbeiten.

datascience

Python kollaboriert besser

Ein Data Scientist kommt selten allein, denn Data Science ist Teamarbeit. Und wo Teams ein gemeinsames Ziel erreicht sollen, werden besondere Anforderungen an die Arbeitsumgebung gestellt. Python gilt als eine syntaktisch leicht verständliche Programmiersprache, die manchmal sogar als “executable Pseudocode” bezeichnet wird (was allerdings dann doch leicht übertrieben ist…). Es ist also für alle Teammitglieder eine relativ einfach zu erlernende Sprache. Dabei muss Python nicht von allen Teammitgliedern favorisiert werden, denn eigene lokale Prototypen können in R, Octave oder was auch immer erstellt werden, lassen sich dann aber auch einfach in Python integrieren. Für richtig schnelle Anwendungen sind Python und R als Interpretersprachen sowieso zu langsam, solche Anwendungen werden am Ende in C/C++ umgesetzt werden müssen, aber selbst dann bietet Python nicht zu unterschätzende Vorteile: Der Erfolg von Python im wissenschaftlichen Rechnen beruht nämlich auch auf der unkomplizierten Integration von Quellcode der Programmiersprachen C, C++ und Fortran.

Neue Spieler auf dem Feld: Scala und Julia

Leider kann ich zu den beiden Programmiersprachen Scala und Julia (noch) nicht viel sagen. Scala scheint sich meiner Einschätzung nach als eine neue Alternative für Python zu entwickeln. Scala ist ein Produkt aus dem Java-Universum und war als eine Programmiersprache für unterschiedlichste Zwecke gedacht. Die Sprache setzt sich im Big Data Science immer weiter durch, einige Tools für Big Data Analytics (Apache Spark, Apache Flink) sind auf Scala ausgelegt und basieren selbst auf dieser Programmiersprache. Was Scala als eine stark von Java inspirierte Sprache sehr sympathisch macht, ist der enorm kompakte Code. Ein MapReduce-Algorithmus lässt sich in Scala mit einem Bruchteil an Code erstellen, als es in Java der Fall wäre, wie es auch die Code-Beispiele der Spark-Webseite eindrücklich zeigen: (Was ist eigentlich Apache Spark?)

Text Search in Python (Apache Spark)

textFile = sc.textFile("hdfs://...")

# Creates a DataFrame having a single column named "line"
df = textFile.map(lambda r: Row(r)).toDF(["line"])
errors = df.filter(col("line").like("%ERROR%"))
# Counts all the errors
errors.count()
# Counts errors mentioning MySQL
errors.filter(col("line").like("%MySQL%")).count()
# Fetches the MySQL errors as an array of strings
errors.filter(col("line").like("%MySQL%")).collect()
Text Search in Scala (Apache Spark)

val textFile = sc.textFile("hdfs://...")

// Creates a DataFrame having a single column named "line"
val df = textFile.toDF("line")
val errors = df.filter(col("line").like("%ERROR%"))
// Counts all the errors
errors.count()
// Counts errors mentioning MySQL
errors.filter(col("line").like("%MySQL%")).count()
// Fetches the MySQL errors as an array of strings
errors.filter(col("line").like("%MySQL%")).collect()
Text Search in Java (Apache Spark)

// Creates a DataFrame having a single column named "line"
JavaRDD<String> textFile = sc.textFile("hdfs://...");
JavaRDD<Row> rowRDD = textFile.map(
  new Function<String, Row>() {
    public Row call(String line) throws Exception {
      return RowFactory.create(line);
    }
  });
List<StructField> fields = new ArrayList<StructField>();
fields.add(DataTypes.createStructField("line", DataTypes.StringType, true));
StructType schema = DataTypes.createStructType(fields);
DataFrame df = sqlContext.createDataFrame(rowRDD, schema);

DataFrame errors = df.filter(col("line").like("%ERROR%"));
// Counts all the errors
errors.count();
// Counts errors mentioning MySQL
errors.filter(col("line").like("%MySQL%")).count();
// Fetches the MySQL errors as an array of strings
errors.filter(col("line").like("%MySQL%")).collect();

Julia wurde (ähnlich wie R) explizit für den Zweck der statistischen Datenanalyse entwickelt, wird auf Grund des aktuellen Beta-Status noch kaum produktiv eingesetzt. Da Julia auf sehr schnelle Anwendungen ausgerichtet ist, liegt in Julia die neue Hoffnung für jene, für die R und Python zu langsame Interpretersprachen sind.

Buchempfehlungen zum Einstieg in R oder Python

Es versteht sich von selbst, dass ich alle Bücher auch selbst besitze und mehr als nur das Vorwort gelesen habe…

Was ist Eure Erfahrung? Ihr seid gefragt!

Schreibt Eure Meinung einfach als Kommentar zu diesem Artikel! Wer meint, den Vergleich logischer, “richtiger” und nachvollziehbarer aufs digitale Papier bringen zu können, darf einen Artikelvorschlag übrigens gerne an redaktion@data-science-blog.com senden!