Die Familie der apply-Funktionen in R Teil 3: Weitere mit apply() verwandte Funktionen

Einordnung des Artikels

• Einführung in die Informatik
  • Einführung in die Programmiersprache R
    • Funktionen in R
      • Die Funktion apply() in R: Iteration über die Zeilen oder Spalten einer Matrix
      • Die Familie der apply-Funktionen in R Teil 1: Verarbeitung von Listen mit lapply(), sapply(), vapply() und rapply()
      • Die Familie der apply-Funktionen in R Teil 2: Die Verarbeitung mehrerer Listen mit mapply(), Map() und outer()
      • Die Familie der apply-Funktionen in R Teil 3: Weitere mit apply() verwandte Funktionen

Die in diesem Artikel vorgestellten Funktionen setzen Kenntnisse über Dataframes und Faktoren voraus, die hier nicht erklärt werden; sie finden sich in:

1. Dataframes in R: der Datentyp data frame
2. Dataframes in R: Anwendungen
3. Faktoren in R: der Datentyp factor
4. Faktoren in R: Anwendungen. Hier wird insbesondere tapply() ausführlich erklärt.

Einführung

Die folgende Tabelle zeigt, welche Funktionen in diesem Kapitel besprochen werden, die Eingabewerte der Funktionen sind nicht angegeben:

Die Tabelle zeigt 2 Besonderheiten:

1. Die Funktionen zerfallen in 2 Gruppen:
   • Funktionen, die Objekte oder Anweisungen wiederholen.
   • Funktionen, die zuerst Daten gruppieren und darauf eine Funktion anwenden
2. Alle Funktionen der zweiten Gruppe scheinen identische Aufgaben zu erfüllen. Im Folgenden werden ihre Unterschiede besprochen, die meist den Rückgabewert betreffen.

Für alle Funktionen werden ihre Ein- und Rückgabewerte sowie die Unterschiede in der Verarbeitung der eingegebenen Daten besprochen. Dazu werden typische Anwendungen vorgestellt.

Die Funktion rep()

Die Funktion rep()

rep(x, times = 1, length.out = NA, each = 1)

wurde bei Vektoren beschrieben; am häufigsten wird sie eingesetzt, um Kopien eines Objektes zu erzeugen, wie in rep(x = 17, times = 5) .

Man kann rep() aber auch im folgenden Sinn zur Familie der apply()-Funktionen zählen: Derartige Funktionen sollen Schleifen ersetzen beziehungsweise deren Konfiguration erleichtern. Und oft werden Schleifen eingesetzt, um bei jedem Durchlauf eine Kopie eines Objektes zu erzeugen; oder es wird abhängig vom Zählindex eine Auswahl aus mehreren Objekten getroffen. Alle in den den Schleifendurchläufen erzeugten Objekte werden zu einem Objekt (meist einer Liste) zusammengefasst.

Um nicht nur triviale Wiederholungen zu erzeugen, wie durch rep(x = 15, times = 5) , sondern komplexe Strukturen aufzubauen, kann man

• die Argumente times und each gleichzeitig verwenden,
• für das zu wiederholende Objekt x nicht nur elementare Datentypen sondern Listen einsetzen.

Die folgenden Beispiele sollen einen Eindruck vermitteln, welche Einsatzmöglichkeiten rep() bietet.

Ausgangspunkt ist eine Liste lst mit drei Komponenten, wobei jede Komponente ein Vektor der Länge drei ist:

lst <- vector(mode = "list", length = 3)
lst[[1]] <- (1:3)
lst[[2]] <- (4:6)
lst[[3]] <- (7:9)
str(lst)
# List of 3
# $ : int [1:3] 1 2 3
# $ : int [1:3] 4 5 6
# $ : int [1:3] 7 8 9

Wird nur das Argument times gesetzt, werden die Wiederholungen in eine neue Liste verpackt:

lst.2 <- rep(x = lst, times = 2)
str(lst.2)
# List of 6
# $ : int [1:3] 1 2 3
# $ : int [1:3] 4 5 6
# $ : int [1:3] 7 8 9
# $ : int [1:3] 1 2 3
# $ : int [1:3] 4 5 6
# $ : int [1:3] 7 8 9

Das Argument times gibt somit an, wie oft die Liste als Ganzes wiederholt wird. Das Argument each gibt an, wie oft die einzelnen Komponenten der Liste wiederholt werden. Im folgenden Skript wird zuerst nur das Argument each gesetzt, anschließend, times und each:

lst.3 <- rep(x = lst, each = 2)
str(lst.3)
# List of 6
# $ : int [1:3] 1 2 3
# $ : int [1:3] 1 2 3
# $ : int [1:3] 4 5 6
# $ : int [1:3] 4 5 6
# $ : int [1:3] 7 8 9
# $ : int [1:3] 7 8 9

lst.2.2 <- rep(x = lst, times = 2, each = 2)
str(lst.2.2)
# List of 12
# $ : int [1:3] 1 2 3
# $ : int [1:3] 1 2 3
# $ : int [1:3] 4 5 6
# $ : int [1:3] 4 5 6
# $ : int [1:3] 7 8 9
# $ : int [1:3] 7 8 9
# $ : int [1:3] 1 2 3
# $ : int [1:3] 1 2 3
# $ : int [1:3] 4 5 6
# $ : int [1:3] 4 5 6
# $ : int [1:3] 7 8 9
# $ : int [1:3] 7 8 9

Die Beispiele zeigen noch nicht die volle Flexibilität, die rep() ermöglicht: Da die Liste lst drei Komponenten hat, kann man dem Argument times drei Werte vorgeben, die bestimmen, wie oft jede Komponente wiederholt werden soll. Das folgende Skript setzt im Argument times einen Vektor mit den drei Komponenten (2, 3, 4), wodurch

• die erste Komponente von lst zweimal wiederholt wird,
• die zweite Komponente von lst dreimal wiederholt wird und
• die dritte Komponente von lst viermal wiederholt wird.

lst.4 <- rep(x = lst, times = c(2, 3, 4))
str(lst.4)
# List of 9
# $ : int [1:3] 1 2 3
# $ : int [1:3] 1 2 3
# $ : int [1:3] 4 5 6
# $ : int [1:3] 4 5 6
# $ : int [1:3] 4 5 6
# $ : int [1:3] 7 8 9
# $ : int [1:3] 7 8 9
# $ : int [1:3] 7 8 9
# $ : int [1:3] 7 8 9

Fügt man dieses Beispiel mit dem zusammen, was zu Beginn der Beschreibung von rep() über Schleifen gesagt wurde, so muss man im Schleifendurchlauf lediglich den Vektor erzeugen, der an das Argument times übergeben wird, das Erzeugen der Objekte wird anschließend durch rep() erledigt.

Bisher wurden Funktionen aus der Familie der apply()-Funktionen dadurch erläutert, dass ihre Arbeitsweise in die drei Phasen split-apply-combine unterteilt wurden. Versucht man dies hier anzuwenden, erkennt man:

1. Die split-Phase besteht lediglich in der Unterteilung der Liste x in ihre Komponenten.
2. Die apply-Phase trifft die Auswahl der Komponente gemäß den Eingabewerten times und each und erstellt die entsprechenden Kopien.
3. In der combine-Phase werden die Komponenten zu einem Objekt (meist Liste) zusammengefügt.

Die Funktion replicate()

Die Funktion replicate() besitzt keinen eigenen Eintrag in der R-Dokumentation, sie wird unter lapply() erklärt.

Betrachtet man nur die Eingabewerte von replicate(), könnte man meinen, dass wie mit rep() lediglich Kopien erzeugt werden:

• n: Anzahl der Wiederholungen.
• expr: Ein Ausdruck (expression), der bei jeder Wiederholung ausgewertet wird.
• simplify: falls möglich wird der Rückgabewert in ein Feld verwandelt (also zu Vektor oder Matrix im Fall einer oder zwei Dimensionen).

replicate(n, expr, simplify = "array")

Die Funktionalität von replicate() geht aber über die von rep() hinaus, wenn man mit expr Zufallszahlen erzeugen lässt.

Das folgende Skript zeigt zunächst ein typisches Beispiel für eine Simulation des Würfelns: Es wird 60 mal gewürfelt und anschließend der Mittelwert der Augenzahl berechnet.

v <- sample(x = (1:6), size = 60, prob = rep(1/6, times = 6), replace = TRUE)
mean(v)
# 3.416667

Wenn man die Einzelergebnisse, die in v abgespeichert werden, nicht benötigt, kann man die Funktionsaufrufe auch verschachteln:

mean( sample(x = (1:6), size = 60, prob = rep(1/6, times = 6), replace = TRUE) )

Möchte man diese Simulation 10 mal hintereinander ausführen, also 10 Mittelwerte erzeugen, so könnte man den letzten Befehl in einer for-Schleife aufrufen. Dies kann aber auch durch replicate() erledigt werden: mit dem Eingabewert n bestimmt man, wie oft ein Befehl ausgeführt wird, mit expr legt man fest, welcher Befehl ausgeführt wird.

replicate( n = 10, 
           expr = mean(sample(x = (1:6), size = 60, prob = rep(1/6, times = 6), replace = TRUE)) )
# [1] 3.783333 3.933333 3.250000 3.183333 3.200000 3.366667 3.400000 3.133333 3.500000 3.750000

Da der Rückgabewert einer einzelnen Simulation hier eine Zahl ist, bildet replicate() einen Vektor der Länge 10.

Ist dagegen der Rückgabewert von expr ein Vektor, gibt replicate() eine Matrix zurück. Das folgende Beispiel zeigt dies.

Die Simulation besteht jetzt darin, dass ein Würfel 6 mal geworfen wird. Mit Hilfe von replicate() wird dies 10 mal wiederholt. Der Rückgabewert von replicate() ist eine 6 × 10 Matrix, in der noch alle Würfel-Ergebnisse erkennbar sind:

replicate( n = 10, 
           expr = sample(x = (1:6), size = 6, prob = rep(1/6, times = 6), replace = TRUE) )
#      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,]    5    2    3    3    5    5    3    3    6     3
# [2,]    4    4    4    6    2    1    6    5    6     5
# [3,]    1    2    4    2    1    1    6    1    4     2
# [4,]    3    2    1    2    1    4    6    1    3     4
# [5,]    4    3    6    6    2    6    3    4    2     3
# [6,]    6    3    1    5    6    5    3    6    6     6

Versucht man wieder das Verhalten von replicate() mit den drei Phasen split-apply-combine zu beschreiben, so erhält man:

1. Es gibt keine split-Phase; die Funktion replicate() hat keinen Eingabewert für ein Objekt, das verarbeitet wird.
2. In der apply-Phase wird der an expr übergebene Ausdruck n-mal ausgewertet.
3. In der combine-Phase werden die Ergebnisse aus der apply-Phase zusammengesetzt. Man kann als Rückgabewert entweder ein Feld (also auch Vektor oder Matrix) oder eine Liste erzeugen, was über den Eingabewert simplify gesteuert wird.

Die Funktion ave()

Der Name ave() steht natürlich als Abkürzung für average und ave() wird in der default-Konfiguration verwendet, um Mittelwerte eines Vektors bezüglich eines Faktors zu berechnen. Dazu können auch mehrere Faktoren eingegeben werden; dann wird aus allen möglichen Level-Kombinationen ein neuer Faktor berechnet.

Die Funktion ave() befindet sich im Paket stats. Sie hat viele Gemeinsamkeiten – aber auch einige wichtige Unterschiede – mit der Funktion tapply(), die in Faktoren in R: Anwendungen erklärt wurde. Darauf wird weiter unten eingegangen.

Das Beispiel einer Temperatur-Messreihe, das hier ausführlich beschrieben wird, wurde schon in Faktoren in R: Anwendungen verwendet, um die Funktion tapply() einzuführen. Weiter unten wird das Beispiel fortgeführt, um die Funktion by() zu erklären.

Beschreibung der Funktion ave() und ihre Unterschiede zu tapply()

Die Eingabewerte der Funktion tapply() lauten:

tapply(X, INDEX, FUN = NULL, ..., default = NA, simplify = TRUE)

Die Funktion ave() besitzt drei Eingabewerte:

ave(x, ..., FUN = mean)

1. Ein numerischer Vektor x.
2. Das Argument ... , mit dem beliebig viele Faktoren eingegeben werden können.
3. Eine Funktion FUN, die per default gleich mean() ist.

Die Arbeitsweise von ave() kann man folgendermaßen beschreiben:

1. Der Vektor x wird gemäß der Faktoren in Gruppen eingeteilt (wobei auch einige Gruppen unbesetzt sein können).
2. Anschließend wird über die Komponenten von x iteriert:
   • Es wird festgestellt, zu welcher Gruppe x[i] gehört.
   • Auf die entsprechende Gruppe wird die Funktion FUN angewendet.
   • Dieser Wert wird in einem Vektor abgespeichert.
3. Der Rückgabewert ist dann dieser Vektor mit der Länge des Eingabevektors x (es wurde ja für jede Komponente von x ein Funktionswert berechnet).

Man beachte die Unterschiede zur Funktion tapply():

1. Der wichtigste Unterschied betrifft die Länge des Rückgabewertes: Bei tapply() ist es die Anzahl der Gruppen, bei ave() die Länge des Eingabevektors x. Zudem lässt sich bei tapply() konfigurieren, ob die Rückgabewerte in einem Feld (array) oder einer Liste verpackt werden sollen.
2. Bei tapply() gibt es keinen default-Wert für FUN.
3. Werden die Gruppen durch mehrere Faktoren gebildet, werden bei ave() die Faktoren im Argument ... der Reihe nach angegeben; bei tapply() ist das Argument ... für weitere Eingabewerte an FUN reserviert, daher müssen mehrere Faktoren im Argument INDEX zu einer Liste zusammengefasst werden.

Eine typische Anwendung von ave()

Das folgende Skript zeigt ein einfaches Beispiel, wie ave() eingesetzt werden kann. Dazu wird eine Messreihe für Temperaturen verwendet (Vektor temp in Zeile 2). Die Temperaturen werden in drei Intervalle eingeteilt (Zeile 5), dafür wird ein Faktor f.temp.coarse erzeugt (Zeile 10). Die Ausgabe mit table() zeigt, wie viele Messwerte sich in den drei Intervallen befinden (Zeile 15). Eine graphische Darstellung der Messreihe und der Intervallgrenzen ist in Abbildung 2 unten zu sehen.

# Temperaturen:
temp <- c(12.4, 13.8, 15.0, 18.2, 20.3, 22.8, 21.1, 21.7, 18.6, 16.4)

# Intervallgrenzen
brk <- c(-20, 15, 20, 40)

# Vereinbarung von Namen für die Temperatur-Intervalle
lev <- c("cold", "mild", "warm")

f.temp.coarse <- cut(x = temp, breaks = brk, labels = lev, ordered_result = TRUE)
f.temp.coarse
# [1] cold cold cold mild warm warm warm warm mild mild
# Levels: cold < mild < warm

table(f.temp.coarse)
# cold mild warm 
# 3    3    4 

print(ave(x = temp, f.temp.coarse), digits = 5)
# [1] 13.733 13.733 13.733 17.733 21.475 21.475 21.475 21.475 17.733 17.733

Beim Aufruf von ave() in Zeile 19 wird für x die Temperatur-Messreihe und für ... der Faktor f.temp.coarse eingegeben; FUN wird nicht gesetzt (der default-Wert ist mean). In der Ausgabe (Zeile 20) erkennt man 10 Komponenten, die aber nur 3 unterschiedliche Werte annehmen, nämlich die Mittelwerte der drei von f.temp.coarse gebildeten Gruppen. (Die print()-Funktion wird nur eingesetzt, um die Nachkommastellen zu kontrollieren.)

Ruft man anstelle von ave() die Funktion tapply() mit FUN = mean auf, so werden die drei Mittelwerte zurückgegeben und die Namen der Faktor-Levels werden weitergereicht:

means <- tapply(X = temp, INDEX = f.temp.coarse, FUN = mean)
means
# cold     mild     warm 
# 13.73333 17.73333 21.47500

str(means)
# num [1:3(1d)] 13.7 17.7 21.5
# - attr(*, "dimnames")=List of 1
# ..$ : chr [1:3] "cold" "mild" "warm"

Der Rückgabewert ist ein eindimensionales Feld, wie man an der Struktur erkennt (Zeile 6).

Aufgabe:

Definieren Sie zu obigem Beispiel einen weiteren Faktor und testen Sie den Aufruf von ave() mit 2 Faktoren in ... .

Beschreiben Sie die Unterschiede zum entsprechenden Aufruf von tapply().

Weitere Beispiele für die Anwendung von ave()

Im letzten Unterabschnitt wurden die Funktionen

ave(x = temp, f.temp.coarse, FUN = mean)
tapply(X = temp, INDEX = f.temp.coarse, FUN = mean)

verglichen. Dass man an den Mittelwerten der Temperaturen in den drei Intervallen interessiert ist, kann man sich leicht vorstellen. Aber es ist naheliegend zu fragen, wozu der Vektor dienen soll, der durch ave() erzeugt wird.

Dazu sollte man sich vorstellen, dass die Temperatur-Messreihe eigentlich einer Zuordnung

t ↦ T

entspricht, das soll heißen, dass die Temperatur T zu verschiedenen Zeitpunkten t gemessen wird und eine Funktion wie in Abbildung 2 aufgetragen werden kann. Und wenn man diese Funktion T(t) glätten möchte, so zeigt die Anwendung von ave() die Vorgehensweise: man definiert Intervalle, innerhalb derer der exakte Temperatur-Verlauf nicht interessiert und berechnet innerhalb dieser Intervalle die entsprechenden Temperatur-Mittelwerte. Es gibt natürlich ausgefeiltere Methoden zum Glätten von Funktionen, aber eine Vorstufe ist hier erkennbar. (Ein echtes Glätten findet mit ave() natürlich nicht statt, da zwischen den Intervallen große Sprünge auftreten.)

Und wenn man die Funktion T(t) nicht nur glätten, sondern zugleich nach oben abschätzen möchte, kann man anstelle von FUN = mean mit FUN = max arbeiten (Zeile 2). Oder mit FUN = min , wenn man die Funktion nach unten abschätzen möchte (Zeile 6). Das folgende Skript zeigt dies für die Temperatur-Messreihe.

# Abschätzung durch Maximum:
ave(x = temp, f.temp.coarse, FUN = max)
# [1] 15.0 15.0 15.0 18.6 22.8 22.8 22.8 22.8 18.6 18.6

# Abschätzung durch Minimum:
ave(x = temp, f.temp.coarse, FUN = min)
# [1] 12.4 12.4 12.4 16.4 20.3 20.3 20.3 20.3 16.4 16.4

# Berechnung der Spannweite der x-Werte:
ave(x = temp, f.temp.coarse, FUN = function(x){ return(max(x) - min(x)) })
# [1] 2.6 2.6 2.6 2.2 2.5 2.5 2.5 2.5 2.2 2.2

Das letzte Beispiel zeigt den Einsatz einer selbst definierten Funktion: jetzt wird für jeden x-Wert die Spannweite (oder Variation) des Intervalls berechnet, in dem er sich befindet (Zeile 10). Die tatsächlichen Variationen sind natürlich deutlich kleiner als die Intervall-Längen, mit denen der Faktor f.temp.coarse definiert wurde.

Die Implementierung von ave()

Da die Implementierung der Funktion ave() sehr lehrreich für den Umgang mit Vektoren, Faktoren und Funktionen ist, soll sie hier gezeigt und erläutert werden:

# Implementierung von ave(), aufrufbar mit print(ave):
function (x, ..., FUN = mean) 
{
    if (missing(...)) 
        x[] <- FUN(x)
    else {
        g <- interaction(...)
        split(x, g) <- lapply(split(x, g), FUN)
    }
    x
}

Zeile 1: Die Ausgabe des Quelltextes erfolgt mit der print()-Funktion für Funktionen.

Zeile 2: Man erkennt die Eingabewerte (wie sie oben besprochen wurden).

Zeile 4 und 5 (if-Zweig): Falls in ... keine Faktoren eingegeben wurden, wird der neue x-Wert direkt durch die Funktion FUN berechnet.

Zeile 6 bis 8 (else-Zweig): Die Faktoren aus ... werden an interaction() weitergereicht, wodurch derjenige Faktor g berechnet wird, der alle möglichen Level-Kombinationen der eingegebenen Faktoren enthält.

Schwer verständlich ist vielleicht Zeile 8: Auf der rechten Seite wird lapply() mit split(x, g) aufgerufen. Die Funktion split() sorgt hier dafür, dass der Vektor x gemäß dem Faktor g in Teilvektoren zerlegt wird. Durch lapply() wird jetzt über diese Teilvektoren iteriert und jedem Teilvektor ein Wert zugeordnet, der durch die Anwendung von FUN auf den Teilvektor entsteht (im default-Fall sind dies Mittelwerte).

Der Rückgabewert von lapply() ist eine Liste und daher werden die soeben berechneten Werte in eine Liste verpackt.

Auf der linken Seite von Zeile 8 steht die replacement-Version von split(). Jetzt wird für jede Komponente von x festgestellt, zu welcher Gruppe sie gehört (bezüglich des Faktors g) und der x-Wert wird durch den entsprechenden Wert aus der Liste ersetzt, die von lapply() erzeugt wurde.

Zeile 10: Der Rückgabewert von ave() ist dann der neu berechnete Wert von x.

Aufgabe:

Erstellen Sie eine "Testversion" von ave(), indem Sie obigen Quelltext mit debug-Ausgaben anreichern und testen Sie die oben besprochenen Beispiele.

Die Phasen split-apply-combine für ave()

Die Arbeitsweise der Funktion ave() wird – wie bei allen Funktionen der apply-Familie – sehr treffend durch die Aufspaltung in die drei Phasen split-apply-combine beschrieben; Abbildung 1 versucht dies für den Aufruf mit einem Vektor x, einem Faktor f und einer Funktion FUN darzustellen:

1. In der Phase split wird der Vektor x zuerst in seine Komponenten zerlegt; anschließend werden die Komponenten von x gemäß der Levels von f in Gruppen eingeteilt (in Abbildung 1 entspricht dies den beiden ersten Pfeilen, die mit f beschriftet sind).
2. In der Phase apply wird auf jede dieser Gruppen die Funktion FUN angewendet. Dies sollen in Abbildung 1 durch die Pfeile ganz rechts dargestellt werden, die mit FUN beschriftet sind. Allerdings ist so die Anzahl der Funktionswerte gleich der Anzahl der Levels von f. Die Funktion f ordnet aber jeder Komponente von x einen entsprechenden Funktionswert zu. Dies soll durch die Pfeile angedeutet werden, die nach unten zeigen.
3. In der Phase combine werden dann nur noch diese Funktionswerte zu einem Vektor zusammengefasst und bilden den Rückgabewert von ave().

Abbildung 1: Die drei Phasen split-apply-combine beim Aufruf von ave() mit einem Vektor x, einem Faktor f und der Funktion FUN. Weitere Erklärungen im Text.

Aufgabe: Versuchen Sie eine ähnliche Beschreibung der drei Phasen split-apply-combine für die Funktion tapply() zu geben.

Die Funktion by()

Beschreibung der Funktion by()

Die Funktion

tapply(X, INDEX, FUN = NULL, ..., default = NA, simplify = TRUE)

wird meist angewendet, um einen Vektor X mit Hilfe eines Faktors INDEX in Gruppen einzuteilen und dann eine Funktion FUN auf die einzelnen Gruppen anzuwenden. Die Funktion by() ist ein Wrapper für tapply(), um dies mit Dataframes auszuführen. Zum Verständnis der Arbeitsweise von by() sind daher Kenntnisse über Dataframes und Faktoren nötig, die hier nicht erklärt werden.

Die Eingabewerte von by() ähneln denen von tapply():

by(data, INDICES, FUN, ..., simplify = TRUE)

Dabei soll der Name data daran erinnern, dass by() üblicherweise auf Dataframes angewendet wird.

Eine typische Anwendung von by()

Im Folgenden wird an einem möglichst einfachen Beispiel gezeigt, wie by() arbeitet, wenn lediglich ein Faktor an das Argument INDEX übergeben wird. Im Allgemeinen kann man eine Liste von Faktoren übergeben. Dies wurde im Zusammenhang mit tapply() ausführlich erklärt und kann leicht auf by() übertragen werden.

Da Beispiel lehnt sich an das erste Beispiel für tapply() in Faktoren in R: Anwendungen an und besteht in folgenden Schritten:

1. Für eine Temperaturmessreihe wird ein Dataframe df.temp angelegt.
2. Die Temperaturen werden in drei Intervalle eingeteilt, dafür wird ein Faktor f.temp.coarse erzeugt.
3. Die Arbeitsweise von by() kann man am Besten nachvollziehen, wenn man für FUN die Funktion str() einsetzt.
4. Auswertung der Temperaturmessreihe: Anwendung von by() mit der Funktion summary().

1. Dataframe df.temp für eine Temperaturmessreihe anlegen

Zur Simulation einer Auswertung wird eine Temperatur-Messreihe angelegt, in der stündlich eine Temperatur gemessen wird. Dazu werden zwei Vektoren (für die Uhrzeit und die Temperatur) in einem Dataframe verpackt. Die Temperatur-Messreihe ist unten in Abbildung 2 links dargestellt.

temp <- c(12.4, 13.8, 15.0, 18.2, 20.3, 22.8, 21.1, 21.7, 18.6, 16.4)
df.temp <- data.frame(time = (9:18), temp = temp)

#    time temp
# 1     9 12.4
# 2    10 13.8
# 3    11 15.0
# 4    12 18.2
# 5    13 20.3
# 6    14 22.8
# 7    15 21.1
# 8    16 21.7
# 9    17 18.6
# 10   18 16.4

Abbildung 2: Links: Die Darstellung der Temperatur-Messreihe. Dazu wird stündlich die Temperatur aufgenommen. Modelliert wird dies durch ein Dataframe mit zwei Spalten (Zeit und Temperatur). Rechts: Die Temperaturen werden drei willkürlich gewählten Intervallen zugeordnet, was durch einen Faktor modelliert wird. Die vier Intervallgrenzen sind gelb bis rot gekennzeichnet.

Mit Hilfe von summary() kann man eine kleine statistische Auswertung des Dataframes erstellen lassen:

summary(df.temp)
#     time            temp      
# Min.   : 9.00   Min.   :12.40  
# 1st Qu.:11.25   1st Qu.:15.35  
# Median :13.50   Median :18.40  
# Mean   :13.50   Mean   :18.03  
# 3rd Qu.:15.75   3rd Qu.:20.90  
# Max.   :18.00   Max.   :22.80

2. Faktor erzeugen: Die Temperaturen werden in drei Intervalle eingeteilt

Um zu demonstrieren, wie die Temperatur-Messwerte und damit das Dataframe gruppiert werden, werden 4 willkürliche Intervallgrenzen gewählt (brk in Zeile 2), die drei Temperatur-Intervalle festlegen; ihre Namen lauten cold, mild und warm (lev in Zeile 5). Aus den Intervallen wird ein Faktor f.temp.coarse gebildet (mit Hilfe der Funktion cut() in Zeile 7).

# Intervallgrenzen
brk <- c(-20, 15, 20, 40)

# Vereinbarung von Namen für die Temperatur-Intervalle
lev <- c("cold", "mild", "warm")

f.temp.coarse <- cut(x = temp, breaks = brk, labels = lev, ordered_result = TRUE)
f.temp.coarse
# [1] cold cold cold mild warm warm warm warm mild mild
# Levels: cold < mild < warm

table(f.temp.coarse)
# cold mild warm 
# 3    3    4

Mit table() lässt sich leicht feststellen, wie viele Temperatur-Messwerte in den drei Intervallen liegen (Zeile 12).

In Abbildung 2 rechts ist nochmal die Temperatur-Messreihe zu sehen; aber jetzt ist die y-Achse so skaliert, dass alle Intervallgrenzen sichtbar sind.

3. Anwendung von str() auf die Teil-Dataframes

Wie so oft ist die Funktion str() am Besten geeignet, um Befehle nachzuvollziehen. Die Funktion by() wird auf das Dataframe angewendet, als Faktor wird obige Gruppierung der Temperaturen verwendet und als Funktion FUN die Struktur str():

by(data = df.temp, INDICES = f.temp.coarse, FUN = str)
# 'data.frame': 3 obs. of  2 variables:
#   $ time: int  9 10 11
# $ temp: num  12.4 13.8 15
# 'data.frame': 3 obs. of  2 variables:
#   $ time: int  12 17 18
# $ temp: num  18.2 18.6 16.4
# 'data.frame': 4 obs. of  2 variables:
#   $ time: int  13 14 15 16
# $ temp: num  20.3 22.8 21.1 21.7

Man erkennt, dass aus dem Dataframe df.temp drei Teil-Dataframes gebildet werden (mit 3, 3 beziehungsweise 4 Observablen, siehe Zeile 2, 5, 8). Und man kann leicht nachvollziehen, dass sie gerade den Temperatur-Intervallen entsprechen.

4. Auswertung der Messreihe mit der Funktion summary()

Übergibt man an das Argument FUN die Funktion summary(), werden für die drei Teil-Dataframes die statistischen Auswertungen vorgenommen, die oben für das gesamte Dataframe erzeugt wurden:

temp.summary <- by(data = df.temp, INDICES = f.temp.coarse, FUN = summary)
temp.summary
# f.temp.coarse: cold
# time           temp      
# Min.   : 9.0   Min.   :12.40  
# 1st Qu.: 9.5   1st Qu.:13.10  
# Median :10.0   Median :13.80  
# Mean   :10.0   Mean   :13.73  
# 3rd Qu.:10.5   3rd Qu.:14.40  
# Max.   :11.0   Max.   :15.00  
# ------------------------------------------------------------------------------------------------------------------------

#   f.temp.coarse: mild
# time            temp      
# Min.   :12.00   Min.   :16.40  
# 1st Qu.:14.50   1st Qu.:17.30  
# Median :17.00   Median :18.20  
# Mean   :15.67   Mean   :17.73  
# 3rd Qu.:17.50   3rd Qu.:18.40  
# Max.   :18.00   Max.   :18.60  
# ------------------------------------------------------------------------------------------------------------------------ 

#   f.temp.coarse: warm
# time            temp      
# Min.   :13.00   Min.   :20.30  
# 1st Qu.:13.75   1st Qu.:20.90  
# Median :14.50   Median :21.40  
# Mean   :14.50   Mean   :21.48  
# 3rd Qu.:15.25   3rd Qu.:21.98  
# Max.   :16.00   Max.   :22.80

Der Rückgabewert von by()

Welchen Rückgabewert by() erzeugt, hängt von zwei Eingabewerten ab:

1. Vom Rückgabewert der Funktion FUN.
2. Davon, ob simplify mit dem default-Wert TRUE verwendet wird oder gleich FALSE gesetzt wird. Im ersten Fall wird versucht (wie bei tapply()) ein Feld zu erzeugen. Im zweiten Fall wird eine Liste erzeugt.

Untersucht man etwa die Struktur des Objektes, das durch by(data = df.temp, INDICES = f.temp.coarse, FUN = str) erzeugt wird, erhält man:

str(by(data = df.temp, INDICES = f.temp.coarse, FUN = str))
# List of 3
# $ cold: NULL
# $ mild: NULL
# $ warm: NULL
# - attr(*, "dim")= int 3
# - attr(*, "dimnames")=List of 1
# ..$ f.temp.coarse: chr [1:3] "cold" "mild" "warm"
# - attr(*, "call")= language by.data.frame(data = df.temp, INDICES = f.temp.coarse, FUN = str)
# - attr(*, "class")= chr "by"

Daran erkennt man:

1. Da die Funktion str() keinen Rückgabewert hat, entsteht eine Liste, deren 3 Komponenten gleich NULL sind.
2. Es sind weitere Attribute gesetzt, die von den Eingabewerten data und INDICES stammen.
3. Zusätzlich ist das Attribut class gesetzt und besitzt den Wert "by" . Dadurch ist das Objekt der Klasse by zugehörig.

Die Phasen split-apply-combine für by()

Vergleicht man die Arbeitsweise von by() mit Abbildung 1 oben, wo die Phasen split-apply-combine für ave() dargestellt sind, so findet man nur wenige Unterschiede:

1. Verarbeitet wird anstelle eines Vektors ein Dataframe. Durch den eingegebenen Faktor wird es in Teil-Dataframes zerlegt (anstelle der Teilvektoren).
2. Die Funktion FUN wird auf diese Teil-Dataframes angewendet.
3. Wird der Rückgabewert gebildet, so wird für jedes Teil-Dataframe eine Listen-Komponente erzeugt; die Funktion ave() hatte für jede Komponente des Eingabevektors x einen Funktionswert berechnet und diese zu einem Vektor zusammengefasst.
4. In der Liste sind weitere Attribute gesetzt, die im letzten Unterabschnitt gezeigt wurden.

Die Funktion aggregate() für Dataframes

Die Funktion aggragate() befindet sich im Paket stats und ist generisch implementiert. Hier wird allerdings nur die Version für Dataframes besprochen:

## S3 method for class 'data.frame'
aggregate(x, by, FUN, ..., simplify = TRUE, drop = TRUE)

Sie ist sehr eng verwandt mit tapply() und wendet eine Funktion FUN auf Daten x an, die mit Hilfe von by, einer Liste von Faktoren, gruppiert werden. Der Unterschied zu tapply() besteht im Rückgabewert wie in den Beispielen unten besprochen wird.

Die genaue Bedeutung der Eingabewerte ist:

• x ist ein Dataframe oder wird in ein Dataframe verwandelt.
• by ist eine Liste von Vektoren (oft Spalten eines Dataframes) oder Faktoren, mit denen die Daten x gruppiert werden. Ist by kein Faktor, wird es in einen Faktor verwandelt. Auch wenn nur ein Objekt eingegeben wird, muss es in eine Liste verpackt werden. (Da das Argument ... hierfür nicht zur Verfügung steht, kann man mehrere Faktoren nicht einfach der Reihe nach eingeben, sondern muss sie in eine Liste verpacken.) Die Komponenten der Liste by müssen dieselbe Länge besitzen wie x.
• Die Funktion FUN, die auf die gruppierten Daten angewendet wird.
• Weitere Eingabewerte für FUN können in ... eingegeben werden.
• Mit simplify wird festgelegt, ob der Rückgabewert ein Dataframe ist (default-Wert) oder eine Liste (simplify = FALSE ).
• Der logische Wert drop gibt an, ob nicht besetzte Gruppen aus dem Rückgabewert entfernt werden (siehe Beispiel unten).

Der Rückgabewert ist ein nach speziellen Regeln aufgebautes Dataframe, in dem sowohl alle möglichen Level-Kombinationen als auch die mit Hilfe von FUN berechneten Werte enthalten sind. Wie das Dataframe aufgebaut ist, wird unten bei den Beispielen besprochen.

Ein einfaches Beispiel für die Anwendung von aggregate()

Um zunächst ein möglichst einfaches Beispiel für die Anwendung von aggragate() zu zeigen, wird wieder die Temperatur-Messreihe von oben definiert:

# Dataframe:
temp <- c(12.4, 13.8, 15.0, 18.2, 20.3, 22.8, 21.1, 21.7, 18.6, 16.4)
df.temp <- data.frame(time = (9:18), temp = temp)
df.temp
#    time temp
# 1     9 12.4
# 2    10 13.8
# 3    11 15.0
# 4    12 18.2
# 5    13 20.3
# 6    14 22.8
# 7    15 21.1
# 8    16 21.7
# 9    17 18.6
# 10   18 16.4

# Faktor:
# Intervallgrenzen
brk <- c(-20, 15, 20, 40)

# Vereinbarung von Namen für die Temperatur-Intervalle
lev <- c("cold", "mild", "warm")

f.temp.coarse <- cut(x = temp, breaks = brk, labels = lev, ordered_result = TRUE)
f.temp.coarse
# [1] cold cold cold mild warm warm warm warm mild mild
# Levels: cold < mild < warm

Dieses Dataframe df.temp und der Faktor f.temp.coarse sind bereits aus früheren Beispielen bekannt. Jetzt wird aggragate() aufgerufen, wobei das Dataframe, der Faktor f.temp.coarse und die Funktion mean als Eingabewerte gesetzt werden; für die Ausgabe wird print() eingesetzt, um die Nachkommastellen zu kontrollieren:

df.aggr <- aggregate(x = df.temp, by = list(f.temp.coarse), FUN = mean)
print(df.aggr, digits = 4)
#   Group.1  time  temp
# 1    cold 10.00 13.73
# 2    mild 15.67 17.73
# 3    warm 14.50 21.48

str(df.aggr)
# 'data.frame': 3 obs. of  3 variables:
#   $ Group.1: Ord.factor w/ 3 levels "cold"<"mild"<..: 1 2 3
# $ time   : num  10 15.7 14.5
# $ temp   : num  13.7 17.7 21.5

Man erkennt an den Ausgaben:

1. Das ursprüngliche Dataframe df.temp hatte 10 Zeilen (für die 10 Zeiten beziehungsweise Temperaturen), das durch aggragate() erzeugte Dataframe df.aggr besitzt nur noch drei Zeilen. Dies entspricht genau der Anzahl der Levels des Faktors f.temp.coarse, der an das Argument by übergeben wurde.
2. Das Dataframe df.aggr hat 3 Spalten. Neu hinzugekommen ist die Spalte Group.1, die die drei Levels des Faktors beinhaltet. Werden mehrere Faktoren mit by gesetzt, werden hier alle Level-Kombinationen gebildet, siehe Beispiel weiter unten. An der Ausgabe der Struktur liest man ab, dass die erste Spalte kein Vektor sondern ein Faktor ist.
3. Die Namen der Spalten time und temp werden aus dem ursprünglichen Dataframe übernommen; sie werden jetzt mit den Mittelwerten der Gruppen belegt (die Zahlenwerte sollten aus früheren Beispielen bekannt sein).

An diesem Beispiel erkennt man somit den wichtigsten Unterschied zwischen aggragate() und den anderen Funktionen, mit denen auch schon die Mittelwerte der Gruppen berechnet wurden. Der Rückgabewert von aggragate() ist ein Dataframe und zwar ein spezielles Dataframe: Es wird eine neue Spalte gebildet, in der die Levels des verwendeten Faktors stehen und somit ist es leicht die Zuordnung zwischen den Levels und den mit FUN berechneten Werte herzustellen. Für eine weitere Verarbeitung der Daten kann dies vorteilhaft sein gegenüber einer Berechnung, die nur die Mittelwerte liefert (beziehungsweise die von FUN berechneten Werte).

Vergleich von aggregate() mit anderen Funktionen

Um den zuletzt genannten Vorteil von aggregate() zu verdeutlichen, sollen nochmal im Überblick die anderen Funktionsaufrufe für das Beispiel zur Mittelwert-Berechnung bei der Temperaturmessreihe gezeigt werden; dazu wird jeweils der Rückgabewert ausgegeben (mit einer Beschränkung der Nachkommastellen) und die Struktur str() des Rückgabewertes:

# temp und df.temp wie oben

############################################
#      aggregate():
df.aggr <- aggregate(x = df.temp, by = list(f.temp.coarse), FUN = mean)
print(df.aggr, digits = 4)
#   Group.1  time  temp
# 1    cold 10.00 13.73
# 2    mild 15.67 17.73
# 3    warm 14.50 21.48

str(df.aggr)
# 'data.frame': 3 obs. of  3 variables:
#   $ Group.1: Ord.factor w/ 3 levels "cold"<"mild"<..: 1 2 3
# $ time   : num  10 15.7 14.5
# $ temp   : num  13.7 17.7 21.5

############################################
#      tapply():
mw.tapply <- tapply(X = temp, INDEX = f.temp.coarse, FUN = mean)
print(mw.tapply, digits = 4)
#  cold  mild  warm 
# 13.73 17.73 21.48 

str(mw.tapply)
# num [1:3(1d)] 13.7 17.7 21.5
# - attr(*, "dimnames")=List of 1
# ..$ : chr [1:3] "cold" "mild" "warm"

############################################
#      ave():
mw.ave <- ave(x = temp, f.temp.coarse)
print(mw.ave, digits = 4)
#  [1] 13.73 13.73 13.73 17.73 21.48 21.48 21.48 21.48 17.73 17.73

str(mw.ave)
#  num [1:10] 13.7 13.7 13.7 17.7 21.5 ...

############################################
#      by():
mw.by <- by(data = temp, INDICES = f.temp.coarse, FUN = mean)
print(mw.by, digits = 4)
# f.temp.coarse: cold
# [1] 13.73
# ------------------------------------------------------------------------------------------------------------------------ 
#   f.temp.coarse: mild
# [1] 17.73
# ------------------------------------------------------------------------------------------------------------------------ 
#   f.temp.coarse: warm
# [1] 21.48

str(mw.by)
# 'by' num [1:3(1d)] 13.7 17.7 21.5
# - attr(*, "dimnames")=List of 1
# ..$ f.temp.coarse: chr [1:3] "cold" "mild" "warm"
# - attr(*, "call")= language by.default(data = temp, INDICES = f.temp.coarse, FUN = mean)

Die folgende Tabelle gibt die Übersicht, welche Rückgabewerte erzeugt werden:

Aufgaben:

1. Die Funktionen aggragate(), tapply() und by() besitzen das Argument simplify = TRUE . Untersuchen Sie das Verhalten der Funktionen im Beispiel oben, wenn man simplify = FALSE setzt.

2. Geben Sie zu obiger Tabelle an, welche Länge (beziehungsweise Dimensionen) die Rückgabewerte haben.

Die Funktion aggregate() bei mehreren Faktoren

Um die Arbeitsweise und den Rückgabewert von aggragate() kennenzulernen und um sie mit anderen Funktionen zu vergleichen, wurde zunächst ein einfaches Beispiel mit nur einem Faktor gewählt. Wie das Dataframe für den Rückgabewert im Allgemeinen aufgebaut wird, erkennt man erst, wenn man mehrere Faktoren verwendet. Dazu wird das Beispiel mit der Turniertabelle aufgegriffen, mit dem in Faktoren in R: Anwendungen die Funktion tapply() erklärt wurde.

Das folgende Skript definiert das Dataframe result, in dem alle relevanten Daten über das Turnierergebnis gespeichert sind; ein Spieler ist eindeutig durch seine Punktzahl identifizierbar, zusätzlich ist angegeben, für welchen Verein er angetreten ist und zu welcher Altersklasse er gehört:

# Definition der Turnier-Tabelle als Dataframe:

score <- c(28, 27, 27, 24, 23, 23, 21, 18, 18, 18, 17 , 17, 17, 16, 14, 14, 13, 13, 13, 11, 11, 10, 9, 9, 8, 8, 4, 3, 1, 0)
club <- c("A", "C", "G", "B", "A", "B", "D", "A", "F", "E", "C", "D", "B", "G", "B", "C", "E", "D", "A", "B", "G", "E", "F", "F", "D", "C", "F", "A", "B", "G")
age <- c("J", "S", "S", "S", "S", "J", "S", "J", "J", "S", "J", "S", "S", "J", "S", "J", "J", "S", "J", "S", "J", "S", "S", "S", "J", "S", "S", "J", "S", "S")

length(score)           # 30
length(club)            # 30
length(age)             # 30
sum(score)              # 435 = 15*29 = 29 + 28 + ... + 1

result <- data.frame(score, club, age)

str(result)
# 'data.frame': 30 obs. of  3 variables:
#   $ score: num  28 27 27 24 23 23 21 18 18 18 ...
# $ club : Factor w/ 7 levels "A","B","C","D",..: 1 3 7 2 1 2 4 1 6 5 ...
# $ age  : Factor w/ 2 levels "J","S": 1 2 2 2 2 1 2 1 1 2 ...

Es gibt 7 Vereine und 2 Altersklassen, womit sich 14 Kombinationen (oder Gruppen) bilden lassen. Es ist naheliegend zu fragen:

Wie viele Punkte (in der Summe oder im Mittel) haben die Vertreter der 14 Gruppen erzielt?

Je nachdem, welchen Rückgabewert man benötigt, kann man für die Beantwortung der Frage eine der oben besprochenen Funktionen auswählen. Mit tapply() wurde die Frage in Faktoren in R: Anwendungen bereits beantwortet:

p.age.club <- tapply(X = result$score, INDEX = list(result$age, result$club), FUN = sum)
p.age.club
#    A  B  C  D  E  F  G
# J 62 23 31  8 13 18 27
# S 23 67 35 51 28 22 27

str(p.age.club)
# num [1:2, 1:7] 62 23 23 67 31 35 8 51 13 28 ...
# - attr(*, "dimnames")=List of 2
# ..$ : chr [1:2] "J" "S"
# ..$ : chr [1:7] "A" "B" "C" "D" ...

Berechnet man diese Punkte-Summen mit by() muss man nur die Namen der Argumente anpassen und erhält eine andere Form der Ausgabe.

Interessanter ist hier die Anwendung von aggragate() (wieder mit FUN = sum ):

df.age.club <- aggregate(x = result$score, by = list(result$age, result$club), FUN = sum)
df.age.club
#    Group.1 Group.2  x
# 1        J       A 62
# 2        S       A 23
# 3        J       B 23
# 4        S       B 67
# 5        J       C 31
# 6        S       C 35
# 7        J       D  8
# 8        S       D 51
# 9        J       E 13
# 10       S       E 28
# 11       J       F 18
# 12       S       F 22
# 13       J       G 27
# 14       S       G 27

str(df.age.club)
# 'data.frame': 14 obs. of  3 variables:
# $ Group.1: Factor w/ 2 levels "J","S": 1 2 1 2 1 2 1 2 1 2 ...
# $ Group.2: Factor w/ 7 levels "A","B","C","D",..: 1 1 2 2 3 3 4 4 5 5 ...
# $ x      : num  62 23 23 67 31 35 8 51 13 28 ...

Zeile 1: Die beiden Faktoren, nach denen die berechneten Summen gruppiert werden sollen, müssen zu einer Liste zusammengefasst werden.

Zeile 19 und 20: An der Struktur des Rückgabewertes von aggragate() erkennt man, dass ein Dataframe mit 14 Zeilen erzeugt wurde (es gibt hier 14 Level-Kombinationen).

Zeile 3 bis 17: Die Ausgabe des Dataframes zeigt zwei führende Spalten Group.1 und Group.2, die in dem ursprünglichen Dataframe nicht enthalten waren. Sie beschreiben alle 14 möglichen Level-Kombinationen der verwendeten Faktoren, also alle Kombinationen von Altersklasse und Vereinszugehörigkeit. Die dritte Spalte des Dataframes mit Namen x beinhaltet die berechneten Punkte-Summen der 14 Gruppen.

Nach diesem Beispiel sollte auch klar sein, welche Bedeutung der Eingabewert drop = TRUE von aggragate() hat: Sind Level-Kombinationen nicht besetzt, werden sie nicht in den Rückgabewert aufgenommen. Da hier alle Kombinationen besetzt sind, erhält man tatsächlich 14 Zeilen.

Und es sollte auch klar sein: Werden n Faktoren dem Argument by übergeben, werden im Rückgabewert n Spalten erzeugt. (Sie werden mit Hilfe von expand.grid() aus den Faktoren erzeugt.)

Um die Ergebnisse leichter zu interpretieren, kann man das Dataframe nach der Spalte x sortieren:

df.age.club.sort <- df.age.club[order(df.age.club$x), ]
df.age.club.sort
#    Group.1 Group.2  x
# 7        J       D  8
# 9        J       E 13
# 11       J       F 18
# 12       S       F 22
# 2        S       A 23
# 3        J       B 23
# 13       J       G 27
# 14       S       G 27
# 10       S       E 28
# 5        J       C 31
# 6        S       C 35
# 8        S       D 51
# 1        J       A 62
# 4        S       B 67

Zur Veranschaulichung bietet sich die Funktion interaction.plot() aus dem Paket stats an, die hier nur kurz erklärt werden soll:

interaction.plot(x.factor = result$club, trace.factor = result$age, 
                 response = result$score, fun = sum,
                 col = c("blue", "red"), lty = 3,
                 main = "Punkte je Altersklasse und Verein",
                 xlab = "Verein", ylab = "Summe der Punkte", trace.label = "Alter")

Die wichtigsten Argumente sind:

1. Der Faktor x.factor, der auf der x-Achse aufgetragen wird.
2. Der Faktor trace.factor, der verwendet wird, um die beiden Graphen zu zeichnen.
3. Das Argument response gibt an, aus welchen Daten der Wert auf der y-Achse berechnet wird.
4. Die Funktion fun gibt an, mit welcher Funktion die y-Werte berechnet werden. Da es 14 Kombinationen der Faktoren-Level gibt, werden insgesamt 14 y-Werte berechnet. Sie werden zu zwei "Spuren" für die beiden Altersklassen zusammengefasst.

Die restlichen Argumente dienen der Beschriftung und Darstellung.

Abbildung 3: Veranschaulichung des mit Hilfe von aggregate() berechneten Dataframes durch die Funktion interaction.plot(). Aufgetragen sind die Punkte-Summen, die von den 14 Gruppen (Vereinszugehörigkeit und Altersklasse) erzielt werden. Die 14 Level-Kombinationen werden in 2 Spuren für die beiden Altersklassen dargestellt.

Aufgabe:

Berechnen Sie das entsprechende Dataframe, das nicht die Punkte-Summen der 14 Gruppen sondern ihre mittlere Punktzahl angibt.

Erzeugen Sie das entsprechende Diagramm.

Die Phasen split-apply-combine für die Funktion aggregate()

Wie bei allen Funktionen der apply-Familie kann man ihre Arbeitsweise treffend durch die Aufteilung in die drei Phasen split-apply-combine beschreiben. Für aggragate() lautet die Beschreibung der Phasen:

1. In der split-Phase wird der Eingabewert x – falls es noch keines ist – in ein Dataframe x verwandelt und gemäß der Faktoren in by in Teil-Dataframes aufgeteilt.
2. Die apply-Phase wird die Funktion FUN auf die Teil-Dataframes angewendet. Falls drop = TRUE , werden nicht besetzte Level-Kombinationen ignoriert und auch nicht in den Rückgabewert aufgenommen.
3. In der combine-Phase wird ein neues Dataframe aufgebaut: Wenn by aus n Faktoren besteht, enthält es gegenüber x, dem eingegebenen Dataframe, n zusätzliche Spalten, die alle möglichen Level-Kombinationen aus by wiedergeben. Die weiteren Spalten enthalten die in der apply-Phase berechneten Werte. Die Anzahl der Zeilen stimmt mit der Anzahl der Teil-Dataframes überein. Ist simplify = FALSE gesetzt, wird anstelle eines Dataframes eine Liste erzeugt.

Anwendung: Berechnung der Faltung zweier diskreter Wahrscheinlichkeitsverteilungen

Was ist eine Faltung?

Es ist nicht das Ziel hier eine mathematisch exakte Definition der Faltung zu formulieren und ihre Eigenschaften zu besprechen. Es soll lediglich an einem einfachen Beispiel gezeigt werden, wie die Faltung motiviert ist und wie man sie in "einfachen Fällen" berechnet. Ausführliche Erläuterungen zur Faltung finden sich in Einführung des Begriffs der Faltung von Wahrscheinlichkeitsmaßen.

Mit "einfachen Fällen" ist hier folgendes gemeint:

Es werden Zufallsvariablen X betrachtet, die endlich viele Werte x₁, x₂, ..., x_n annehmen können. Die Wahrscheinlichkeiten

p_i = P(X = x_i), i = 1, 2, ..., n

müssen sich dann zu 1 addieren. Etwas salopp werden die Menge der p_i als die Wahrscheinlichkeitsverteilung der Zufallsvariable X bezeichnet.

Nicht behandelt werden die Fälle, dass eine Zufallsvariable abzählbar viele Werte oder sogar überabzählbar viele Werte annimmt. Bei abzählbar vielen Werten muss die Folge der Wahrscheinlichkeiten p_i eine Nullfolge sein, deren Summe immer noch gegen 1 konvergiert. Bei überabzählbar vielen Werten muss man die Wahrscheinlichkeit eines Ereignisses der Art

a < X < b, mit a < b,

mit Hilfe einer sogenannten Wahrscheinlichkeitsdichte beschreiben.

Wird ein Zufallsexperiment zweimal hintereinander und zwar unabhängig voneinander durchgeführt, so ist oft die Summe der Zufallsvariablen, die die Einzel-Experimente beschreiben, eine relevante Größe. Paradebeispiel ist das zweimalige Werfen eines Würfels, wobei die Summe der Augenzahlen gesucht ist.

Ist X die Zufallsvariable, die die Augenzahl beim ersten Wurf beschreibt, und Y die Zufallsvariable für den zweiten Wurf, so beschreibt die Zufallsvariable Z mit

Z = X + Y

die Summe der Augenzahlen.

In diesem Fall sagt man auch, dass X und Y identisch verteilt sind, da für beide Würfe der identische Würfel verwendet wird. Und X und Y sind unabhängig voneinander, wenn man annehmen kann, dass der erste Wurf das Ergebnis des zweiten Wurfes nicht beeinflusst (und umgekehrt).

Eine mathematisch naheliegende Frage ist jetzt:

Wie kann man aus den Wahrscheinlichkeitsverteilungen von X und Y die Wahrscheinlichkeitsverteilung von Z = X + Y herleiten?

Für das zweimalige Würfeln ist die Rechnung ganz einfach, wenn man einen Laplace-Würfel verwendet:

• Jede Kombination (x, y) kommt mit Wahrscheinlichkeit 1/36 vor.
• Zu einem gegebenen Wert von z der Zufallsvariable Z muss man nur abzählen, wie viele Kombinationen (x, y) es gibt mit x + y = z.
• Die Anzahl dieser Kombinationen multipliziert mit 1/36 ist die gesuchte Wahrscheinlichkeit P(Z = z).

Und jetzt muss man nur z durch alle möglichen Werte laufen lassen:

• Ist z = 2, so gibt es nur eine Kombination (x, y) = (1, 1), also ist P(Z = 2) = 1/36.
• Für z = 3 gibt es 2 Kombinationen: (1, 2), (2, 1), also ist P(Z = 3) = 2/36 = 1/18.
• ...
• Für z = 7 gibt es 6 Kombinationen: (1, 6), (2, 5), ..., (6, 1), also ist P(Z = 7) = 6/36 = 1/6.
• ...
• Für z = 12 gibt es wieder nur eine Kombination (x, y) = (6, 6), also ist P(Z = 12) = 1/36.

Etwas aufgebläht schreibt man die Berechnung der Wahrscheinlichkeit, etwa für P(Z = 7):

P(Z = 7) = P(X = 1) · P(Y = 6) + P(X = 2) · P(Y = 5) + ... + P(X = 6) · P(Y = 1).

Aber das Beispiel von P(Z = 7) zeigt wie man im Allgemeinen vorgehen muss, um die Wahrscheinlichkeiten der Werte von Z zu berechnen:

• Man gibt sich einen Wert für z vor.
• Jetzt durchläuft x alle möglichen Werte (also die Wertemenge von X).
• Damit Z den gegebenen Wert z annimmt, muss y den Wert y = z - x annehmen, denn nur so ist x + y = x + (z - x) = z.
• Die Wahrscheinlichkeit P(Z = z) berechnet sich jetzt aus einer Summe von Produkten, nämlich:

P(Z = z) = ∑_i P(X = x_i) · P(Y = z - x_i).

Führt man diese Berechnung für alle z aus der Wertemenge von Z durch, hat man die Wahrscheinlichkeitsverteilung von Z berechnet.

Die Begründung für die Berechnung von P(Z = z) kann man leicht angeben: Die Ereignisse, die in den einzelnen Summanden vorkommen, sind disjunkt, daher können ihre Wahrscheinlichkeiten addiert werden. Und die einzelnen Summanden beschreiben jeweils Verbund-Wahrscheinlichkeiten der Art P(X = x_i, Y = z - x_i), die aufgrund der Unabhängigkeit als Produkt P(X = x_i) · P(Y = z - x_i) geschrieben werden.

Derartige Summen von Produkten nennt man die Faltung der Wahrscheinlichkeitsverteilungen von X und Y.

Die folgenden Abschnitte entwickeln eine Funktion, die die Faltung für diskrete Zufallsvariablen mit endlich vielen Werten berechnet; sie beruht darauf, dass eine Zufallsvariable durch ein Dataframe modelliert wird.

Modellierung einer Zufallsvariable mit Hilfe eines Dataframes

In Grundbegriffe der Wahrscheinlichkeitsrechnung: Die Zufallsvariable wurde diskutiert, wie man Zufallsvariablen mit Hilfe von Dataframes modellieren kann. Dazu werden die möglichen Werte der Zufallsvariable in einer Spalte namens value und die zugehörigen Wahrscheinlichkeiten in einer Spalte namens prob abgespeichert.

Für den Fall eines Laplace-Würfels definiert man das Dataframe durch:

X <- data.frame(value = (1:6), prob = rep(x = 1/6, times = 6))
X
#   value      prob
# 1     1 0.1666667
# 2     2 0.1666667
# 3     3 0.1666667
# 4     4 0.1666667
# 5     5 0.1666667
# 6     6 0.1666667

In Grundbegriffe der Wahrscheinlichkeitsrechnung: Die Zufallsvariable wurde dann auch diskutiert, welche Anforderungen man an ein Dataframe stellen wird, damit es eine Zufallsvariable beschreibt. Dort wurde eine Funktion is.randomVariable() implementiert, die ein gegebenes Dataframe untersucht – es lässt sich lange diskutieren, ob dies alle Anforderungen an eine Zufallsvariable sind. Auch die Fehlermeldungen in den cat()-Befehlen lassen sich aussagekräftiger formulieren.

is.randomVariable <- function(rv, eps = 1e-8){
  stopifnot(is.data.frame(rv))
  result <- TRUE
  
  if( !( identical(names(rv), c("prob", "value")) || identical(names(rv), c("value", "prob")) ) ){
    cat("names \n")
    result <- FALSE
  }
  if( !identical(length(rv$value), length(rv$prob)) ){
    cat("length != \n")
    result <- FALSE
  }
  if(any(rv$prob < 0)){
    cat("neg. Prob. \n")
    result <- FALSE
  }
  if( abs( sum(rv$prob) - 1 ) > eps ){
    cat("sum Prob \n")
    result <- FALSE
  }
  if( !identical(rv$value, unique(rv$value)) ){
    cat("Eind. \n")
    result <- FALSE
  }
  return(result)
}

Die Berechnung der Faltung

Am Paradebeispiel des zweimaligen Würfelns kann man sich jetzt leicht überlegen, wie die Faltung zweier Wahrscheinlichkeitsverteilungen berechnet wird.

Der Laplace-Würfel wird wie oben gezeigt als Dataframe modelliert; X steht für die Zufallsvariable, die die Augenzahl beschreibt. Ebenso kann eine zweite Zufallsvariable Y modelliert werden, die hier aber mit X übereinstimmt.

Betrachtet man jetzt die Zufallsvariable Z, die durch

Z = X + Y

definiert ist, so bildet man zunächst alle möglichen Kombinationen der Werte von X mit den Werten von Y. Dies kann mit expand.grid() oder outer() geschehen. Hier wird outer() verwendet.

Die Wahrscheinlichkeit einer Kombination wird durch das entsprechende Produkt der Einzelwahrscheinlichkeiten berechnet, was wiederum mit outer() realisiert wird. Da outer() eigentlich eine Matrix erzeugt, wird in beiden Fällen as.vector() angewendet.

Das folgende Skript zeigt wie eine Vorstufe der Zufallsvariable Z berechnet wird:

Z <- data.frame( value = as.vector( outer(X = X$value, Y = X$value, FUN = "+") ),
                 prob = as.vector( outer(X = X$prob, Y = X$prob, FUN = "*") ) )
Z
#    value       prob
# 1      2 0.02777778
# 2      3 0.02777778
# 3      4 0.02777778
# 4      5 0.02777778
# 5      6 0.02777778
# 6      7 0.02777778
# 7      3 0.02777778
# 8      4 0.02777778
# 9      5 0.02777778
# 10     6 0.02777778
# 11     7 0.02777778
# 12     8 0.02777778
# 13     4 0.02777778
# 14     5 0.02777778
# 15     6 0.02777778
# 16     7 0.02777778
# 17     8 0.02777778
# 18     9 0.02777778
# 19     5 0.02777778
# 20     6 0.02777778
# 21     7 0.02777778
# 22     8 0.02777778
# 23     9 0.02777778
# 24    10 0.02777778
# 25     6 0.02777778
# 26     7 0.02777778
# 27     8 0.02777778
# 28     9 0.02777778
# 29    10 0.02777778
# 30    11 0.02777778
# 31     7 0.02777778
# 32     8 0.02777778
# 33     9 0.02777778
# 34    10 0.02777778
# 35    11 0.02777778
# 36    12 0.02777778

Betrachtet man die Werte (Spalte value), die Z annehmen kann, fällt sofort auf, dass nahezu alle Werte mehrfach vorkommen (außer 2 und 12), was bei der Modellierung einer Zufallsvariable als Dataframe nicht vorkommen darf. Daher wird Z im folgenden Sinn bereinigt: mehrfach vorkommende Werte werden entfernt, aber dabei die zugehörigen Wahrscheinlichkeiten addiert – dies entspricht genau der Summenbildung bei der Faltung:

Z.aggr <- aggregate(x = Z$prob, by = list(Z$value), FUN = sum)
names(Z.aggr) <- c("value", "prob")
Z.aggr
#    value       prob
# 1      2 0.02777778
# 2      3 0.05555556
# 3      4 0.08333333
# 4      5 0.11111111
# 5      6 0.13888889
# 6      7 0.16666667
# 7      8 0.13888889
# 8      9 0.11111111
# 9     10 0.08333333
# 10    11 0.05555556
# 11    12 0.02777778

Z.aggr$prob * 36
# [1] 1 2 3 4 5 6 5 4 3 2 1
is.randomVariable(Z.aggr)    # TRUE
sum(Z.aggr$x)    # 1

Zeile 1: Die Werte der "Vorstufe" von Z werden als Faktor verwendet (Argument by), um mehrfach vorkommende Werte zusammenzufassen. Dazu werden die Wahrscheinlichkeiten addiert durch x = Z$prob und FUN = sum in aggragate().

Zeile 2: Da das von aggragate() erzeugte Dataframe die falschen Namen hat, werden sie gesetzt, um wieder eine Zufallsvariable zu modellieren.

Zeile 17: Die Wahrscheinlichkeiten sind leichter zu lesen, wenn man sie mit 36 multipliziert.

Zeile 19 und 20: Man kann testen, ob das Dataframe Z.aggr die Bedingungen an eine Zufallsvariable erfüllt.

Implementierung der Faltung als Funktion

Nach diesen Vorüberlegungen ist es nicht mehr schwer eine Funktion composition() zu implementieren, die zu zwei als Dataframe gegebenen Zufallsvariablen X und Y das Dataframe für die Zufallsvariable Z = X + Y berechnet.

composition <- function(X, Y = X, FUN = "+"){
  stopifnot(is.randomVariable(X), is.randomVariable(Y))
  Z <- data.frame( value = as.vector( outer(X = X$value, Y = Y$value, FUN = FUN) ), 
                   prob = as.vector( outer(X = X$prob, Y = Y$prob, FUN = "*") ) )
  # Beachte: Z$value und Z$prob haben identische Länge, 
  # müssen aber gekürzt werden: Werte können mehrfach vorkommen 
  Z.aggr <- aggregate(x = Z$prob, by = Z[1], FUN = sum)
  names(Z.aggr) <- c("value", "prob")
  return(Z.aggr)
}

# Test:

X2 <- composition(X = X)
X2
#    value       prob
# 1      2 0.02777778
# 2      3 0.05555556
# 3      4 0.08333333
# 4      5 0.11111111
# 5      6 0.13888889
# 6      7 0.16666667
# 7      8 0.13888889
# 8      9 0.11111111
# 9     10 0.08333333
# 10    11 0.05555556
# 11    12 0.02777778

Zeile 1: Oben wurden die Zufallsvariablen X und Y immer addiert, es spricht nichts dagegen sie anders zu verknüpfen (etwa Produkt), daher wird heir eine beliebige Funktion zugelassen, deren default-Wert gleich der Addition gesetzt wird.

Zeile 3 und 4: Mit outer() werden die Kombinationen der Werte (mit der beliebigen Funktion FUN) und die Produkte der Wahrscheinlichkeiten (mit FUN = "*" ) berechnet.

Zeile 7: Schreibt man im Argument by nicht Z$value sondern Z[1] , muss man nicht in eine Liste verpacken.

Zeile 8 und 9: Die Namen der Spalten müssen ausdrücklich gesetzt werden, ansonsten ist das von aggragate() erzeugte Dataframe der Rückgabewert.

Zeile 14: Um die Zufallsvariable für die Augensumme beim zweimaligen Würfeln zu berechnen, muss man nur das Dataframe X setzen und man erhält das bekannte Ergebnis.

Abbildung 4: Die Wahrscheinlichkeiten für die Augensumme beim zweimaligen Würfeln berechnet mit der Faltung der Wahrscheinlichkeiten für das einmalige Würfeln mit Hilfe der Funktion ''composition()''.

Aufgabe:

Berechnen Sie Wahrscheinlichkeiten für die Augensumme bei 4, 8, 16, 32 Würfen eines Laplace-Würfels. Die Ergebnisse sind in Abbildung 5 gezeigt; um die Wahrscheinlichkeiten zu vergleichen, ist die Skalierung der y-Achse zu beachten

Abbildung 5: Entsprechend zu Abbildung 4 werden weitere Faltungen berechnet, die die Wahrscheinlichkeitsverteilung für das 4, 8-, 16, 32-malige Würfeln beschreiben, und als Stabdiagramm dargestellt. Beim Vergleich der Diagramme ist auf die Skalierung der y-Achse zu achten.

Zusammenfassung

Die folgende Tabelle zeigt eine Kurzbeschreibung der in diesem Kapitel vorgestellten Funktionen. Man beachte dabei aber, dass die gezeigten Argumente meist nicht die komplette Liste der Eingabewerte darstellt. Es werden immer nur diejenigen Eingabewerte gezeigt, die hier auch besprochen wurden. Für die allgemeine Verwendung der Funktionen ist es erforderlich die Dokumentation zu Rate zu ziehen.