Tutorial Archives

R Data Frames meistern mit dplyr – Teil 1

November 21, 2016/in Data Mining, Data Science, Data Science Hack, R Statistics, Statistics, Tutorial/by Dr. Uwe Schneider

Dieser Artikel ist Teil 1 von 2 aus der Artikelserie R Data Frames meistern mit dplyr.

Data Frames sind das Arbeitspferd von R, wenn Daten in eine Struktur gepackt werden sollen, um sie einzulesen, zu säubern, zu transformieren, zu analysieren und zu visualisieren. Abstrakt gesprochen sind Data Frames nichts anderes als Relationen, also Mengen von Tupels, gebildet aus Elementen von geeigneten Mengen.

Dieses Konzept hat sich auch außerhalb des R-Universums bestens bewährt, umzusammengesetzte Daten, Beobachtungen oder Geschäftsobjekte zu repräsentieren. Der beste Beleg für diese Aussage sind die allgegenwärtigen Relationalen Datenbanksysteme (RDBMS). Dort werden Relationen als Tabellen (Tables) oder Sichten (Views) bezeichnet, und darauf wirkt eine mächtige, imperative Abfrage- und Manipulationssprache namens Structured Query Language, kurz:
SQL.

SQL ist in meiner Wahrnehmung die Lingua Franca der Datenverarbeitung, da sie im Kern über sehr viele Softwareprodukte gleich ist und nach erstaunlich geringem Lernaufwand mächtige Auswerte- und Manipulationsoperationen an den Daten ermöglicht. Hier eine SQL-Anweisung, um eine fiktive Tabelle aller Verkäufe (SALES) nach den Top-10-Kunden in diesem Jahr zu untersuchen:

SELECT customer_id, SUM(sales_amt) AS amount,

COUNT(*) AS n_sales, COUNT(DISTINCT product_id) AS n_products

FROM sales

WHERE sales_date LIKE ’2016-%’

GROUP BY customer_id

ORDER BY SUM(sales_amt) DESC

LIMIT 10

Dieser selbsterklärliche Code aus sieben Zeilen hat einen enormen Effekt: Er fast alle Verkäufe des Jahres 2016 auf Basis der Kundennummer zusammen, berechnet dabei die Summe aller Verkaufsbeträge, zählt die Anzahl der Transaktionen und der verschiedenen vom Kunden gekauften Produkte. Nach Sortierung gemäß absteigenden Umsatzes schneidet der Code nach dem 10. Kunden ab.

SQL kann aber mit der gleichen Eleganz noch viel mehr: Beispielsweise verbinden Joins die Daten mehrerer Tabellen über Fremschlüsselbeziehungen oder analytische Funktionen bestimmen Rankings und laufende Summen. Wäre es nicht toll, wenn R ähnlich effektiv mit Data Frames analoger Struktur umgehen könnte? Natürlich! Aber schon der Versuch, obige SQL-Query auf einem R Data Frame mit den althergebrachten Bordmitteln umzusetzen (subset, aggregate, merge, …), führt zu einem unleserlichen, uneleganten Stück Code.

Genau in diese Bresche springt der von vielen anderen Bibliotheken bekannte Entwickler Hadley Wickham mit seiner Bibliothek dplyr: Sie standardisiert Operationen auf Data Frames analog zu SQL-Operationen und führt zu einer wirklich selbsterklärlichen Syntax, die noch dazu sehr performant abgearbeitet wird. Ganz analog zu ggplot2, das sich an der Grammar of Graphics orientiert, spricht Wickham bei dplyr von einer Grammar of Data Manipulation. Die Funktionen zur Manipulation nennt er folgerichtig Verben.

Dabei treten naturgemäß eine Reihe von Analogien zwischen den Teilen eines SELECT-Statements und dplyr-Funktionen auf:

SELECT-Operation	dplyr-Funktion
Bildung der Spaltenliste	select()
Bildung eines Ausdrucks	mutate()
WHERE-Klausel	filter()
GROUP BY Spaltenliste	group_by()
Bildung von Aggregaten wie sum() etc.	summarise()
HAVING-Klausel	filter()
ORDER BY Spaltenliste	arrange()
LIMIT-Klausel	slice()

Die ersten Schritte

Ich möchte die Anwendung von dplyr mithilfe des Standard-Datensatzes Cars93
aus dem Paket MASS demonstrieren:

> install.packages("dplyr")

> library(dplyr)

> cars <- MASS::Cars93

> str(cars)

’data.frame’: 93 obs. of 27 variables:

$ Manufacturer : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 ...

$ Model : Factor w/ 93 levels "100","190E","240",..: 49 56 ...

Die erste Aufgabe soll darin bestehen, aus dem Data Frame alle Autos zu selektieren, die vom Hersteller “Audi” stammen und nur Model und Anzahl Passagiere auszugeben. Hier die Lösung in Standard-R und mit dplyr:

> # Standard R: Zu allen Audis das Model, Anzahl Insassen

> subset(cars,Manufacturer=="Audi")[,c("Model","Passengers")]

Model Passengers

3 90 5

4 100 6

> # Das gleiche mit dplyr

> select(filter(cars,Manufacturer=="Audi"),Model,Passengers)

Model Passengers

1 90 5

2 100 6

Man sieht, dass die neue Funktion filter() der Zeilenselektion, also der Funktion subset() entspricht. Und die Auswahl der Ergebnisspalten, die in Standard-R durch Angabe einer Spaltenliste zwischen [ und ] erfolgt, hat in dplyr das Pendant in der Funktion select().

select() ist sehr mächtig in seinen Möglichkeiten, die Spaltenliste anzugeben. Beispielsweise funktioniert dies über Positionslisten, Namensmuster und ggf. das auch noch negiert:

1 2	select(cars, -starts_with("L"))

Die obige Abfrage projiziert aus dem Data Frame sämtliche Spalten, die nicht mit “L” beginnen. Das scheint zunächst ein unscheinbares Feature zu sein, zahlt sich aber aus, wenn analytische Data Frames Dutzende oder Hunderte von Spalten haben, deren Bezeichnung sich nach einem logischen Namensschema richtet.
Soweit ist das noch nicht spektakulär. dplyr hilft uns in obigem Beispiel, als erstes bestimmte Datensätze zu selektieren und als zweites die interessierenden Spalten zu projizieren. dplyr ist aber bezüglich der Verarbeitung von Data Frames sehr intuitiv und funktional, sodass wir früher oder später viele Operationen auf unserem Data Frame verketten werden. So erreichen wir die Mächtigkeit von SQL und mehr. Die funktionale Syntax aus dem letzten Beispiel wird dann ganz schnell unleserlich, da die Verabeitungsreihenfolge (zuerst filter(), dann select()) nur durch Lesen des Codes von innen nach außen und von rechts nach links ersichtlich wird.

Daher geht dplyr einen Schritt weiter, indem es den eleganten Verkettungsoperator %>% aus dem magrittr-Paket importiert und zur Verfügung stellt. Dadurch werden die verschachtelten Ausdrücke in Sequenzen von Operationen gewandelt und somit sehr viel lesbarer und wartbarer:

select(filter(cars,Manufacturer=="Audi"),Model,Passengers)

Model Passengers

1 90 5

2 100 6

> cars %>%

+ filter(Manufacturer=="Audi") %>%

+ select(Model,Passengers)

Model Passengers

1 90 5

2 100 6

Diese in meinen Augen geniale Syntax durch den neuen Operator %>% erlaubt einen sequenziellen Aufbau der Operationen auf einem Data Frame. Benutzer der Unix-Kommandozeile werden hier leicht die Analogie zu Pipes erkennen. Ganz abstrakt kann man sagen, dass damit folgende Operationen äquivalent sind:

Traditioneller Funktionsaufruf	Verkettung mit %>%
f(a,b)	a %>% f(b)
f(a,b,c)	a %>% f(b,c)
g(f(a,b),c)	a %>% f(b) %>% g(c)

Weiteres erklärt die Dokumentation zum %>%-Operator im Paket magrittr mithilfe
des Befehls ?magrittr::‘%>%‘.

Neue Variablen

Durch die Funtionen select() und filter() können wir aus Data Frames Spalten projizieren und Zeilen selektieren. Ergebnisse neuer Ausdrücke entstehen hingegen mit dem Verb mutate():

> # Mit Umrechnung der Mileage in Verbrauch und des

> # Kaufpreises von USD in EUR

> cars2 <-

+ cars %>%

+ filter(Manufacturer=="Audi") %>%

+ mutate(l_100km = 235 / MPG.city,

+ eur = Min.Price * 1000 / 1.1) %>%

+ select(Model,Passengers,l_100km,eur)

> cars2

Model Passengers l_100km eur

1 90 5 11.75000 23545.45

2 100 6 12.36842 28000.00

> class(cars2)

[1] "data.frame"

Im obigen Beispiel wird zunächst auf den Hersteller Audi selektiert und danach auf einen Streich zwei neue Spalten eingeführt, l_100km und eur. Durch Zuweisen auf eine neue Variable wird das fertige Ergebnis dauerhaft gespeichert. Hierbei handelt es sich wieder um ein natives Data Frame-Objekt. Die Operation transmute() arbeitet analog zu mutate(), verwirft aber nach Bildung der Ausdrücke alle nicht genannten Spalten. Somit können wir obiges Beispiel auch wie folgt schreiben:

cars3 <-

+ cars %>%

+ filter(Manufacturer=="Audi") %>%

+ transmute(Model = Model, Passengers = Passengers,

+ l_100km = 235 / MPG.city,

+ eur = Min.Price * 1000 / 1.1)

cars3

Model Passengers l_100km eur

1 90 5 11.75000 23545.45

2 100 6 12.36842 28000.00

Aggregate

Neben der Selektion von Zeilen und Spalten sowie der Bildung abgeleiteter Ausdrücke ist bei Datenbanktabellen die Gruppierung und Aggregation mit GROUP BY eine sehr wichtige Operation. Dies gilt auch für Data Frames in R, wenngleich hier der Funktionsumfang über diverse Funktionen wie table() oder aggregate() verteilt ist und wenig intuitiv ist.

Hier bringt dplyr ebenfalls eine großartige Verbesserung mit. Das entsprechende Verb heißt group_by(). Diese Operation wird zusammen mit einer Spaltenliste auf ein Data Frame angewendet:

grp_cars <- cars %>% group_by(Manufacturer)

> class(grp_cars)

[1] "grouped_df" "tbl_df" "tbl" "data.frame"

Das Ergebnis von group_by() ist ein Objekt, das “mehr” ist als ein Data Frame, sondern auch noch einige spezifische Strukturinformationen von dplyr enthält. In unserem Beispiel sind dies Indizes von Zeilen, die zum gleichen Hersteller gehören. Das ursprüngliche Data Frame wird hierbei nicht kopiert, sondern nur eingebettet.

Nach Anwenden einer group_by()-Operation ist das Data Frame optimal vorbereitet für die eigentliche Aggregation mit summarise():

> sum_cars <- cars %>%

+ group_by(Manufacturer) %>%

+ summarise(nCars = n(), avg_price = mean(Min.Price))

> str(sum_cars)

Classes ‘tbl_df’, ‘tbl’ and ’data.frame’: 32 obs. of 3 variables:

$ Manufacturer: Factor w/ 32 levels "Acura","Audi",..: 1 2 3 4 5 6 7 8 9 10 ...

$ nCars : int 2 2 1 4 2 8 1 2 6 2 ...

$ avg_price : num 21.1 28.4 23.7 20.8 35.2 ...

Das Resultat von summarise() ist wieder ein Data Frame, das neben den ursprünglichen Gruppierungskriterien nur noch die Aggregate enthält.

Daten in Reih’ und Glied

Zwischen Relationalen Datenbanken und R-Data Frames besteht ein wesentlicher konzeptioneller Unterschied: Die Ergebnisse eines SELECT-Befehls haben keine definierte Reihenfolge, so lange die Zeilen nicht mit der Klausel ORDER BY festgelegt wird. Im Gegensatz dazu haben die Zeilen von Data Frames eine konstante Reihenfolge, die sich aus der Anordnung derWerte in den Spaltenvektoren ergibt.

Dennoch ist es manchmal wünschenswert, Data Frames umzusortieren, um eine fachliche Reihenfolge abzubilden. Hierzu dient in dplyr das Verb arrange(), das im Standard-R weitgehend der Indizierung eines Data Frames mit Ergebnissen der order()-Funktion entspricht, aber syntaktisch eleganter ist:

cars %>%

+ arrange(desc(Horsepower),Manufacturer) %>%

+ select(Manufacturer, Model, Horsepower,Min.Price) %>%

+ slice(1:5)

Manufacturer Model Horsepower Min.Price

1 Chevrolet Corvette 300 34.6

2 Dodge Stealth 300 18.5

3 Cadillac Seville 295 37.5

4 Infiniti Q45 278 45.4

5 Mazda RX-7 255 32.5

Dieses Beispiel hat zum Ziel, die fünf PS-stärksten Autos zu selektieren. Die arrange()-Funktion sortiert hier zunächst absteigend nach der PS-Stärke, dann aufsteigend nach Herstellername. Die Selektion der 5 ersten Zeilen erfolgt mit der hilfreichen Funktion slice(), die aus einem Data Frame Zeilen anhand ihrer Reihenfolge selektiert.

Fazit und Ausblick

Mit dplyr wird die Arbeit mit Data Frames stark verbessert: Im Vergleich zu “nacktem” R bringt das Paket eine klarere Syntax, abgerundete Funktionalität und bessere Performance. In der Kürze dieses Artikels konnte ich dies nur oberflächlich anreissen. Daher verweise ich auf die vielen Hilfe-Seiten, Vignetten und Internet-Videos zum Paket. Im zweiten Teil dieses Artikels werde ich auf einige fortgeschrittene Features von dplyr eingehen, z.B. die Verknüpfung von Data Frames mit Joins, die Window-Funktionen und die Verwendung von Datenbanken als Backend.

Weiter zu R Data Frames meistern mit dplyr – Teil 2.

ABC-XYZ-Analyse

October 20, 2016/in Business Analytics, Business Intelligence, Data Mining, Data Science, Data Science Hack, Python, Statistics, Tutorial/by Benjamin Aunkofer

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):

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

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

SAP.EKPO

Wir benötigen daraus nur folgende Zeilen:

1 2	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.

1 2	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):

1	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:

1 2	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.

Materialklasse

# 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()

Materialnummer

MATNR_MENGEN = EKPO_X.MENGE.groupby(EKPO_X.MATNR).sum()

MATNR_PREISEINHEIT_MENGE = EKPO_X.Preiseinheitsmenge.groupby(EKPO_X.MATNR).sum()

MATNR_VOLUMEN = EKPO_X.NETWR.groupby(EKPO_X.MATNR).sum()

MATNR_Preise = EKPO_X.NETPR.groupby(EKPO_X.MATNR).mean()MATNR_MENGEN = MATNR_MENGEN.to_frame()

MATNR_PREISEINHEIT_MENGE = MATNR_PREISEINHEIT_MENGE.to_frame()

MATNR_Preise = MATNR_Preise.to_frame()

MATNR_VOLUMEN = MATNR_VOLUMEN.to_frame()EKPO_G = MATNR_MENGEN.join(MATNR_PREISEINHEIT_MENGE, how='left')

EKPO_G = EKPO_G.join(MATNR_Preise, how='left')

EKPO_G = EKPO_G.join(MATNR_VOLUMEN, how='left')

EKPO_G = EKPO_G.sort_values(['NETWR'], ascending=False)

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:

1 2	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:

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

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):

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

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:

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

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:

1 2	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.

Benford-Analyse

September 21, 2016/in Audit Analytics, Mathematics, Python, Statistics, Tutorial/by Benjamin Aunkofer

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

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…

1 2	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):

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):

1 2	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.

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)

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.

Die Rastrigin-Funktion

September 15, 2016/in Mathematics, optimization, Python, Tutorial, Visualization/by Benjamin Aunkofer

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(2\pi 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:

1	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()

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()

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.

Einführung in WEKA

July 20, 2016/in Business Analytics, Business Intelligence, Data Mining, Data Science, Main Category, Statistics, Tool Introduction, Tutorial, Visualization/by Victoria Ivanova

Waikato Environment for Knowledge Analysis, kurz WEKA, ist ein quelloﬀenes, umfangreiches, plattformunabhängiges Data Mining Softwarepaket. WEKA ist in Java geschrieben und wurde an der WAIKATO Iniversität entwickelt. In WEKA sind viele wichtige Data Mining/Machine Learning Algorithmen implementiert und es gibt extra Pakete, wie z. B. LibSVM für Support Vector Machines, welches nicht in WEKA direkt implementiert wurde. Alle Einzelheiten zum Installieren und entsprechende Download-Links ﬁndet man unter auf der Webseite der Waikato Universität. Zusammen mit der Software wird ein Manual und ein Ordner mit Beispiel-Datensätzen ausgeliefert. WEKA arbeitet mit Datensätzen im sogenannten attribute-relation file format, abgekürzt arﬀ. Das CSV-Format wird aber ebenfalls unterstützt. Eine Datei im arﬀ-Format ist eine ASCI-Textdatei, welche aus einem Header- und einem Datateil besteht. Im Header muss der Name der Relation und der Attribute zusammen mit dem Typ stehen, der Datenteil beginnt mit einem @data-Schlüsselwort. Als Beispiel sei hier ein Datensatz mit zwei Attributen und nur zwei Instanzen gegeben.

@relation my_relation

@attribute first_attribute numeric

@attribute second_attribute numeric

@attribute class {-1,1}

@data

2.5 3.8 1

1.2 1.5 -1

WEKA unterstützt auch direktes Einlesen von Daten aus einer Datenbank (mit JDBC) oder URL. Sobald das Tool installiert und gestartet ist, landet man im Hauptmenü von WEKA – WEKA GUI Chooser 1.

Abbildung 1: WEKA GUI Chooser

Der GUI Chooser bietet den Einstieg in WEKA Interfaces Explorer, Experimenter, KnowledgeFlow und simple CLI an. Der Explorer ist ein graphisches Interface zum Bearbeiten von Datensätzen, Ausführen von Algorithmen und Visualisieren von den Resultaten. Es ist ratsam, dieses Interface als Erstes zu betrachten, wenn man in WEKA einsteigen möchte. Beispielhaft führen wir jetzt ein paar Algorithmen im Explorer durch.

Der Explorer bietet mehrere Tabs an: Preprocess, Classify, Cluster, Associate, Select attributes und Visualize. Im Preprocess Tab hat man die Möglichkeit Datensätze vorzubereiten. Hier sind zahlreiche Filter zum Präprozessieren von Datensätzen enthalten. Alle Filter sind in supervised und unsupervised unterteilt, je nachdem, ob das Klassenattribut mitbetrachtet werden soll oder nicht. Außerdem kann man entweder Attribute oder Instanzen betrachten, mit Attributen lässt man Filter spaltenweise arbeiten und bei Instanzen reihenweise. Die Auswahl der Filter ist groß, man kann den ausgewählten Datensatz diskretisieren, normalisieren, Rauschen hinzufügen etc. Unter Visualize können z. B. die geladenen Datensätze visualisert werden. Mit Select attributes kann man mithilfe von Attribut Evaluator und Search Method ein genaueres Ergebnis erzielen. Wenn man im Preprocess den Datensatz lädt, erhält man einen Überblick über den Datensatz und dessen Visualisierung. Als Beispiel wird hier der Datensatz diabetes.arﬀ genommen, welcher mit WEKA zusammen ausgeliefert wird. Dieser Datensatz enthält 768 Instanzen mit je 9 Attributen, wobei ein Attribut das Klassenattribut ist. Die Attribute enthalten z. B. Informationen über die Anzahl der Schwangerschaften, diastolischer Blutdruck, BMI usw. Alle Attribute, außer dem Klassenattribut, sind numerisch. Es gibt zwei Klassen tested negativ und tested positiv, welche das Resultat des Testens auf diabetes mellitus darstellen. über Preprocess -> Open File lädt man den Datensatz in WEKA und sieht alle relevanten Informationen wie z. B. Anzahl und Name der Attribute. Nach dem Laden kann der Datensatz klassiﬁziert werden.

Abbildung 2: Diabetes.arﬀ Datensatz geladen in WEKA

Hierzu einfach auf Classify klicken und unter Choose den gewünschten Algorithmus auswählen. Für diesen Datensatz wählen wir jetzt den Algorithmus kNN (k-Nearest Neighbour). Der Algorithmus klassiﬁziert das Testobjekt anhand der Klassenzugehörigkeit von den k Nachbarobjekten, die am nähsten zu dem Testobjekt liegen. Die Distanz zwischen den Objekten und dem Testobjekt wird mit einer Ähnlichkeitsmetrik bestimmt, meistens als euklidische oder Manhattan-Distanz. In WEKA ist der Algorithmus unter lazy iBk zu ﬁnden. Wenn man auf das Feld neben dem Algorithmusnamen in WEKA mit rechter Maustaste klickt, kann man unter show properties die Werte für den ausgewählten Algorithmus ändern, bei iBk kann man u.A. den Wert für k ändern. Für den ausgewählten Datensatz diabetes.arﬀ stellen wir beispielsweise k = 3 ein und führen die 10-fache Kreuzvalidierung durch, indem wir unter Test Options die Cross Validation auswählen. Nach der Klassiﬁkation werden die Ergebnisse in einer Warhheitsmatrix präsentiert. In unserem Fall sieht diese wie folgt aus:

a b <-- classified as

410 90 | a = tested_negative

120 148 | b = tested_positive

Die Anzahl der richtig klassiﬁzierten Instanzen beträgt 72.6563 %. Wenn man in der Result list auf den entsprechenden Algorithmus einen Rechtsklick macht, kann man z. B. noch den Fehler der Klassiﬁzierung visualisieren. Entsprechend lassen sich im Explorer unter Cluster Clustering-Algorithmen und unter Associate Assoziationsalgorithmen auf einen ausgewählten Datensatz anwenden. Die restlichen Interfaces von WEKA bieten z. T. die gleiche Funktionalität oder erweitern die Möglichkeiten des Experimentierens, fordern aber mehr Erfahrung und Wissen von dem User. Das Experimenter Interface dient dazu, mehrere Datensätze mit mehreren Algorithmen zu analysieren. Mit diesem Interface kann man groß-skalierte Experimente durchführen. Simple CLI bietet dem User eine Kommandozeile, statt einem graphischen Interface, an.

Eine Hadoop Architektur mit Enterprise Sicherheitsniveau

June 1, 2016/in Big Data, Data Science Hack, Data Security, Data Warehousing, Database, Datacenter, Hacking, Hadoop, Hadoop Framework, Main Category, Tutorial/by Georgios Gkekas

Dies ist Teil 3 von 3 der Artikelserie zum Thema Eine Hadoop-Architektur mit Enterprise Sicherheitsniveau.

Die ideale Lösung

Man denkt, dass die Integration einer sehr alten Technologie, wie ActiveDirectory oder LDAP zusammen mit einem etablierten und ausgereiften Framework wie Hadoop reibungslos funktionieren würde. Leider sind solche Annahmen in der IT Welt zu gut um wahr zu sein. Zum Glück gibt es bereits erste Erfahrungsberichte von Unternehmen, die ihre Hadoop Infrastruktur an ein zentrales IMS gekoppelt haben.

Da die meisten Unternehmen Active Directory als IMS benutzen, werden die im Folgenden dargestellte Bilder und Architekturen dies ebenfalls tun. Die vorgeschlagene Architektur ist jedoch derartig flexibel und technologieunabhängig, dass man das Active Directory auf den Bildern problemlos gegen LDAP austauschen könnte. Vielmehr ist die Integration eines Hadoop Clusters mit LDAP einfacher, da beide Technologien nativ zu Linux sind.

Schritt Eins – Integration von Hadoop mit Active Directory

Der erste Schritt, um Hadoop in dasActive Directory zu integrieren, ist ein sogenannter One-Way Trust von der Linux Welt hin zur Windows Welt . Dabei ist das Vertrauen des Authentisierungsmechanismuses von Hadoop zum Active Directory gemeint. Alle Identity Management Systeme bieten diese Funktionalität an, um sich gegenseitig vertrauen zu können und User aus anderen Domänen (Realms) zu akzeptieren. Das ermöglicht z.B. globalen Firmen mit vielen Standorten und unterschiedlichen IT Infrastrukturen und Identity Management Systemen diese zu verwalten und miteinander kommunizieren zu lassen.

Das Key Distribution Center (KDC) von Kerberos ist das Herz des Kerberos Systems im Hadoop. Hier werden die User und ihre Passwörter oder Keytabs geschützt und verwaltet. Dabei brauchen wir lediglich den One Way Trust von KDC zu Active Directory. Allerdings gibt es eine vielversprechendere Technologie, die FreeIPA. Diese hat laut Wikipedia das Ziel, ein einfach zu verwaltendes Identity,-Policy-and-Audit-System (IPA) zur Verfügung zu stellen. Seit der Version 3.0.0 kann sich FreeIPA in das Active Directory integrieren. Die aussagekräftigen Vorteile von FreeIPA sind folgende:

Reibungslose Integration mit Active Directory
Es wird zusammen mit der Technologie SSSD geliefert, die das temporäre Speichern von Rechten und Passwörtern erlaubt. Das erlaubt auch offline den Zugriff auf Fähigkeiten und Unabhängigkeit vom zentralen IPA, dem unterliegenden System.
Integrierte Kerberos und Single Sign On (SSO) Funktionalitäten.

Wir lassen dann FreeIPA die Verwaltung von Kerberos und die primäre Authentisierung unseres Clusters übernehmen. Sowohl das Active Directory, als auch FreeIPA erlauben eine kinderleichte Umsetzung des One Way Trusts mithilfe von Web Tools. Im Prinzip muss man beim One Way Trust lediglich die öffentlichen Zertifikate jedes Tools mit denen der anderen bekannt machen.

Schritt Zwei – Synchronisation der Rechte & Rollen von Active Directory

Jetzt sind alle User, die sich im Active Directory befinden, unserem Hadoop Cluster bekannt. Ein User kann sich mithilfe des kinit Kommandos und nach Eingabe seines Usernames und Passwortes einloggen. Aber man braucht auch die im Active Directory definierten Rollen und Gruppen, um eine Autorisierung mithilfe von Ranger oder Sentry zu ermöglichen. Ohne die Provisionierung der Rollen haben wir bei der Autorisierung ein ähnliches Problem, wie es bei der Authentisierung aufgetreten ist. Man müsste die Rollen selber verwalten, was nicht ideal ist.

Zum Glück gibt es verschiedene Ansätze um eine regelforme Synchronisierung der Gruppen von Active Directory in Ranger oder Sentry zu implementieren. Ranger kommt mit einem LDAP Plugin namens uxugsync, das sowohl mit LDAP als auch mit dem Active Directory kommunizieren kann. Leider hat die aktuelle Version dieses Plugins einige Nachteile:

Leistungsprobleme, weil es defaultsmäßig versucht, den ganzen Hierarchiebaum von Active Directory zu synchronisieren. Das kann zu einem großen Problem für große Firmen werden, die mehrere tausend User haben. Außerdem müssen nicht alle User Zugriff auf Hadoop haben.
Man kann bestimmte User syncen lassen, indem man ihren Gruppename im Gruppenfeld vom Plugin einträgt. Nachteil dabei ist, dass diese Abfrage nicht rekursiv funktioniert und alle Gruppe die im Ranger sein sollen einzeln abgefragt werden müssen, Das wiederum skaliert nicht sonderlich gut.
Massive und regelmäßige Abfragen des Plugins können sogar zu einem DDoS Angriff auf den zentralen Active Directory führen.

Eine bessere Lösung wäre es, wenn wir die schönen Features des SSSD Deamons (der wie oben beschrieben zusammen mit FreeIPA kommt) ausnutzen könnten. Mithilfe von SSSD werden alle User und ihre entsprechenden Gruppen dem unterliegenden Linux Betriebssystem bekannt gemacht. Das bedeutet, dass man ein einfaches Script schreiben könnte, das die User und ihre Gruppen vom System direkt abfragt und zu Ranger oder Sentry über ihre entsprechende REST APIs überträgt. Dabei schont man sowohl das Active Directory vor regelmäßigen und aufwändigen Abfragen und schafft sogar ein schnelleres Mapping der Rollen zwischen Hadoop und Betriebssystem, auch wenn Active Directory nicht erreichbar ist. Es gibt derzeit Pläne, ein solches Plugin in den nächsten Versionen von Ranger mitzuliefern.

Schritt Drei – Anlegen und Verwaltung von technischen Usern

Unser System hat jedoch neben personalisierten Usern, die echten Personen in einem Unternehmen entsprechen, auch technische User. Die technischen Users (Nicht Personalisierte Accounts – NPA), sind die Linux User mit denen die Hadoop Dienste gestartet werden. Dabei hat HDFS, Ambari usw. jeweils seinen eigenen User mit demselben Namen. Rein theoretisch könnten diese User auch im Active Directory einen Platz finden.

Meiner Meinung nach gehören diese User aber nicht dorthin. Erstens, weil sie keine echten User sind und zweitens, weil die Verwaltung solcher User nach Upgrades oder Neuinstallation des Clusters schwierig sein kann. Außerdem müssen solche User nicht den gleichen Sicherheitspolicies unterliegen, wie die normalen User. Am besten sollten sie kein Passwort besetzen, sondern lediglich ein Kerberos Keytab, das sich nach jedem Upgrade oder Neuinstallierung des Clusters neu generiert und in FreeIPA angelegt ist. Deswegen neige ich eher dazu, die NPAs in IPA anzulegen und zu verwalten.

High Level Architektur

Das folgende Bild fasst die Architektur zusammen. Hadoop Dienste, die üblicherweise in einer explorativen Umgebung benutzt werden, wie Hive und HBase, werden mit dargestellt. Es ist wichtig zu beachten, dass jegliche Technologie, die ein Ausführungsengine für YARN anbietet, wie Spark oder Storm, von dieser Architektur ebenfalls profitiert. Da solche Technologien nicht direkt mit den unterliegenden Daten interagieren, sondern diese immer über YARN und die entsprechenden Datanodes erhalten, benötigen sie auch keine besondere Darstellung oder Behandlung. Der Datenzugriff aus diesen 3rd Party Technologien respektiert die im Ranger definierten ACLs und Rollen des jeweiligen Users, der sie angestoßen hat.

Architektur in einer Mehrclusterumgebung

Wir haben schon das Argument untermauert, warum unsere technischen User direkt im IPA liegen sollten. Das kann jedoch insofern Probleme verursachen, wenn man mit mehreren Clustern arbeitet, die alle die gleichen Namen für ihre technischen User haben. Man merkt sofort, dass es sich hier um eine Namenskollision handelt. Es gibt zwei Lösungsansätze hierfür:

Man fügt den Namen Präfixen, die als kurze Beschreibungen der jeweiligen Umgebung dienen, wie z.B. ada, proj1, proj2 hinzu. Dadurch haben die User unterschiedliche Namen, wie proj1_hdfs für die proj1 Umgebung und ada_hdfs für die ada Umgebung. Man kann diese Lösung auch bei Kerberos KDCs benutzen, die in jeder Umgebung dediziert sind und die technischen User der jeweiligen Umgebung beibehalten.
Man benutzt einen separaten Realm für jede Umgebung und damit auch eine separate IPA Instanz. Hier gibt es wiederum zwei verschiedene Ansätze. Ich muss jedoch zugeben, dass ich die Zweite nie ausprobiert habe und daher für ihre Durchführbarkeit nicht garantieren kann:
1. Man bindet jede Umgebung einzeln über ihre FreeIPA per One Way Trust an das zentrale Active Directory. Das hat natürlich den Nachteil einer uneinheitlichen User Management Infrastruktur für alle Umgebungen, da Jede ihre eigene IPA Infrastruktur verwaltet und wartet.
2. Man baut einen Hierarchiebaum von unterschiedlichen IPA Instanzen, so wie man es bei Forests von Active Directory Instanzen macht.

Das folgende Bild stellt den letzten Ansatz dar. Im Prinzip haben wir hier einen hierarchischen IPA Cluster mit mehreren One Way Trusts von den lokalen IPA Instanzen zu der zentralen IPA.

Zusammenfassung

Wie Sie vielleicht von der gesamten Diskussion her abgeleitet haben, ist die Umsetzung einer unternehmerisch-konformen und personenbasierten Sicherheitsarchitektur innerhalb von Hadoop keine einfache Sache. Man muss mit unterschiedlichen Architekturen und Ansätzen spielen, bevor man einen relativ vernünftigen oder sogar idealen Zustand erreicht hat. Die Berücksichtigung der jeweiligen IT Architektur spielt dabei eine sehr große Rolle. Ich hoffe, ich konnte die wichtigsten Merkmalen einer solchen Architektur und die Punkte, die ein Architekt besonders beachten muss, klar darstellen.

Als Zusammenfassung habe ich Ihnen am Ende eine Art Shoppingliste aller Komponenten zusammengestellt, die wichtig für den personalisierten Zugriff im Hadoop sind:

Kerberos – Authentisierung
FreeIPA – Authentisierung, Integration mit Active Directory
Active Directory oder LDAP
Ranger oder Sentry
1. Plugin für Rollen/Gruppen Mapping zwischen AD und dem Betriebssystem
Optional SSSD für schnellere Abfrage der Gruppen und Rollen des Betriebssystems

Zurück zu Teil 2 von 3 – Sicherheitstechnologie in Hadoop

Eine Hadoop Architektur mit Enterprise Sicherheitsniveau

May 20, 2016/in Big Data, Business Analytics, Business Intelligence, Cloud, Data Science Hack, Data Security, Data Warehousing, Database, Datacenter, Hacking, Hadoop, Hadoop Framework, Main Category, NoSQL, Tools, Tutorial/by Georgios Gkekas

Dies ist Teil 2 von 3 der Artikelserie zum Thema Eine Hadoop-Architektur mit Enterprise Sicherheitsniveau.

Der aktuelle Stand der Technologie

Zum Glück ist Hadoop heutzutage ein bisschen reifer, als es noch vor zehn Jahren war. Es gibt viele Tools, einige davon OpenSource und einige lizenziert, die den Sicherheitsmangel im Hadoop zu lösen versuchen. Die Tabelle unten zeigt eine Auswahl der am meisten genutzten Sicherheitstools. Da jedes Tool von einer anderen Hadoop Distribution bevorzugt wird, habe ich diese Parameter mit berücksichtigt.

Es ist zu beachten, dass die zwei populärsten Hadoop Distributions (Hortonworks und Cloudera) kaum Unterschiede aufweisen, wenn man sie auf funktionaler Ebene vergleicht. Der größte Unterschied besteht darin, dass Hortonworks ein Open Source und Cloudera ein kommerzielles Produkt ist. Abgesehen davon hat jeder Vendor den einen oder anderen Vorteil, ein ausführlicher Vergleich würde jedoch den Rahmen dieses Artikels sprengen.

Hadoop kommt von der Stange ohne aktivierte Authentisierung. Die Hadoop Dienste vertrauen jedem User, egal als was er oder sie sich ausgibt. Das sieht folgendermaßen aus:

Angenommen Mike arbeitet an einer Maschine, die ihm Zugriff auf den Hadoop Cluster erlaubt und Sudo-Rechte gibt. Aber Mike hat das Passwort für den hdfs Superuser nicht. Er kann sich jetzt einfach als der hdfs User ausgeben, indem er die folgenden Kommandos ausführt. Dabei bekommt er fatalerweise alle Rechten des hdfs Superusers und ist in der Lage das gesamte HDFS Filesystem zu löschen. Es würde sogar bereits der Environment variabel USER ausreichen, um einen anderen User umzuwandeln.

Kerberos ist im Moment der einzige Weg um Authentisierung im Hadoop zu gewährleisten. Kein Weg führt daran vorbei, es sei denn, man ist verrückt genug, um ein hochkompliziertes System auf Linux basierter ACLs auf jeder Maschine zu installieren und zu verwalten, um User daran zu hindern sich falsch zu authentifizieren. Es ist zudem wichtig zu beachten, dass Kerberos als einziges Sicherheitsmerkmal zur Authentifizierung dient, aber ohne richtige Authentisierung gibt es auch keine richtige Autorisierung. Wenn User jetzt selbst in der Lage sind, sich beliebig als jemand anderes auszugeben, können sie so selbst zu den sensibelsten Daten unbefugten Zugriff erlangen.

Apache Ranger oder Sentry erlauben die Definition und Verwaltung von Access Control Lists (ACLs). Diese Listen legen fest, welche User Zugriff auf welchen Bereich des HDFS Filesystems haben Der gleiche Effekt kann auch ohne diese Tools, durch einfache Hadoop ACLs erreicht werden, die den normalen Linux ACLs ähneln. Es empfiehlt sich jedoch die neuesten Tools zu benutzen, wegen a) ihrer Benutzerfreundlichkeit, b) ihrer ausgearbeiteten APIs, die einem Administrator erlauben die Listen ohne GUI zu verwalten und beim Programmieren sogar zu automatisieren, und c) wegen ihrer Auditingfähigkeiten, die das Nachverfolgen von Zugriffen und Aktionen ermöglichen.

Anbei ist das Bild einer Ranger Policy, die der Gruppe der User rekursiv Lese- und Ausführungsrechte auf das Verzeichnis /projects/autonomous_driving gibt.

Alle einzelne Stücke des Puzzles kommen zusammen

Nachdem wir ermittelt haben, welche Technologien es gibt, die uns zu einem sicheren Cluster verhelfen, müssen diese im nächsten Schritt zusammengesetzt werden. Zum Glück hat jeder Vendor seine eigene Technologie, um Tools aus dem Hadoop Ecosystem zu integrieren und zu verwalten. Cloudera beispielsweise bietet den sehr wirksamen Cloudera Manager und Hortonworks das Apache Ambari an. Die beiden Tools kümmern sich um das Anlegung der technischen Hadoop User (hdfs, hadoop, hive, ranger, e.t.c.) und der entsprechenden Kerberos Keytabs, die den technischen Usern erlauben, sich gegenüber Hadoop zu authentisieren. Am Ende der Installation hat man sämtliche Konfigurationen zentral platziert und kann neue personalisierte Accounts anlegen. Man kann sich dann im Ranger oder Sentry Web UI anmelden und ACLs für die User und Gruppen definieren.

Das ist allerdings nicht der Idealzustand. Jedes Unternehmen verwaltet ihre User bereits in bestimmten Verwaltungssystemen, die sich innerhalb der IT Infrastruktur befinden. Diese Systeme (oder auch Identity Management Systems) sind ein wichtiges vertikales, abteilungsübergreifendes Element der unternehmerischen IT Architektur. Jedes EDS Tool im Unternehmen ist an ein Identity Management System, wie Active Directory oder LDAP, gekoppelt und muss damit die User nicht selbst verwalten.

Der Stellenwert solcher Tools wird sofort erkennbar, wenn man die strengen Sicherheitsregeln eines modernen Unternehmens betrachtet: Passwörter müssen bestimmte Kriterien erfüllen und alle 30 Tagen gewechselt werden. Außerdem darf niemand eins seiner letzten zehn Passwörter benutzen.

Eine IT Architektur, die die Implementierung solcher unternehmensbreiten Anforderungen in jeder einzelne Applikation fördert ist der Alptraum jedes Applikationsentwicklers und zeigt das Versagen des IT-Architekten.

Aber lassen Sie uns zurück zu unserem Hauptthema kommen. Wie können wir ein System wie Active Directory oder LDAP in Hadoop integrieren? Der nächste Abschnitt gibt die Antwort auf diese Frage.

Weiter zu Teil 3 von 3 – Eine Einterprise Hadoop Architektur für beste Sicherheit

Zurück zu Teil 1 von 3 – Motivation und Anforderungen einer Data Science Plattform

Machine Learning mit Python – Minimalbeispiel

April 26, 2016/in Artificial Intelligence, Data Mining, Data Science Hack, Machine Learning, Mathematics, optimization, Predictive Analytics, Python, Tutorial, Visualization/by Benjamin Aunkofer

Maschinelles Lernen (Machine Learning) ist eine Gebiet der Künstlichen Intelligenz (KI, bzw. AI von Artificial Intelligence) und der größte Innovations- und Technologietreiber dieser Jahre. In allen Trendthemen – wie etwa Industrie 4.0 oder das vernetzte und selbstfahrende Auto – spielt die KI eine übergeordnete Rolle. Beispielsweise werden in Unternehmen viele Prozesse automatisiert und auch Entscheidungen auf operativer Ebene von einer KI getroffen, zum Beispiel in der Disposition (automatisierte Warenbestellungen) oder beim Festsetzen von Verkaufspreisen.

Aufsehen erregte Google mit seiner KI namens AlphaGo, einem Algortihmus, der den Weltmeister im Go-Spiel in vier von fünf Spielen besiegt hatte. Das Spiel Go entstand vor mehr als 2.500 Jahren in China und ist auch heute noch in China und anderen asiatischen Ländern ein alltägliches Gesellschaftsspiel. Es wird teilweise mit dem westlichen Schach verglichen, ist jedoch einfacher und komplexer zugleich (warum? das wird im Google Blog erläutert). Machine Learning kann mit einer Vielzahl von Methoden umgesetzt werden, werden diese Methoden sinnvoll miteinander kombiniert, können durchaus äußerst komplexe KIs erreicht werden. Der aktuell noch gängigste Anwendungsfall für Machine Learning ist im eCommerce zu finden und den meisten Menschen als die Produktvorschläge von Amazon.com bekannt: Empfehlungsdienste (Recommender System).

Klassifikation via K-Nearest Neighbour Algorithmus

Ein häufiger Zweck des maschinellen Lernens ist, technisch gesehen, die Klassifikation von Daten in Abhängigkeit von anderen Daten. Es gibt mehrere ML-Algorithmen, die eine Klassifikation ermöglichen, die wohl bekannteste Methode ist der k-Nearest-Neighbor-Algorithmus (Deutsch:„k-nächste-Nachbarn”), häufig mit “kNN” abgekürzt. Das von mir interviewte FinTech StartUp Number26 nutzt diese Methodik beispielsweise zur Klassifizierung von Finanztransaktionen.

Um den Algorithmus Schritt für Schritt aufbauen zu können, müssen wir uns

Natürlich gibt es in Python, R und anderen Programmiersprachen bereits fertige Bibliotheken, die kNN bereits anbieten, denen quasi nur Matrizen übergeben werden müssen. Am bekanntesten ist wohl die scikit-learn Bibliothek für Python, die mehrere Nächste-Nachbarn-Modelle umfasst. Mit diesem Minimalbeispiel wollen wir den grundlegenden Algorithmus von Grund auf erlernen. Wir wollen also nicht nur machen, sondern auch verstehen.

Vorab: Verwendete Bibliotheken

Um den nachstehenden Python-Code (Python 3.x, sollte allerdings auch mit Python 2.7 problemlos funktionieren) ausführen zu können, müssen folgende Bibliotheken eingebunden werden:

import numpy as numpy

import matplotlib.pyplot as pyplot

from mpl_toolkits.mplot3d import Axes3D #Erweiterung für die Matplotlib - siehe: http://matplotlib.org/mpl_toolkits/

Übrigens: Eine Auflistung der wohl wichtigsten Pyhton-Bibliotheken für Datenanalyse und Datenvisualisierung schrieb ich bereits hier.

Schritt 1 – Daten betrachten und Merkmale erkennen

Der erste Schritt ist tatsächlich der aller wichtigste, denn erst wenn der Data Scientist verstanden hat, mit welchen Daten er es zu tun hat, kann er die richtigen Entscheidungen treffen, wie ein Algorithmus richtig abgestimmt werden kann und ob er für diese Daten überhaupt der richtige ist.

In der Realität haben wir es oft mit vielen verteilten Daten zu tun, in diesem Minimalbeispiel haben wir es deutlich einfacher: Der Beispiel-Datensatz enthält Informationen über Immobilien über vier Spalten.

Quadratmeter: Größe der nutzbaren Fläche der Immobilie in der Einheit m²
Wandhoehe: Höhe zwischen Fußboden und Decke innerhalb der Immobilie in der Einheit m
IA_Ratio: Verhältnis zwischen Innen- und Außenflächen (z. B. Balkon, Garten)
Kategorie: Enthält eine Klassifizierung der Immobilie als “Haus”, “Wohnung” und “Büro”

[box]Hinweis für Python-Einsteiger: Die Numpy-Matrix ist speziell für Matrizen-Kalkulationen entwickelt. Kopfzeilen oder das Speichern von String-Werten sind für diese Datenstruktur nicht vorgesehen![/box]

def readDataSet(filename):

fr = open(filename) # Datei-Stream vorbereiten

numberOfLines = len(fr.readlines()) # Anzahl der Zeilen ermitteln

returnMat = numpy.zeros((numberOfLines-1,3)) # Eine Numpy-Matrix in Höhe der Zeilenanzahl (minus Kopfzeile) und in Breite der drei Merkmal-Spalten

classLabelVector = [] # Hier werden die tatsächlichen Kategorien (Haus, Wohnung, Büro) vermerkt

classColorVector = [] # Hier werden die Kategorien über Farben vermerkt (zur späteren Unterscheidung im 3D-Plot!)

#print(returnMat) # Ggf. mal die noch die ausge-null-te Matrix anzeigen lassen (bei Python 2.7: die Klammern weglassen!)

fr = open(filename) # Datei-Stream öffnen

index = 0

for line in fr.readlines(): # Zeile für Zeile der Datei lesen

if index != 0: # Kopfzeile überspringen

line = line.strip()

listFromLine = line.split('\t') # Jede Zeile wird zur temporären Liste (Tabulator als Trennzeichen)

returnMat[index-1,:] = listFromLine[1:4] #Liste in die entsprechende Zeile der Matrix überführen

classLabel = listFromLine[4] # Kategorie (Haus, Wohnung, Büro) für diese Zeile merken

if classLabel == "Buero":

color = 'yellow'

elif classLabel == "Wohnung":

color = 'red'

else:

color = 'blue'

classLabelVector.append(classLabel) # Kategorie (Haus, Wohnung, Büro) als Text-Label speichern

classColorVector.append(color) # Kategorie als Farbe speichern (Büro = gelb, Wohnung = rot, Haus = Blau)

index += 1

return returnMat,classLabelVector, classColorVector

Aufgerufen wird diese Funktion dann so:

1 2	dataSet, classLabelVector, classColorVector = readDataSet("K-Nearst_Neighbour-DataSet.txt")

Die Matrix mit den drei Spalten (Quadratmeter, Wandhohe, IA_Ratio) landen in der Variable “dataSet”.

Schritt 2 – Merkmale im Verhältnis zueinander perspektivisch betrachten

Für diesen Anwendungsfall soll eine Klassifizierung (und gewissermaßen die Vorhersage) erfolgen, zu welcher Immobilien-Kategorie ein einzelner Datensatz gehört. Im Beispieldatensatz befinden sich vier Merkmale: drei Metriken und eine Kategorie (Wohnung, Büro oder Haus). Es stellt sich zunächst die Frage, wie diese Merkmale zueinander stehen. Gute Ideen der Datenvisualisierung helfen hier fast immer weiter. Die gängigsten 2D-Visualisierungen in Python wurden von mir bereits hier zusammengefasst.

[box]Hinweis: In der Praxis sind es selten nur drei Dimensionen, mit denen Machine Learning betrieben wird. Das Feature-Engineering, also die Suche nach den richtigen Features in verteilten Datenquellen, macht einen wesentlichen Teil der Arbeit eines Data Scientists aus – wie auch beispielsweise Chief Data Scientist Klaas Bollhoefer (siehe Interview) bestätigt.[/box]

fig = pyplot.figure()

ax = fig.add_subplot(111)

ax.scatter(dataSet[:,0], dataSet[:,1], marker='o', color=classColorVector)

ax.set_xlabel("Raumflaeche in Quadratmeter")

ax.set_ylabel("Wandhohe")

ax.set_xlim(xmin=0)

ax.set_ylim(ymin=0)

pyplot.show()

fig = pyplot.figure()

ax = fig.add_subplot(111)

ax.scatter(dataSet[:,0], dataSet[:,2], marker='o', color=classColorVector)

ax.set_xlabel("Raumflaeche in Quadratmeter")

ax.set_ylabel("IA_Ratio")

ax.set_xlim(xmin=0)

ax.set_ylim(ymin=0)

pyplot.show()

Die beiden Scatter-Plots zeigen, das Häuser (blau) in allen Dimensionen die größte Varianz haben. Büros (gelb) können größer und höher ausfallen, als Wohnungen (rot), haben dafür jedoch tendenziell ein kleineres IA_Ratio. Könnten die Kategorien (blau, gelb, rot) durch das Verhältnis innerhalb von einem der beiden Dimensionspaaren in dem zwei dimensionalen Raum exakt voneinander abgegrenzt werden, könnten wir hier stoppen und bräuchten auch keinen kNN-Algorithmus mehr. Da wir jedoch einen großen Überschneidungsbereich in beiden Dimensionspaaren haben (und auch Wandfläche zu IA_Ratio sieht nicht besser aus),

Eine 3D-Visualisierung eignet sich besonders gut, einen Überblick über die Verhältnisse zwischen den drei Metriken zu erhalten: (die Werte wurden hier bereits normalisiert, liegen also zwischen 0,00 und 1,00)

fig = pyplot.figure()

ax = fig.add_subplot(111, projection='3d')

ax.scatter(dataSet[:,0], dataSet[:,2], dataSet[:,1], marker='o', color=classColorVector)

ax.set_xlabel("Raumflaeche in Quadratmeter")

ax.set_ylabel("IA_Ratio")

ax.set_zlabel("Wandhoehe in Meter")

ax.set_xlim(xmin=0)

ax.set_ylim(ymin=0)

ax.set_zlim(zmin=0)

pyplot.show()

Es zeigt sich gerade in der 3D-Ansicht recht deutlich, dass sich Büros und Wohnungen zum nicht unwesentlichen Teil überschneiden und hier jeder Algorithmus mit der Klassifikation in Probleme geraten wird, wenn uns wirklich nur diese drei Dimensionen zur Verfügung stehen.

Schritt 3 – Kalkulation der Distanzen zwischen den einzelnen Punkten

Bei der Berechnung der Distanz in einem Raum hilft uns der Satz des Pythagoras weiter. Die zu überbrückende Distanz, um von A nach B zu gelangen, lässt sich einfach berechnen, wenn man entlang der Raumdimensionen Katheten aufspannt.

$c = \sqrt{a^2+ b^2}$

Die Hypotenuse im Raum stellt die Distanz dar und berechnet sich aus der Wurzel aus der Summe der beiden Katheten im Quadrat. Die beiden Katheten bilden sich aus der Differenz der Punktwerte (q, p) in ihrer jeweiligen Dimension.Bei mehreren Dimensionen gilt der Satz entsprechend:

$Distanz = \sqrt{(q_1-p_1)^2+(q_2-p_2)^2+…+(q_n-p_n)^2}$

Um mit den unterschiedlichen Werte besser in ihrer Relation zu sehen, sollten sie einer Normalisierung unterzogen werden. Dabei werden alle Werte einer Dimension einem Bereich zwischen 0.00 und 1.00 zugeordnet, wobei 0.00 stets das Minimum und 1.00 das Maximum darstellt.

$NormWert = \frac{Wert - Min}{Wertspanne} = \frac{Wert - Min}{Max - Min}$

def normalizeDataSet(dataSet):

dataSet_n = numpy.zeros(numpy.shape(dataSet)) #[[ 0. 0. 0.]

# [ 0. 0. 0.]

# ...,

# [ 0. 0. 0.]

# [ 0. 0. 0.]]

minValues = dataSet.min(0) # [ 10. 2.6 0.]

ranges = dataSet.max(0) - dataSet.min(0) # [ 1775. 2.4 68.]

minValues = dataSet.min(0) # [ 10. 2.6 0.]

maxValues = dataSet.max(0) # [ 1785. 5. 68.]

ranges = maxValues - minValues # [ 1775. 2.4 68.]

rowCount = dataSet.shape[0] # 1039

# numpy.tile() wiederholt Sequenzen (hier: [[ 10. 2.6 0. ], ..., [ 10. 2.6 0. ]]

dataSet_n = dataSet - numpy.tile(minValues, (rowCount, 1)) #[[ 2.56000000e+02 9.00000000e-01 1.80000000e+01]

# [ 6.60000000e+01 2.00000000e-01 5.40000000e+01]

# [ 3.32000000e+02 1.50000000e-01 1.00000000e+01]

# ...,

# [ 1.58000000e+02 6.00000000e-01 0.00000000e+00]

# [ 5.70000000e+01 1.00000000e-01 5.20000000e+01]

# [ 1.68000000e+02 2.00000000e-01 0.00000000e+00]]

dataSet_n = dataSet_n / numpy.tile(ranges, (rowCount, 1)) #[[ 0.14422535 0.375 0.26470588]

# [ 0.0371831 0.08333333 0.79411765]

# [ 0.18704225 0.0625 0.14705882]

# ...,

# [ 0.08901408 0.25 0.]

# [ 0.03211268 0.04166667 0.76470588]

# [ 0.09464789 0.08333333 0.]]

#print(dataSet_n)

return dataSet_n, ranges, minValues

Die Funktion kann folgendermaßen aufgerufen werden:

1 2	dataSet_n, ranges, minValues = normalizeDataSet(dataSet)

Schritt 4 & 5 – Klassifikation durch Eingrenzung auf k-nächste Nachbarn

Die Klassifikation erfolgt durch die Kalkulation entsprechend der zuvor beschriebenen Formel für die Distanzen in einem mehrdimensionalen Raum, durch Eingrenzung über die Anzahl an k Nachbarn und Sortierung über die berechneten Distanzen.

def classify(inX, dataSet, labels, k):

rowCount = dataSet.shape[0] # Anzahl an Zeilen bestimmen

diffMat = numpy.tile(inX, (rowCount,1)) - dataSet # Berechnung der Katheten

# (über tile() wird der Eingangsdatensatz über die Zeilenanzahl des dataSet vervielfacht,

# der dataSet davon substrahiert)

sqDiffMat = diffMat**2 # Quadrat der Katheten

sqDistances = sqDiffMat.sum(axis=1) # Aufsummieren der Differenzpaare

distances = sqDistances**0.5 # Quadratwurzel über alle Werte

sortedDistIndicies = distances.argsort() # Aufsteigende Sortierung

classCount = {}

#print("inX = %s, k = %s" % (inX, k))

#print(sortedDistIndicies)

for i in range(k): # Eingrenzung auf k-Werte in der sortierten Liste

closest = labels[sortedDistIndicies[i]] # Label (Kategorie [Büro, Wohnung, Haus] entsprechend der Sortierung aufnehmen

classCount[closest] = classCount.get(closest, 0) + 1 # Aufbau eines Dictionary über die

sortedClassCount = sorted(classCount, key = classCount.get, reverse=True) # Absteigende Sortierung der gesammelten Labels in k-Reichweite

# wobei die Sortierung über den Count (Value) erfolgt

#print(classCount)

#print(sortedClassCount[0])

return sortedClassCount[0] # Liefere das erste Label zurück

# also das Label mit der höchsten Anzahl innerhalb der k-Reichweite

Über folgenden Code rufen wir die Klassifikations-Funktion auf und legen die k-Eingrenzung fest, nebenbei werden Fehler gezählt und ausgewertet. Hier werden der Reihe nach die ersten 30 Zeilen verarbeitet:

errorCount = 0

k = 5 # k-Eingrenzung (hier: auf 5 Nachbarn einschränken)

rowCount = dataSet_n.shape[0] # Anzahl der Zeilen im gesamten Datensatz

numTestVectors = 30 # Datensätze 0 - 29 werden zum testen von k verwendet,

# die Datensätze ab Zeile 30 werden zur Klassifikation verwendet

for i in range(0, numTestVectors): # Aufruf des Klassifikators von 0 bis 29

result = classify(dataSet_n[i,:], dataSet_n[numTestVectors:rowCount,:], classLabelVector[numTestVectors:rowCount], k)

print("%s - the classifier came back with: %s, the real answer is: %s" %(i, result, classLabelVector[i]))

if (result != classLabelVector[i]):

errorCount += 1.0

print("Error Count: %d" % errorCount)

Nur 30 Testdatensätze auszuwählen ist eigentlich viel zu knapp bemessen und hier nur der Übersichtlichkeit geschuldet. Besser ist für dieses Beispiel die Auswahl von 100 bis 300 Datensätzen. Die Ergebnisse sind aber bereits recht ordentlich, allerdings fällt dem Algorithmus – wie erwartet – noch die Unterscheidung zwischen Wohnungen und Büros recht schwer.

0 – klassifiziert wurde: Buero, richtige Antwort: Buero
1 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
2 – klassifiziert wurde: Buero, richtige Antwort: Buero
3 – klassifiziert wurde: Buero, richtige Antwort: Buero
4 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
5 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
6 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
7 – klassifiziert wurde: Wohnung, richtige Antwort: Buero
8 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
9 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
10 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
11 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
12 – klassifiziert wurde: Buero, richtige Antwort: Buero
13 – klassifiziert wurde: Wohnung, richtige Antwort: Buero
14 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
15 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
16 – klassifiziert wurde: Buero, richtige Antwort: Buero
17 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
18 – klassifiziert wurde: Haus, richtige Antwort: Haus
19 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
20 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
21 – klassifiziert wurde: Buero, richtige Antwort: Buero
22 – klassifiziert wurde: Buero, richtige Antwort: Buero
23 – klassifiziert wurde: Buero, richtige Antwort: Buero
24 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
25 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
26 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
27 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
28 – klassifiziert wurde: Wohnung, richtige Antwort: Wohnung
29 – klassifiziert wurde: Buero, richtige Antwort: Buero
Error Count: 2

Über weitere Tests wird deutlich, dass k nicht zu niedrig und auch nicht zu hoch gesetzt werden darf.

Datensätze	k	Fehler
150	1	25
150	3	23
150	5	21
150	20	26

Ein nächster Schritt wäre die Entwicklung eines Trainingprogramms, dass die optimale Konfiguration (k-Eingrenzung, Gewichtung usw.) ermittelt.

Fehlerraten herabsenken

Die Fehlerquote ist im Grunde niemals ganz auf Null herabsenkbar, sonst haben wir kein maschinelles Lernen mehr, sondern könnten auch feste Regeln ausmachen, die wir nur noch einprogrammieren (hard-coding) müssten. Wer lernt, macht auch Fehler! Dennoch ist eine Fehlerquote von 10% einfach zu viel für die meisten Anwendungsfälle. Was kann man hier tun?

Den Algorithmus verbessern (z. B. optimale k-Konfiguration und Gewichtung finden)
mehr Merkmale finden (= mehr Dimensionen)
mehr Daten hinzuziehen (gut möglich, dass alleine dadurch z. B. Wohnungen und Büros besser unterscheidbar werden)
einen anderen Algorithmus probieren (kNN ist längst nicht für alle Anwendungen ideal!)

Das Problem mit den Dimensionen

Theoretisch kann kNN mit undenklich vielen Dimensionen arbeiten, allerdings steigt der Rechenaufwand damit auch ins unermessliche. Der k-nächste-Nachbar-Algorithmus ist auf viele Daten und Dimensionen angewendet recht rechenintensiv.

In der Praxis hat nicht jedes Merkmal die gleiche Tragweite in ihrer Bedeutung für die Klassifikation und mit jeder weiteren Dimension steigt auch die Fehleranfälligkeit, insbesondere durch Datenfehler (Rauschen). Dies kann man sich bei wenigen Dimensionen noch leicht bildlich vorstellen, denn beispielsweise könnten zwei Punkte in zwei Dimensionen nahe beieinander liegen, in der dritten Dimension jedoch weit auseinander, was im Ergebnis dann eine lange Distanz verursacht. Wenn wir beispielsweise 101 Dimensionen berücksichtigen, könnten auch hier zwei Punkte in 100 Dimensionen eng beieinander liegen, läge jedoch in der 101. Dimension (vielleicht auch auf Grund eines Datenfehlers) eine lange Distanz vor, wäre die Gesamtdistanz groß. Mit Gewichtungen könnten jedoch als wichtiger einzustufenden Dimensionen bevorzugt werden und als unsicher geltende Dimensionen entsprechend entschärft werden.

Je mehr Dimensionen berücksichtigt werden sollen, desto mehr Raum steht zur Verfügung, so dass um wenige Datenpunkte viel Leerraum existiert, der dem Algorithmus nicht weiterhilft. Je mehr Dimensionen berücksichtigt werden, desto mehr Daten müssen zu Verfügung gestellt werden, im exponentiellen Anstieg – Wo wir wieder beim Thema Rechenleistung sind, die ebenfalls exponentiell ansteigen muss.

Weiterführende Literatur

Machine Learning in Action

Introduction to Machine Learning with Python

Einführung in Data Science: Grundprinzipien der Datenanalyse mit Python

KNN: Rückwärtspass

April 12, 2016/in Artificial Intelligence, Big Data, Data Mining, Data Science, Data Science Hack, Machine Learning, Mathematics, Octave, Tutorial/by Nico Hezel

Im letzten Artikel der Serie haben wir gesehen wie bereits trainierte Netzwerke verwendet werden können. Als Training wird der Prozess bezeichnet der die Gewichte in einen Netzwerk so anpasst, dass bei einem Vorwärtspass durch ein Netzwerk zu einen festgelegten Eingangsdatensatz ein bestimmtes Ergebnis in der Ausgangsschicht ausgegeben wird. Im Umkehrschluss heißt das auch, dass wenn etwas anderes ausgeliefert wurde als erwartet, das Netzwerk entweder noch nicht gut genug oder aber auf ein anderes Problem hin trainiert wurde.

Training

Das Training selbst findet in drei Schritten statt. Zunächst werden die Gewichte initialisiert. Üblicherweise geschieht das mit zufälligen Werten, die aus einer Normalverteilung gezogen werden. Je nachdem wie viele Gewichte eine Schicht hat, ist es sinnvoll die Verteilung über den Sigma Term zu skalieren. Als Daumenregeln kann dabei eins durch die Anzahl der Gewichte in einer Schicht verwendet werden.

Im zweiten Schritt wird der Vorwärtspass für die Trainingsdaten errechnet. Das Ergebnis wird beim ersten Durchlauf alles andere als zufrieden stellend sein, es dient aber dem Rückwärtspass als Basis für dessen Berechnungen und Gewichtsänderungen. Außerdem kann der Fehler zwischen der aktuellen Vorhersage und dem gewünschten Ergebnis ermittelt werden, um zu entscheiden, ob weiter trainiert werden soll.

Der eigentliche Rückwärtspass errechnet aus der Differenz der Vorwärtspassdaten und der Zieldaten die Steigung für jedes Gewicht aus, in dessen Richtung dieses geändert werden muss, damit das Netzwerk bessere Vorhersagen trifft. Das klingt zunächst recht abstrakt, die genauere Mathematik dahinter werde ich in einem eigenen Artikel erläutern. Zur besseren Vorstellung betrachten wir die folgende Abbildung.

visuelle Darstellung aller Gewichtskombinationen und deren Vorhersagefehler

Das Diagramm zeigt in blau zu allen möglichen Gewichtskombinationen eines bestimmten, uns unbekannten, Netzwerks und Problems den entsprechenden Vorhersagefehler. Die Anzahl der Kombinationen hängt von der Anzahl der Gewichte und der Auflösung des Wertebereiches für diese ab. Theoretisch ist die Menge also unendlich, weshalb die blaue Kurve eine von mir ausgedachte Darstellung aller Kombinationen ist. Der erste Vorwärtspass liefert uns eine Vorhersage die eine normalisierte Differenz von 0.6 zu unserem eigentlichen Wunschergebnis aufweist. Visualisiert ist das Ganze mit einer schwarzen Raute. Der Rückwärtspass berechnet aus der Differenz und den Daten vom Vorwärtspass einen Änderungswunsch für jedes Gewicht aus. Da die Änderungen unabhängig von den anderen Gewichten ermittelt wurden, ist nicht bekannt was passieren würde wenn alle Gewichte sich auf einmal ändern würden. Aus diesem Grund werden die Änderungswünsche mit einer Lernrate abgeschwächt. Im Endeffekt ändert sich jedes Gewicht ein wenig in die Richtung, die es für richtig erachtet. In der Hoffnung einer Steigerung entlang zu einem lokalen Minimum zu folgen, werden die letzten beiden Schritte (Vor- und Rückwärtspass) mehrfach wiederholt. In dem obigen Diagramm würde die schwarze Raute der roten Steigung folgen und sich bei jeder Iteration langsam auf das linke lokale Minimum hinzubewegen.

Anwendungsbeispiel und Programmcode

Um den ganzen Trainingsprozess im Einsatz zu sehen, verwenden wir das Beispiel aus dem Artikel “KNN: Vorwärtspass”. Die verwendeten Daten kommen aus der Wahrheitstabelle eines X-OR Logikgatters und werden in ein 2-schichtiges Feedforward Netzwerk gespeist.

XOR Wahrheitstabelle

X1	X2	Y = X1 ⊻ X2
0	0	0
0	1	1
1	0	1
1	1	0

Der Programmcode ist in Octave geschrieben und kann zu Testzwecken auf der Webseite von Tutorialpoint ausgeführt werden. Die erste Hälfte von dem Algorithmus kennen wir bereits, der Vollständigkeit halber poste ich ihn noch einmal, zusammen mit den Rückwärtspass. Hinzugekommen sind außerdem ein paar Konsolenausgaben, eine Lernrate- und eine Iterations-Variable die angibt wie viele Trainingswiederholungen durchlaufen werden sollen.

%--------------------- Daten -----------------------

X = [0 0; % Eingangsdaten

0 1;

1 0;

1 1]

Y = [0;1;1;0] % erwartete XOR Ausgangsdaten

theta1 = normrnd(0, 1/(3*2), 3, 2); % 3x2 Gewichtsmatrix

theta2 = normrnd(0, 1/(3*1), 3, 1); % 3x1 Gewichtsmatrix

m = length(X) % Anzahl der Eingangsdaten

iteration = 10000 % Anzahl der Trainingsiterationen

alpha = 0.8 % lernrate

printf("\n\nStarte Training ... ")

for(i = 1:iteration)

%--------------------- Vorwärtspass -----------------------

V = X; % anlegen der Eingangsdaten an die Eingangsschicht

% 1. berechne die Aktivierungen der verborgenen Schicht

Vb = [ones(m,1) V]; % hinzufügen der Bias Units (sind immer 1)

Zv = Vb * theta1; % Summe aus den Eingangswerten multipliziert mit deren Gewichten

H = 1 ./ (1 .+ e.^-Zv); % anwenden der Sigmoid Funktion auf die Aktivierungsstärke Zv

% 2. berechne die Aktivierungen der Ausgangsschicht

Hb = [ones(m,1) H]; % hinzufügen der Bias Units an die verborgene Schicht

Zh = Hb * theta2; % Produkt aus den Aktivierungen der Neuronen in H und Theta2

O = 1 ./ (1 .+ e.^-Zh); % Vorhersage von dem Netzwerk

% 3. berechne die Vorhersageungenauigkeit

loss = (O .- Y) .^ 2; % quadratischer Fehler von der Vorhersage und der Zielvorgabe Y

mse = sum(loss) / m; % durchschnittlicher quadratischer Fehler aller Vorhersagen

%--------------------- Rückwärtspass -----------------------

% 1. Ableitung der Fehlerfunktion

d = O .- Y; % Differenzmatrix zwischen der Vorhersage und der Zielvorgabe Y

% 2. berechne die Änderungen für Theta2 und die Ableitung der Ausgangsschicht

OMO = ones(size(O)) .- O; % Zwischenvariable: 1-Minus-Vorhersage

Zhd = d .* O .* OMO; % Ableitung der Sigmoid Funktion

theta2c = Hb' * Zhd; % Änderunswunsch für Theta2

Hd = Zhd * theta2'; % Ableitung von der Ausgangsschicht

Hd(:,[1]) = []; % Ableitung von der Bias Unit

% 3. berechne die Änderungen für Theta1 und die Ableitung der verborgenen Schicht

HMO = ones(size(H)) .- H; % Zweischenvariable: 1 Minus Aktivierung der verborgenen Schicht

Zvd = Hd .* H .* HMO; % Ableitung der Sigmoid Funktion von der Aktivierungsstärke Zv

theta1c = Vb' * Zvd; % Änderunswunsch für Theta1

% weitere Ableitungen sind nicht notwendig

theta1 -= theta1c .* alpha; % ändere die Gewichte von Theta1 und Theta2

theta2 -= theta2c .* alpha; % der Änderungswunsch wird von der Lernrate abgeschwächt

endfor

% Ausgabe von der letzten Vorhersage und den Gewichten

printf("abgeschlossen. \n")

printf("Letzte Vorhersage und trainierte Gewichte\n")

theta1

theta2

Zu jeder Zeile bzw. Funktion die wir im Vorwärtspass geschrieben haben, gibt es im Rückwärtspass eine abgeleitete Variante. Dank den Ableitungen können wir die Änderungswünsche der Gewichte in jeder Schicht ausrechnen und am Ende einer Trainingsiteration anwenden. Wir trainieren 10.000 Iterationen lang und verwenden eine Lernrate von 0,8. In komplexeren Fragestellungen, mit mehr Daten, würden diese Werte niedriger ausfallen.

Es ist außerdem möglich den ganzen Programmcode viel modularer aufzubauen. Dazu werde ich im nächsten Artikel auf eine mehr objekt-orientiertere Sprache wechseln. Nichts desto trotz liefert der obige Algorithmus gute Ergebnisse. Hier ist mal ein Ausgabebeispiel:

X =

0 0

0 1

1 0

1 1

Y =

theta1 =

0.114950 0.046125

0.064683 0.139159

-0.164288 -0.094688

theta2 =

0.33607

-0.31128

0.13993

m = 4

iteration = 10000

alpha = 0.80000

Starte Training ... abgeschlossen.

Letzte Vorhersage und trainierte Gewichte

O =

0.014644

0.983308

0.986137

0.013060

theta1 =

3.2162 -3.0431

6.4365 5.6498

-6.3383 -5.8602

theta2 =

4.4759

-9.5057

9.9795

Data Science mit Neo4j und R

February 18, 2016/in Data Mining, Data Science, Data Science Hack, Database, Graph Database, Neo4J, NoSQL, R Statistics, Tools, Tutorial/by Miachael Hunger

Traurig, aber wahr: Data Scientists verbringen 50-80% ihrer Zeit damit, Daten zu bereinigen, zu ordnen und zu bearbeiten. So bleibt nur noch wenig Zeit, um tatsächlich vorausschauende Vorhersagemodelle zu entwickeln. Vor allem bei klassischen Stacks, besteht die Datenanalyse zum Großteil darin, Zeile für Zeile in SQL zu überführen. Zeit zum Schreiben von Modell-Codes in einer statistischen Sprache wie R bleibt da kaum noch. Die langen, kryptischen SQL-Abfragen verlangsamen aber nicht nur die Entwicklungszeit. Sie stehen auch einer sinnvollen Zusammenarbeit bei Analyse-Projekten im Weg, da alle Beteiligten zunächst damit beschäftigt sind, die SQL-Abfragen der jeweils anderen zu verstehen.

Komplexität der Daten steigt

Der Grund für diese Schwierigkeiten: Die Datenstrukturen werden immer komplexer, die Vernetzung der Daten untereinander nimmt immer stärker zu. Zwängt man diese hochgradig verbundenen Datensätze in eine SQL-Datenbank, in der Beziehungen naturgemäß abstrakt über Fremdschlüssel dargestellt werden, erhält man als Ergebnis übermäßig komplizierte Schematas und Abfragen. Als Alternative gibt es jedoch einige NoSQL-Lösungen – allen voran Graphdatenbanken – die solche hochkomplexen und heterogenen Daten ohne Informationsverlust speichern können – und zwar nicht nur die Entitäten an sich, sondern auch besonders die Beziehungen der Daten untereinander.

Datenanalysen zielen immer stärker darauf ab, das Verhalten und die Wünsche von Kunden besser verstehen zu können. Die Fragen lauten z. B.:

Wie hoch ist die Wahrscheinlichkeit, dass ein Besucher auf eine bestimmte Anzeige klickt?
Welcher Kunde sollte in welchem Kontext welche Produktempfehlungen erhalten?
Wie kann man aus der bisherigen Interaktionshistorie des Kunden sein Ziel vorhersagen, bevor er selbst dort ankommt?
In welchen Beziehungen steht Nutzer A zu Nutzer B?

Menschen sind bekanntermaßen von Natur aus sozial. Einige dieser Fragen lassen sich daher beantworten, wenn man weiß, wie Personen miteinander in Verbindung stehen: Unsere Zielperson, Nutzer A ähnelt in seinem Kontext und Verhalten Benutzer B. Und da Benutzer B ein bestimmtes Produkt (z. B. ein Spielfilm) gefällt, empfehlen wir diesen Film auch Nutzer A. In diese Auswertung fließen natürlich auch noch weitere Faktoren mit ein, z. B. die Demographie und der soziale Status des Nutzers, seine Zuordnung zu Peer Groups, vorher gesehene Promotions oder seine bisherigen Interaktionen.

Visualisierung eines Graphen mit RNeo4j

Mit R und Neo4j lassen sich Graphen und Teilgraphen ganz einfach mit RNeo4j, igraph und visNetwork libraries visualisieren.

library(igraph)

library(visNetwork)

library(RNeo4j)

Das folgende Beispiel zeigt wie in einem Graphen Schauspieler und Filme sowie ihre Beziehungen zueinander anschaulich dargestellt werden können, z. B. um Empfehlungen innerhalb eines Filmportals zu generieren. Dabei sind zwei Schauspieler über eine Kante miteinander verbunden, wenn sie beide im gleichen Film mitspielen.

Im ersten Schritt werden dazu in Neo4j die Film-Datensätze importiert (Achtung: Dieser Vorgang löscht die aktuelle Datenbank).

graph = startGraph("http://localhost:7474/db/data/")

importSample(graph, "movies", input=F)

Als nächstes wird mit Cypher eine entsprechende Liste von Beziehungen aus Neo4j gezogen. Wie man sehen kann, ist die Darstellung des gewünschten Graph-Musters innerhalb der Abfrage sehr anschaulich.

query = "

MATCH (p1:Person)-[:ACTED_IN]->(:Movie)<-[:ACTED_IN]-(p2:Person)

WHERE p1.name < p2.name

RETURN p1.name AS from, p2.name AS to, COUNT(*) AS weight

edges = cypher(graph, query)

head(edges)

## from to weight

## 1 Brooke Langton Keanu Reeves 1

## 2 Jack Nicholson Kevin Bacon 1

## 3 Jerry O'Connell Kiefer Sutherland 1

## 4 Oliver Platt Sam Rockwell 1

## 5 John Goodman Susan Sarandon 1

## 6 Gary Sinise Kevin Bacon 1

Die visNetwork Funktion erwartet sowohl Kanten-Dataframes als auch Knoten-Dataframes. Ein Knoten-Dataframe lässt sich daher über die eindeutigen Werte des Kanten-Dataframes generieren.

nodes = data.frame(id=unique(c(edges$from, edges$to)))

nodes$label = nodes$id

head(nodes)

## id label

## 1 Brooke Langton Brooke Langton

## 2 Jack Nicholson Jack Nicholson

## 3 Jerry O'Connell Jerry O'Connell

## 4 Oliver Platt Oliver Platt

## 5 John Goodman John Goodman

## 6 Gary Sinise Gary Sinise

Im Anschluss können die Knoten- und Kanten-Dataframes in das visNetwork übertragen werden.
visNetwork(nodes, edges)

Nun kommt igraph mit ins Spiel, eine Bibliothek von Graph-Algorithmen. Durch Einbindung der Kantenliste lässt sich einfach ein igraph Graph-Objekt erstellen, das den Teilgraphen miteinschließt.

ig = graph_from_data_frame(edges, directed=F)

## IGRAPH UNW- 102 362 --

## + attr: name (v/c), weight (e/n)

## + edges (vertex names):

## [1] Brooke Langton --Keanu Reeves

## [2] Jack Nicholson --Kevin Bacon

## [3] Jerry O'Connell --Kiefer Sutherland

## [4] Oliver Platt --Sam Rockwell

## [5] John Goodman --Susan Sarandon

## [6] Gary Sinise --Kevin Bacon

## [7] J.T. Walsh --Noah Wyle

## [8] Jim Broadbent --Tom Hanks

## + ... omitted several edges

Die Größe der Knoten kann als Funktion der Edge-Betweeness-Centrality definiert werden. In visNetwork entspricht dabei jede “value”-Spalte im Knoten-Dataframe der Größe des Knoten.
nodes$value = betweenness(ig)

head(nodes)

## id label value

## 1 Brooke Langton Brooke Langton 0.000000

## 2 Jack Nicholson Jack Nicholson 511.443714

## 3 Jerry O'Connell Jerry O'Connell 154.815234

## 4 Oliver Platt Oliver Platt 20.643840

## 5 John Goodman John Goodman 1.659259

## 6 Gary Sinise Gary Sinise 33.723499

Mit Einführung der “Value”-Spalte werden die Knoten nun alle unterschiedlich groß dargestellt.
visNetwork(nodes, edges)

Mit Hilfe eines Community-Detection-Algorithmus lassen sich im Graphen nun Cluster finden. In diesem Beispiel wird der „Girvan-Newman”-Algorithmus verwendet, der in igraph als cluster_edge_betweenness bezeichnet wird.

clusters = cluster_edge_betweenness(ig)

clusters[1:2]

## $`1`

## [1] "Brooke Langton" "Liv Tyler" "Charlize Theron"

## [4] "Emil Eifrem" "Dina Meyer" "Diane Keaton"

## [7] "Keanu Reeves" "Gene Hackman" "Ice-T"

## [10] "Al Pacino" "Carrie-Anne Moss" "Clint Eastwood"

## [13] "Orlando Jones" "Takeshi Kitano" "Laurence Fishburne"

## [16] "Richard Harris"

## $`2`

## [1] "Jack Nicholson" "Jerry O'Connell" "J.T. Walsh"

## [4] "Renee Zellweger" "Kiefer Sutherland" "Cuba Gooding Jr."

## [7] "Marshall Bell" "Aaron Sorkin" "Kevin Bacon"

## [10] "Kevin Pollak" "Christopher Guest" "Demi Moore"

## [13] "Regina King" "Kelly Preston" "John Cusack"

## [16] "Danny DeVito" "Bonnie Hunt" "Corey Feldman"

## [19] "Jay Mohr" "James Marshall" "Jonathan Lipnicki"

## [22] "River Phoenix" "Tom Cruise" "Noah Wyle"

## [25] "Wil Wheaton" "John C. Reilly"

In der Liste oben sind alle Schauspieler der ersten zwei Cluster zu sehen. Insgesamt konnten sechs Cluster identifiziert werden.

1 2	length(clusters) ## [1] 6

Durch Hinzufügen einer “Group”-Spalte im Knoten-Dataframe, werden alle Knoten in visNetwork entsprechend ihrer Gruppenzugehörigkeit farblich markiert. Diese Cluster-Zuordnung erfolgt über clusters$membership. Durch Entfernen der “Value”-Spalte lassen sich die Knoten wieder auf eine einheitliche Größe bringen.

nodes$group = clusters$membership

nodes$value = NULL

head(nodes)

## id label group

## 1 Brooke Langton Brooke Langton 1

## 2 Jack Nicholson Jack Nicholson 2

## 3 Jerry O'Connell Jerry O'Connell 2

## 4 Oliver Platt Oliver Platt 3

## 5 John Goodman John Goodman 4

## 6 Gary Sinise Gary Sinise 3

Werden die Knoten- und Kanten-Datenframes erneut in visNetwork übertragen, sind nun alle Knoten eines Clusters in derselben Farbe dargestellt.
visNetwork(nodes, edges)

Mit diesem Workflow lassen sich Teilgraphen in Neo4j einfach abfragen und Cluster-Algorithmen einfach darstellen.

Generell eignen sich Graphdatenbanken wie Neo4j besonders gut, um stark vernetzte und beliebig strukturierte Informationen zu handhaben – egal ob es sich um Schauspieler, Filme, Kunden, Produkte, Kreditkarten oder Bankkonten handelt. Zudem können sowohl den Knoten als auch den Kanten beliebige qualitative und quantitative Eigenschaften zugeordnet werden. Beziehungen zwischen Daten sind also nicht mehr bloße Strukturinformationen, sondern stehen vielmehr im Zentrum des Modells.

Cypher: intuitiv nutzbare Programmiersprache

Die Zeiten, in denen Data Science zum Großteil aus Datenbereinigung und -mapping besteht, sind damit vorbei. Mit dem entsprechenden Ansatz laufen Entwicklungsprozesse deutlich schneller und einfacher ab. Data Scientists kommen mit weniger Code schneller ans Ziel und können mehr Zeit in das tatsächliche Entwickeln von relevanten Modellen investieren. Dabei nutzen sie die Flexibilität einer quelloffenen NoSQL-Graphdatenbank wie Neo4j kombiniert mit der Reife und weiten Verbreitung der Statistiksprache R für statistisches Rechnen und Visualisierung. Programmierer müssen nicht mehr stundenlang komplexe SQL-Anweisungen schreiben oder den ganzen Tag damit verbringen, eine Baumstruktur in SQL zu überführen. Sie benutzen einfach Cypher, eine musterbasierte, für Datenbeziehungen und Lesbarkeit optimierte Abfragesprache und legen los.

Die ersten Schritte

Neue Variablen

Aggregate

Daten in Reih’ und Glied

Fazit und Ausblick

Die Vorarbeit

ABC – Analyse:

Klassifikation:

XYZ – Analyse

Klassifikation:

ABC-XYZ-Analyse

Benford-Algorithmus in Python mit NumPy und Pandas

Analyse der Verteilung der ersten Ziffer in Zahlungsdaten

Abgleich der Ziffernhäufigkeit mit der Erwartung

Gegenüberstellung: Computer-generierte Zufallszahlen

Anwendung in der Praxis

Die Rastrigin-Funktion

Die Rastrigin-Funktion im Standard-Python umsetzen und visualisieren

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

Die ideale Lösung

Schritt Eins – Integration von Hadoop mit Active Directory

Schritt Zwei – Synchronisation der Rechte & Rollen von Active Directory

Schritt Drei – Anlegen und Verwaltung von technischen Usern

High Level Architektur

Architektur in einer Mehrclusterumgebung

Zusammenfassung

Der aktuelle Stand der Technologie

Alle einzelne Stücke des Puzzles kommen zusammen

Klassifikation via K-Nearest Neighbour Algorithmus

Vorab: Verwendete Bibliotheken

Schritt 1 – Daten betrachten und Merkmale erkennen

Schritt 2 – Merkmale im Verhältnis zueinander perspektivisch betrachten

Schritt 3 – Kalkulation der Distanzen zwischen den einzelnen Punkten

Schritt 4 & 5 – Klassifikation durch Eingrenzung auf k-nächste Nachbarn

Fehlerraten herabsenken

Das Problem mit den Dimensionen

Weiterführende Literatur

Training

visuelle Darstellung aller Gewichtskombinationen und deren Vorhersagefehler

Anwendungsbeispiel und Programmcode

Komplexität der Daten steigt

Visualisierung eines Graphen mit RNeo4j

Cypher: intuitiv nutzbare Programmiersprache

Interesting links

Pages

Categories

Archive