Daniela Keller

Warenkorbanalyse in R

/in , , , , , , /by

Was ist die Warenkorbanalyse?

Die Warenkorbanalyse ist eine Sammlung von Methoden, die die beim Einkauf gemeinsam gekauften Produkte oder Produktkategorien aus einem Handelssortiment untersucht. Ziel der explorativen Warenkorbanalyse ist es, Strukturen in den Daten zu finden, so genannte Regeln, die beschreiben, welche Produkte oder Produktkategorien gemeinsam oder eben nicht gemeinsam gekauft werden.

Beispiel: Wenn ein Kunde Windeln und Bier kauft, kauft er auch Chips.

Werden solche Regeln gefunden, kann das Ergebnis beispielsweise für Verbundplatzierungen im Verkaufsraum oder in der Werbung verwendet werden.

Datenaufbau

Die Daten, die für diese Analyse untersucht werden, sind Transaktionsdaten des Einzelhandels. Meist sind diese sehr umfangreich und formal folgendermaßen aufgebaut:

data-bsp

Ausschnitt eines Beispieldatensatzes: Jede Transaktion (= Warenkorb = Einkauf) hat mehrere Zeilen, die mit der selben Transaktionsnummer (Spalte Transaction) gekennzeichnet sind. In den einzelnen Zeilen der Transaktion stehen dann alle Produkte, die sich in dem Warenkorb befanden. In dem Beispiel sind zudem noch zwei Ebenen von Produktkategorien als zusätzliche Informationen enthalten.

Es gibt mindestens 2 Spalten: Spalte 1 enthält die Transaktionsnummer (oder die Nummer des Kassenbons, im Beispielbild Spalte Transaction), Spalte 2 enthält den Produktnamen. Zusätzlich kann es weitere Spalten mit Infos wie Produktkategorie, eventuell in verschiedenen Ebenen, Preis usw. geben. Sind Kundeninformationen vorhanden, z.B. über Kundenkarten, so können auch diese Informationen enthalten sein und mit ausgewertet werden.

Beschreibende Datenanalyse

Die Daten werden zunächst deskriptiv, also beschreibend, analysiert. Dazu werden z.B. die Anzahl der Transaktionen und die Anzahl der Produkte im Datensatz berechnet. Zudem wird die Länge der Transaktionen, also die Anzahl der Produkte in den einzelnen Transaktionen untersucht. Dies wird mit deskriptiven Maßzahlen wie Minimum, Maximum, Median und Mittelwert in Zahlen berichtet sowie als Histogramm grafisch dargestellt, siehe folgende Abbildung.

hist-sizes
Histogramm der Längenverteilung der Transaktionen.

Die häufigsten Produkte werden ermittelt und können gesondert betrachtet werden. Als Visualisierung kann hier ein Balkendiagramm mit den relativen Häufigkeiten der häufigsten Produkte verwendet werden, wie im folgenden Beispiel.

relfreq-items
Relative Häufigkeiten der häufigsten Produkte, hier nach relativer Häufigkeit größer 0,1 gefiltert.

Ähnliche Analysen können bei Bedarf auch auf Kategorien-Ebene oder nach weiteren erhobenen Merkmalen selektiert durchgeführt werden, je nachdem, welche Informationen in den Daten stecken und welche Fragestellungen für den Anwender interessant sind.

Verbundanalyse

Im nächsten Schritt wird mit statistischen Methoden nach Strukturen in den Daten gesucht, auch Verbundanalyse genannt. Als Grundlage werden Ähnlichkeitsmatrizen erstellt, die für jedes Produktpaar die Häufigkeit des gemeinsamen Vorkommens in Transaktionen bestimmen. Solch eine Ähnlichkeitsmatrix ist zum Beispiel eine Kreuztabelle in der es für jedes Produkt eine Spalte und eine Zeile gobt. In den Zellen in der Tabelle steht jeweils die Häufigkeit, wie oft dieses Produktpaar gemeinsam in Transaktionen in den Daten vorkommt, siehe auch folgendes Beispiel.

screenshot-crosstable-ausschnitt

Ähnlichkeitsmatrix oder Kreuztabelle der Produkte: Frankfurter und Zitrusfrüchte werden in 64 Transaktionen zusammen gekauft, Frankfurter und Berries in 22 usw.

Auf Basis solch einer Ähnlichkeitsmatrix wird dann z.B. mit Mehrdimensionaler Skalierung oder hierarchischen Clusteranalysen nach Strukturen in den Daten gesucht und Gemeinsamkeiten und Gruppierungen gefunden. Die hierarchische Clusteranalyse liefert dann ein Dendrogram, siehe folgende Abbildung, in der ähnliche Produkte miteinander gruppiert werden.

dendrogram

Dendrogram als Visualisierung des Ergebnisses der hierarchischen Clusterananlyse. Ähnliche Produkte (also Produkte, die zusammen gekauft werden) werden zusammen in Gruppen geclustert. Je länger die vertikale Verbindungslinie ist, die zwei Gruppen oder Produkte zusammen fasst, um so unterschiedlicher sind diese Produkte bzw. Gruppen.

Assoziationsregeln

Schließlich sollen neben den Verbundanalysen am Ende in den Daten Assoziationsregeln gefunden werden. Es werden also Regeln gesucht und an den Daten geprüft, die das Kaufverhalten der Kunden beschreiben. Solch eine Regel ist zum Beispiel „Wenn ein Kunde Windeln und Bier kauft, kauft er auch Chips.“ Formal: {Windeln, Bier} → {Chips}

Für diese Regeln lassen sich statistische Maßzahlen berechnen, die die Güte und Bedeutung der Regeln beschreiben. Die wichtigsten Maßzahlen sind Support, Confidence und Lift:

Support ist das Signifikanzmaß der Regel. Es gibt an, wie oft die gefundene Regel in den Daten anzuwenden ist. Wie oft also die in der Regel enthaltenen Produkte gemeinsam in einer Transaktion vorkommen. In dem Beispiel oben: Wie oft kommen Windeln, Bier und Chips in einer Transaktion gemeinsam vor?

Confidence ist das Qualitätsmaß der Regel. Es beschreibt, wie oft die Regel richtig ist. In dem oben genanten Beispiel: Wie oft ist in einer Transaktion Chips enthalten, wenn auch Windeln und Bier enthalten sind?

Lift ist das Maß der Bedeutung der Regel. Es sagt aus wie oft die Confidence den Erwartungswert übersteigt. Wie ist die Häufigkeit des gemeinsamen Vorkommens von Windeln, Bier und Chips im Verhältlnis zur erwarteten Häufigkeit des Vorkommens, wenn die Ereignisse stochastisch unabhängig sind?

Algorithmen

In den Daten werden zunächst alle möglichen Regeln gesammelt, die einen Mindestwert an Support und Confidence haben. Die Mindestwerte werden dabei vom Nutzer vorgegeben. Da es sich bei Transaktionsdaten um große Datenmengen handelt und häufig große Anzahlen von Produkten enthalten sind, wird die Suche nach Regeln zu einem komplexen Problem. Es wurden verschiedene effiziente Algorithmen als Suchstrategien entwickelt, z.B. der APRIORI-Algorithmus von Agrawal und Srikant (1994), der auch im weiter unten vorgestellten Paket arules von R verwendet wird.

Sind die Assoziationsregeln gefunden, können Sie vom Nutzer genauer untersucht werden und z.B. nach den oben genannten Kennzahlen sortiert betrachtet werden, oder es werden die Regeln für spezielle Warenkategorien genauer betrachtet, siehe folgendes Beispiel.

screenshot-rules

Beispielausgabe von Regeln, hier die drei Regeln mit dem besten Lift. In der ersten Regel sieht man: Wenn Bier und Wein gekauft wird, wird auch Likör gekauft. Diese Regel hat einen Support von 0,002. Diese drei Produkte kommen also in 0,2 % der Transaktionen vor. Die Confidence von 0,396 zeigt, dass in 39,6 % der Transaktionen auch Likör gekauft wird, wenn Bier und Wein gekauft wird.

Umsetzung mit R

Die hier vorgestellten Methoden zur Warenkorbanalyse lassen sich mit dem Paket arules der Software R gut umsetzen. Im Folgenden gebe ich eine Liste von nützlichen Befehlen für diese Analysen mit dieser Software. Dabei wird mit data hier durchgehend der Datensatz der Transaktionsdaten bezeichnet.

summary(data)

Zusammenfassung des Datensatzes:

  • Anzahl der Transaktionen und Anzahl der Warengruppen
  • die häufigsten Produkte werden genannt mit Angabe der Häufigkeiten
  • Längenverteilung der Transaktionen (Anzahl der Produkte pro Transaktion): Häufigkeiten, deskriptive Maße wie Quartile
  • Beispiel für die Datenstruktur (Levels)
size(data)

Längen der Transaktionen (Anzahl der Produkte pro Transaktion)

hist(size(data))

Histogramm als grafische Darstellung der Transaktionslängen

itemFrequencyPlot(data, support=0.1)

rel. Häufigkeiten der einzelnen Produkte, hier nur die mit mindestens 10 % Vorkommen

crossTable(data)

Äquivalenzmatrix: Häufigkeiten der gemeinsamen Käufe für Produktpaare

dissJacc <- dissimilarity(data[, itemFrequency(data) > 0.05], method = "Jaccard", which = "items")

Unähnlichkeitsmatrix für die hierarchische Clusteranalyse

hcWard <- hclust(dissJacc, method = "ward.D")

Hierarchische Clusteranalyse

plot(hcWard)

Dendrogram der hierarchischen Clusteranalyse

rules <- apriori(data, parameter = list(support = 0.001, confidence = 0.2), control = list(verbose = FALSE))

Assoziationsregeln finden mit APRIORI-Algorithmus, hier Regeln mit mindestens 1% Support und 20 % Confidence

summary(rules)

Zusammenfassung der oben gefundenen Regeln (Anzahl, Eigenschaften Support, Confidence, Lift)

inspect(SORT(rules,by=“lift“)[1:5])

Einzelne Regeln betrachten, hier die laut Lift besten 5 Regeln

Referenzen:

  • Michael Hahsler, Kurt Hornik, Thomas Reutterer: Warenkorbanalyse mit Hilfe der Statistik-Software R, Innovationen in Marketing, S.144-163, 2006.
  • Michael Hahsler, Bettina Grün, Kurt Hornik, Christian Buchta, Introduciton to arules – A computational environment for mining association rules and frequent item sets. (Link zum PDF)
  • Rakesh Agrawal, Ramakrishnan Srikant, Fast algorithms for mining association rules, Proceedings of the 20th VLDB Conference Santiago, Chile, 1994
  • Software R:  R Core Team (2016). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Link: R-Project.org.
  • Paket: arules: Mining Association Rules using R.

Beispieldatensatz: Groceries aus dem Paket arules

Benjamin Aunkofer

Benford-Analyse

/in , , , , /by

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

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

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

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

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

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

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

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

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

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

import math

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

Benford-Algorithmus in Python mit NumPy und Pandas

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

import numpy as np

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

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

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

import pandas as pd

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

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

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

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

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

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

Analyse der Verteilung der ersten Ziffer in Zahlungsdaten

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

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

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

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

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

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

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

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

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

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

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

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

Abgleich der Ziffernhäufigkeit mit der Erwartung

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

import matplotlib.pyplot as pyplot 

fig = pyplot.figure()

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

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

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

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

pyplot.show()

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

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

Gegenüberstellung: Computer-generierte Zufallszahlen

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

from random import randint

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

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

benford-analyse-python-numpy-pseudo-zufallszahlen

Anwendung in der Praxis

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

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

Benjamin Aunkofer

Die Rastrigin-Funktion

/in , , , , /by

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

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

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

Die Rastrigin-Funktion

Mathematisch beschrieben wird die Rastrigin-Funktion wie folgt:

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

-5.12 leq x_i leq 5.12

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

Die Rastrigin-Funktion im Standard-Python umsetzen und visualisieren

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

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

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

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

rastriginValues = []
i = 0


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

rastrigin-line-chart

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

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

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

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

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

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

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

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

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

rastrigin-funktion-python-3d-plot

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

 

 

Nico Hezel

KNN: Rückwärtspass

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

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("nnStarte 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 Gewichten")
 O
 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 =                                                                                                                                                                                                                                                                                                                                                                                                                   
   0                                                                                                                                                                                                              
   1                                                                                                                                                                                                              
   1                                                                                                                                                                                                              
   0                                                                                                                                                                                                              
                                                                                                                                                                                                                  
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

 

Dr. Stefan Kühn

Wie lernen Maschinen?

/in , , , , , , /by

Im dritten Teil meiner Reihe Wie lernen Maschinen? wollen wir die bisher kennengelernten Methoden anhand eines der bekanntesten Verfahren des Maschinellen Lernens – der Linearen Regression – einmal gegenüberstellen. Die Lineare Regression dient uns hier als Prototyp eines Verfahrens aus dem Gebiet der Regression, in weiteren Artikeln werden die Logistische Regression als Prototyp eines Verfahrens aus dem Gebiet der Klassifikation und eine Collaborative-Filtering- bzw. Matrix-Faktorisierungs-Methode als Prototyp eines Recommender-Systems behandelt.

Read more

Benjamin Aunkofer

Wahrscheinlichkeitesrechnung – Grundstein für Predictive Analytics

/in , , , /by

Die Wahrscheinlichkeitsrechnung behandelt die Gesetzmäßigkeiten  des (von außen betrachtet) zufälligen Vorkommens bestimmter Ereignisse aus einer vorgegebenen Ereignismenge. Die mathematische Statistik fasst diese Wahrscheinlichkeitsrechnung zur Stochastik zusammen, der Mathematik des Zufalls. 

Mit diesem Artikel – zu der ich eine Serie plane – möchte ich den Einstieg in Predictive Analytics wagen, zugegebenermaßen ein Themengebiet, in dem man sich sehr schnell verlieren und den Wald vor lauter Bäumen nicht mehr findet. Also belassen wir es erstmal bei einem sanften Einstieg…

Klassische Definition der Wahrscheinlichkeit

Das klassische Verständnis der Wahrscheinlichkeit geht von endlich vielen Ausgängen (Ereignisse) aus, bei denen alle Ausgänge gleich wahrscheinlich sind. Die dafür erdachten Zufallsexperimente wurden von dem französischen Mathematiker Pierre Simon Lapplace (1749 – 1827) zum ersten Mal nachvollziehbar beschrieben. Diese Zufallsexperimente werden daher auch Laplace-Experimente genannt.

Bei einem Laplace Experiment gilt:

Ereignismenge Omega = {omega_1,omega_2,omega_3,…omega_s}
Wahrscheinlichkeit p(w_j)=frac{1}{s}=frac{1}{|Omega|}
(j=1,2,3,…s)

Die Ergebnismenge, das ist die Menge aller möglichen Ereignisse, wird in der Regel mit einem Omega (Omega) gekennzeichnet, ein beliebiges Einzelereignis hingegen als omega (kleines Omega).

Eine typische Laplace-Wahrscheinlichkeitsfrage ist ein bevorstehender Würfelwurf. Wie groß ist die Wahrscheinlichkeit, mit einem echten (unverfälschten) Würfel eine gerade Zahl zu würfeln?

Mit Omega={1,2,3,4,5,6} und A={2,4,6} folgt P(A)=frac{|A|}{|Omega|}=frac{3}{6}=0,5.

Axiomatische Definition der Wahrscheinlichkeit

Jeder Wahrscheinlichkeitsbegriff muss auf denselben äußeren Bedingungen beruhenden Zufallsexperimenten beliebig oft wiederholbar sein. Die axiomatische Definition der Wahrscheinlichkeit P(A) eines Ereignisses A berücksichtigt Axiome. Axiome sind nicht beweisbare Grundpostulate, darunter fallen Gegebenheiten, die gewissermaßen unverstanden sind und deren Vorkommen und Bedeutung in der Regel empirisch belegt werden müssen.
Die Definition der axiomatischen Wahrscheinlichkeit stammt vom russischen Mathematiker Andrej Nikollajewitsch Kolmogorov (1903 – 1987).

In der Realität gibt es keine perfekte Zufälligkeit, denn jedes Ergebnis ist von ganz bestimmten Faktoren abhängig. Auf den Würfelwurf bezogen, hängt das gewürfelte Ergebnis von unüberschaubar vielen Faktoren ab. Wären diese alle bekannt, könnte das Ergebnis exakt berechnet und somit mit einer Sicherheit vorhergesagt werden. Da dafür jedoch in der Praxis unbestimmbar viele Faktoren eine Rolle spielen (beispielsweise die genaue Beschaffenheit des Würfels in Form, Gewicht, Materialwiderstand, der genaue Winkel, die Fallgeschwindigkeit, die Ausgangsposition der Hand und des Würfels) können wir das Ergebnis nur schätzen, indem die Beschreibung des Vorgangs vereinfacht wird. Nur diese Vereinfachung macht es uns möglich, Vorhersagen zu treffen, die dann jedoch nur eine Wahrscheinlichkeit darstellen und somit mit einer Unsicherheit verbunden sind.

In der abstrakten Welt des perfekten Zufalls gäbe es die gleiche Chance, eine “4” zu würfeln, wie jeweils alle anderen Ziffern.

Mit Omega={1,2,3,4,5,6} und A={4} folgt P(A)=frac{|A|}{|Omega|}=frac{1}{6}=0,167.

Das Ergebnis eines Wurfes des Würfels ist in der Realität auch von der Beschaffenheit des Würfels abhängig. Angenommen, der Würfel hat auf Seite der Ziffer “4” bei allen vier Kanten eine Abrundung, die ein Umkippen auf eine andere Seite begünstigen, so bedeutet dies:

  • Die Ziffer “4” hat vier abgerundete Kanten, die Wahrscheinlichkeit eine “4” zu würfeln sinkt stark
  • Die Ziffern “1”, “3”, “5”, “6” haben jeweils eine abgerundete Kante (Berühungskante zur “4”) sinkt
  • Die Ziffer “2” liegt der “4” gegenüber, hat somit keine Berührungskante und keine Abrundung, so steigt ihre Chance gewürfelt zu werden

Nun könnte sich nach einer empirischen Untersuchung mit einer ausreichenden Stichprobe folgende Wahrscheinlichkeit ergeben:

  • p(4) = 0,1
  • p(1) = p(3) = p(5) = p(6) = 0,15
  • p(2) = 0,3
  • P(Omega) = 1,0

Durch die Analyse der bisherigen Wurf-Historie und der Betrachtung der Beschaffenheit der Kanten des Würfels können wir uns somit weit realistischere Wahrscheinlichkeiten über die Wurfergebnisse ermitteln. Wie hoch wäre nun die Wahrscheinlichkeit, nach einem Wurf eine gerade Zahl zu würfeln?

Mit Omega={1,2,3,4,5,6} und A={2,4,6} folgt P(A)=p(2)+p(4)+p(6)=0,55.