Die Familie der apply-Funktionen in R Teil 2: Die Verarbeitung mehrerer Listen mit mapply(), Map() und outer()

Die Funktion lapply() ersetzt eine Iteration über die Komponenten einer Liste, wobei auf jede Komponente eine Funktion FUN angewendet wird; die Rückgabewerte werden wieder zu einer Liste zusammengefasst. Entsprechend wird mit mapply() über mehrere Listen iteriert, wobei in jedem Schritt entsprechende Komponenten ausgewählt werden und darauf wird die Funktion FUN angewendet. Die Funktion Map() ist ein Wrapper für mapply(), der die wichtigsten Anwendungsfälle abdeckt. Die meisten Funktionen in R sind vektorisiert, können also nicht nur auf einen Eingabewert, sondern auf einen Vektor angewendet werden. Die Vektorisierung von Funktionen ist ein in R zentrales Konzept, das ein besseres Verständnis der Funktion mapply() liefert. Zuletzt wird die Funktion outer() mit einigen Anwendungen besprochen. Die Funktion outer() besitzt zwei Vektoren (oder Felder) als Eingabewert und baut daraus ein komplexeres Feld auf.

walter

2 Okt. 2019

split-apply-combine mapply()outer()Vectorize()Vektorisierung Map()Liste Polynom Taylor-Entwicklung

Einordnung des Artikels

Einführung in die Informatik
- Einführung in die Programmiersprache R
  - Funktionen in R
    - 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()

Einführung

In Die Familie der apply-Funktionen in R Teil 1: Verarbeitung von Listen mit lapply(), sapply(), vapply() und rapply() wurden Funktionen vorgestellt, die eine Liste verarbeiten. In diesem zweiten Teil über die Familie der apply-Funktionen werden Funktionen untersucht, die mehrere Listen gleichzeitig verarbeiten; dies sind:

Die Funktion mapply(), die zahlreiche Eingabewerte besitzt und somit mehrere Funktionalitäten bereitstellt.
Die Funktion Map(), die einen Wrapper für mapply() bildet; sie kann nicht so fein konfiguriert werden und soll lediglich diejenigen Einsatzmöglichkeiten von mapply() abdecken, die am häufigsten benötigt werden.
Die Funktion outer(), die zwei Vektoren (oder Felder) als Eingabewerte besitzt, die aber völlig anders abgearbeitet werden als durch mapply() oder Map().

Da Map() auf mapply() zurückgeführt werden kann, wird Letzere sehr ausführlich und Erstere nur sehr kurz besprochen.

Zuletzt wird ein umfangreiches Beispiel diskutiert – die Implementierung von Polynom-Funktionen. Mit diesem Beispiel soll ein besseres Verständnis für die mit den apply-Funktionen verbundenen Konzepte erzielt werden:

vektorisierte Funktionen,
Abarbeitung von Listen (entweder einer einzigen Liste oder mehrerer Listen gleichzeitig),
die Vektorisierung von Funktionen.

Die Funktion mapply()

Charakterisierung der Funktion mapply()

So wie die Funktion lapply() dazu verwendet wird, um über die Komponenten einer Liste (oder eines Vektors) zu iterieren, kann man mit mapply() über mehrere Listen (oder Vektoren) gleichzeitig iterieren; das m in mapply() steht für multivariate.

In Abbildung 1 ist die Arbeitsweise von mapply() dargestellt:

Abbildung 1: Die drei Phasen split-apply-combine bei der Anwendung von mapply() auf zwei Vektoren oder Listen mit der Funktion f(); detaillierte Erklärung im Text.

Phase split:
- Es können beliebig viele Listen (oder Vektoren) als Eingabewerte an mapply() im Argument ... übergeben werden; sie sollten identische Längen haben. In Abbildung 1 sind zwei Listen v und w gezeigt.
- Die Listen werden in ihre Komponenten zerlegt.
- Die jeweils entsprechenden Komponenten (mit gleichem Index) werden zu einem Objekt zusammengefügt (blaue Umrandung in Abbildung 1).
Phase apply:
- Diese zusammengefügten Objekte (v[1], w[1]), (v[2], w[2]), ..., (v[n], w[n]) werden an die Funktion f() übergeben.
- Die Funktion f() muss hier zwei Eingabewerte verarbeiten können. Im Allgemeinen sind es beliebig viele Eingabewerte, da an mapply() beliebig viele Listen übergeben werden können.
- Die Funktion mapply() ersetzt wieder eine Iteration über diese neu zusammengesetzten Objekte: es entstehen n Objekte f(v[1], w[1]), f(v[2], w[2]), ..., f(v[n], w[n]) .
Phase combine: Das Verhalten von mapply() in der Phase combine hängt davon ab, ob das Argument SIMPLIFY auf TRUE (default-Wert) oder FALSE gesetzt ist:
- Für SIMPLIFY = FALSE werden die n Rückgabewerte von f() zu einer Liste zusammengefügt.
- Für SIMPLIFY = TRUE werden sie zu einem Vektor, einer Matrix oder einem Feld zusammengefasst – dies hängt davon ab, welchen Rückgabewert f() erzeugt.
- Ist Letzeres nicht möglich (etwa weil die Rückgabewerte unterschiedliche Länge haben), wird auch bei SIMPLIFY = TRUE eine Liste erzeugt.

Die Eingabewerte von mapply() sind:

mapply(FUN, ..., MoreArgs = NULL, SIMPLIFY = TRUE, USE.NAMES = TRUE)

Das Argument FUN steht für die Funktion, die in der apply-Phase angewendet wird.
In ... werden beliebig viele zu verarbeitende Listen (oder Vektoren) eingegeben.
Das Argument MoreArgs ist dazu da, weitere Argumente aufzunehmen, die an FUN weitergereicht werden sollen; da das Argument ... schon für die zu verarbeitenden Listen eingesetzt wird, kann es nicht nochmals verwendet werden (ein Beispiel dazu folgt weiter unten).
Mit SIMPLIFY = TRUE wird das Verhalten in der combine-Phase festgelegt (siehe Abbildung 1 unten).
Das Argument USE.NAMES gibt an, ob das Attribut names an den Rückgabewert weitergegeben wird.

Einfache Anwendungen der Funktion mapply()

Um die in Abbildung 1 dargestellte Arbeitsweise von mapply() besser zu verstehen, werden zunächst einige einfache Beispiele angeführt; sie sind so einfach gewählt, dass man die Anwendung von mapply() durch eine vektorisierte Funktion ersetzen kann.

1. Beispiel:

An die Funktion mapply() werden – wie in Abbildung 1 – zwei Vektoren übergeben; als Funktion FUN wird die Summe gewählt:

v <- c(1, 3, 5, 7)
w <- (1:4)

mapply(FUN = sum, v, w)
# [1]  2  5  8 11

str( mapply(FUN = sum, v, w) )
# num [1:4] 2 5 8 11

# Kontrolle:
v + w
# [1]  2  5  8 11

Zeile 4 und 5: Die beiden Vektoren v und w haben jeweils die Länge 4. Es werden die beiden ersten Komponenten addiert, die beiden zweiten Komponenten addiert und so weiter. Es entsteht ein Vektor mit 4 Komponenten.

Zeile 7 und 8: Da der default-Wert von SIMPLIFY gleich TRUE ist, werden die 4 Komponenten tatsächlich zu einem Vektor (und keiner Liste) zusammengesetzt.

Zeile 11 und 12: Die Summe der Vektoren kann natürlich direkt mit + gebildet werden, die Anwendung von mapply() ist überflüssig.

2. Beispiel: Bilden von Kombinationen zweier Vektoren

Wieder sind zwei Vektoren gleicher Länge gegeben. Ihre Komponenten mit jeweils gleichem Index sollen zu einem Vektor vereinigt werden:

v <- c(1, 3, 5, 7)
w <- (1:4)

mapply(FUN = c, v, w)
#      [,1] [,2] [,3] [,4]
# [1,]    1    3    5    7
# [2,]    1    2    3    4

str(mapply(FUN = c, v, w, SIMPLIFY = FALSE ) )
# List of 4
# $ : num [1:2] 1 1
# $ : num [1:2] 3 2
# $ : num [1:2] 5 3
# $ : num [1:2] 7 4

Zeile 4 bis 7 Es entstehen 4 Vektoren der Länge zwei, die zu einer Matrix zusammengefügt werden – wieder wegen SIMPLIFY = TRUE . Man beachte, dass die 4 Rückgabewerte von FUN spaltenweise aneinandergehängt werden.

Zeile 9 bis 14: Setzt man dagegen SIMPLIFY = FALSE , entsteht eine Liste (mit 4 Komponenten).

Auch in diesem Beispiel ist die Anwendung von mapply() nicht nötig, da die Aufgabe auch mit rbind() (oder cbind()) erledigt werden kann. Aber man erkennt wieder, wie die Eingabe-Vektoren aufgespalten werden und wie FUN auf die neu gebildeten Objekte angewendet wird.

Typische Anwendungen von mapply()

Schon mehrfach wurde die Funktion weightedMean() verwendet, um Funktionen der apply-Familie zu untersuchen. Auch für mapply() lassen sich damit Anwendungs-Beispiele konstruieren, die ohne mapply() nur sehr umständlich zu realisieren wären.

Verwendet wird wieder die einfache Implementierung von weightedMean():

weightedMean <- function(x, m){
  return(sum(x*m) / sum(m))
}

v <- c(1, 3, 5, 7)
w <- (1:4)

weightedMean(x = v, m = w)
# [1] 5

Dazu muss man sich zunächst überlegen: Wie sieht eigentlich eine Fragestellung aus, so dass sie nur mit mapply() und nicht etwa mit lapply() gelöst werden kann?

Die Funktion weightedMean() hat zwei Eingabewerte: den Vektor x, für den der gewichtete Mittelwert berechnet werden soll, und die (positiven) Gewichtungen m.

Hat man mehrere Vektoren x gegeben und der gewichtete Mittelwert soll immer mit einem einzigen Vektor weights berechnet werden, kann man die Aufgabe mit lapply() erledigen. Angenommen die Vektoren x sind in einer Liste lst gespeichert und der Vektor der Gewichtungsfaktoren ist weights, dann benötigt man den Befehl:

lapply(X = lst, FUN = weightedMean, m = weights)

Hier wird über die Liste lst der x-Vektoren iteriert, wobei immer der Vektor der Gewichtungen m = weights verwendet wird; er wird im Argument ... von lapply() an weightedMean() weitergereicht.

Hat man umgekehrt nur einen Vektor x und verschiedene Vektoren für die Gewichtungen, so tauschen x und m lediglich die Rollen.

Erst wenn es mehrere Vektoren x gibt und zu jedem Vektor x ein eigener Vektor mit Gewichtungen m existiert, kann man das Problem nicht mehr mit lapply() lösen. Wenn man dies mit Abbildung 1 vergleicht, dann ist etwa v eine Liste von Vektoren x und w eine Liste von Vektoren m, die in der split-Phase in ihre n Komponenten zerlegt werden. Jeweils die entsprechenden x- und m-Vektoren werden an die Funktion weightedMean() weitergegeben (apply-Phase), so dass insgesamt n gewichtete Mittelwerte berechnet werden. Diese werden in der combine-Phase zu einem Vektor (oder einer Liste) zusammengefügt.

Bevor man mit den eigentlichen Anwendungen beginnt, sollte man sich klarmachen: Welche Berechnung nimmt mapply() vor, wenn als Argument FUN die Funktion weightedMean() und als Argument ... die beiden Vektoren v und w aus dem Beispiel oben übergeben werden?

Das folgende Skript zeigt die Antwort:

v <- c(1, 3, 5, 7)
w <- (1:4)

mapply(FUN = weightedMean, v, w)
# [1] 1 3 5 7

Man erkennt:

Es ensteht ein Vektor der Länge 4, das heißt dass 4 Mittelwerte berechnet wurden. Vergegenwärtigt man sich die Wirkung von mapply() in der Phase split (siehe Abbildung 1 oben), so wird zur Berechnung der ersten Komponente des Ergebnisses jeweils die erste Komponente der beiden Eingabe-Vektoren v und w ausgewählt. Dabei ist v[1] das Argument x in weightedMean() und entsprechend w[1] das Argument m. (Für die anderen Komponenten analog.) Faktisch werden also keine Gewichtungen vorgenommen, der Mittelwert stimmt mit dem Eingabewert x überein, wodurch als Ausgabewert von mapply() der Eingabe-Vektor v reproduziert wird.

Im Folgenden wird eine Anwendung von mapply() gezeigt, die echte gewichtete Mittelwerte berechnet. Dazu werden 4 x-Vektoren und 4 m-Vektoren jeweils in eine Liste verpackt und anschließend mit mapply() abgearbeitet:

x <- list( a = c(1, 2, 3), b = c(4, 5, 6), c = c(7, 8, 9), d = c(10, 11, 12) )
m <- list( ma = c(1, 1, 0), mb = c(1, 0, 1), mc = c(0, 1, 1), md = c(1, 1, 1) )

mapply(FUN = weightedMean, x, m)
#   a    b    c    d 
# 1.5  5.0  8.5 11.0 

str( mapply(FUN = weightedMean, x, m, SIMPLIFY = FALSE) )
# List of 4
# $ a: num 1.5
# $ b: num 5
# $ c: num 8.5
# $ d: num 11

# Realisierung mit Schleife:
means <- vector(mode = "numeric", length = 4)
for(i in seq_along(x)){
  means[i] <- weightedMean(x = x[[i]], m = m[[i]])
}
means
# [1] 1.5 5.0 8.5 11

Jetzt ist x eine Liste aus vier Vektoren (Zeile 1) und ebenso m (Zeile 2).

Werden die Listen x und m an das Argument ... in weightedMean() übergeben, so werden vier gewichtete Mittelwerte berechnet (Zeile 4 bis 6):

der erste gewichtete Mittelwert verwendet den Vektor a aus x mit den Gewichtungen ma,
der zweite gewichtete Mittelwert verwendet den Vektor b und die Gewichtungen mb, ...

Wird SIMPLIFY = FALSE gesetzt (Zeile 8), werden die vier gewichteten Mittelwerte in eine Liste verpackt.

In beiden Fällen wurden die Namen der Liste x an das Ergebnis weitergereicht (der default-Wert für USE.NAMES ist gleich TRUE).

Zum Vergleich wird die Berechnung der gewichteten Mittelwerte in den Zeilen 16 bis 19 mit einer Schleife vorgenommen: Die Berechnung ist deutlich umständlicher, das man sich selbst um die Phase combine kümmern muss (hier mit einem Vektor realisiert); ebenso müsste man auch selber dafür sorgen, dass die Namen im Endergebnis von x übernommen werden.

Die Situation mit Vektoren, die in Listen verpackt sind, ist etwas ungewöhnlich; naheliegender ist es, die Vektoren aus einem Dataframe zu extrahieren.

Um aber mapply() auf Dataframes anzuwenden, muss man erst klären:

Ist es überhaupt möglich im Argument ... von mapply() Dataframes einzugeben?
Wenn ja: wie werden Dataframes in der split-Phase zerlegt: zeilen- oder spaltenweise zerlegt?

Bei der Diskussion von lapply() und sapply() wurde schon gezeigt, dass man mit diesen Funktionen über Dataframes iterieren kann und dass es spaltenweise abgearbeitet wird. Das folgende Beispiel definiert zwei Dataframes (eines für die x-Vektoren, eines für die m-Vektoren) und führt sie mit mapply() der Funktion weightedMean() zu:

a = c(1, 2, 3); b = c(4, 5, 6); c = c(7, 8, 9); d = c(10, 11, 12)
df.x <- data.frame(a = a, b = b, c = c, d = d)

ma = c(1, 1, 0); mb = c(1, 0, 1); mc = c(0, 1, 1); md = c(1, 1, 1)
df.m <- data.frame(ma, mb, mc, md)

df.x
#   a b c  d
# 1 1 4 7 10
# 2 2 5 8 11
# 3 3 6 9 12

df.m
#   ma mb mc md
# 1  1  1  0  1
# 2  1  0  1  1
# 3  0  1  1  1

means <- mapply(FUN = weightedMean, df.x, df.m)
means
#   a    b    c    d 
# 1.5  5.0  8.5 11.0  

means <- mapply(FUN = weightedMean, df.x, df.m, SIMPLIFY = FALSE)
str(means)
# List of 4
# $ a: num 1.5
# $ b: num 5
# $ c: num 8.5
# $ d: num 11

# Kontrolle mit Schleife:
means <- vector(mode = "numeric", length = 4)
for(i in (1:ncol(df.x))){
  means[i] <- weightedMean(x = df.x[ , i], m = df.m[ , i])
}
means
# [1]  1.5  5.0  8.5 11.0

Zeile 1 und 2: Vier Vektoren werden zu einem Dataframe df.x zusammengefasst (x-Vektoren).

Zeile 4 und 5: Vier Vektoren werden zu einem Dataframe df.m zusammengefasst (m-Vektoren).

Zeile 7 bis 17: Ausgabe der beiden Dataframes.

Zeile 19 bis 22: Mit Hilfe von mapply() werden 4 Mittelwerte berechnet, die in einem Vektor abgespeichert werden. Man kann leicht nachvollziehen, dass dabei jeweils die entsprechenden Spalten der beiden Dataframes ausgewählt und an weightedMean() weitergereicht werden.

Zeile 24 bis 30: Setzt man SIMPLIFY = FALSE , werden die 4 Mittelwerte in einer Liste abgelegt.

Zeile 33 bis 38: Realisierung der Berechnung der Mittelwerte mit einer Schleife.

Das Argument MoreArgs von mapply()

Bei den apply-Funktionen, die bisher besprochen wurden, konnten im Argument ... Eingabewerte an die Funktion FUN weitergereicht werden. Da ... nicht zweimal in der Argument-Liste einer Funktion vorkommen kann und es schon für die zu verarbeitenden Vektoren (oder Listen) reserviert ist, benötigt man einen anderen Mechanismus, um Eingabewerte an FUN weiterzugeben. Diese Eingabewerte müssen als Liste verpackt als MoreArgs gesetzt werden.

Das folgende Skript zeigt eine kleine Abwandlung der Funktion weightedMean(); jetzt besitzt sie das zusätzliche Argument debug, mit dem man Zwischenergebnisse ausgeben kann:

weightedMean <- function(x, m, debug = FALSE){
  z <- sum(x * m)
  n <- sum(m)
  if(debug){
    message("weighted sum = ", z)
    message("sum of weights = ", n)
    message("###########")
  }
  return( z / n )
}

Mit dem default-Wert debug = FAlSE erzeugt weightedMean() die selben Ergebnisse wie die oben angegebene Implementierung. Möchte man aus der Funktion mapply() heraus dafür sorgen, dass die Zwischenergebnisse angezeigt werden, muss debug = TRUE im Argument MoreArgs gesetzt werden:

x <- list( a = c(1, 2, 3), b = c(4, 5, 6), c = c(7, 8, 9), d = c(10, 11, 12) )
m <- list( ma = c(1, 1, 0), mb = c(1, 0, 1), mc = c(0, 1, 1), md = c(1, 1, 1) )

mapply(FUN = weightedMean, x, m, MoreArgs = list(debug = TRUE))
# weighted sum = 3
# sum of weights = 2
# ###########
# weighted sum = 10
# sum of weights = 2
# ###########
# weighted sum = 17
# sum of weights = 2
# ###########
# weighted sum = 33
# sum of weights = 3
# ###########
# a    b    c    d 
# 1.5  5.0  8.5 11.0

Zeile 1 und 2: Es werden wieder die Listen mit den x- und m-Werten definiert.

Zeile 4: Neu im Aufruf von mapply() ist, dass debug = TRUE in eine Liste verpackt an das Argument MoreArgs übergeben wird.

Zeile 5 bis 18: Nach den debug-Ausgaben folgt der Rückgabewert des Aufrufes von mapply().

Die Funktion Map() als Wrapper für mapply()

Die Funktion Map() besitzt nur zwei Eingabewerte:

Map(f, ...)

Mit f wird die anzuwendende Funktion bezeichnet und mit ... werden wie bei mapply() die abzuarbeitenden Listen eingegeben.

An der Implementierung von Map() erkennt man, dass die Eingabewerte tatsächlich nur an mapply() weitergereicht werden:

Map
function (f, ...) 
{
    f <- match.fun(f)
    mapply(FUN = f, ..., SIMPLIFY = FALSE)
}
<bytecode: 0x0000000017c874e0>
<environment: namespace:base>

Insbesondere wird mapply() mit SIMPLIFY = FALSE aufgerufen, das heißt es wird als Rückgabewert eine Liste erzeugt. Möchte man sie in einen Vektor verwandeln, muss man dies selber mit unlist() tun.

Die Dokumentation zu Map() findet man unter ?funprog (funprog steht für functional programming).

Anwendung: Vektorisierung von Funktionen am Beispiel der Polynom-Funktionen

Vektorisierte Funktionen

Es wurde in diesem R-Kurs schon mehrfach betont, dass die meisten Funktionen vektorisiert sind, was aus mathematischer Sichtweise eher ungewöhnlich ist, aber viele praktische Vorteile mit sich bringt. So kann man etwa schnell eine Werte-Tabelle für eine Funktion anlegen und diese zeichnen lassen.

Das folgende Skript definiert einen Vektor von 31 x-Werten und berechnet dazu die Exponentialfunktion. Anschließend werden die Punkte graphisch dargestellt (siehe Abbildung 2):

x <- seq(from = -1, to = 1, length.out = 21)
y <- exp(x)

par(bg = "#FFFFF3")
plot(x, y, axes = TRUE, frame.plot = TRUE, xlim = c(-1, 2), ylim =c(-0.5, 7.5), ylab = "exp(x)" )

lines(x = c(-1, 2), y = c(0, 0))
lines(x = c(0, 0), y = c(-0.5, 7.5))

Abbildung 2: Darstellung von Punkten auf dem Graphen der Exponentialfunktion.

Wenn man Funktionen selber implementiert, kann es leicht passieren, dass sie nicht vektorisiert sind. Oft möchte man aber diese Eigenschaft ausnutzen.

In den folgenden Abschnitten wird eine Problemstellung gezeigt – nämlich Polynom-Funktionen mit gegebenen Koeffizienten –, bei der die naheliegende Implementierung eine nicht vektorisierte Funktion liefert. Anschließend wird diskutiert:

Wie kann man eine Funktion vektorisieren?
Welcher Zusammenhang besteht zwischen Vektorisierung und den apply-Funktionen?

Beispiel: Implementierung von Polynom-Funktionen

Abbildung 3 zeigt in Gleichung 1 die Definition einer Polynom-Funktion p(x): Durch die fest vorgegebenen (reellen) Konstanten (a₀, a₁, ..., a_n) wird ein Polynom vom Grad n erzeugt.

Man beachte, dass sich hier ein Problem mit den Bezeichnungen ergibt: Für ein Polynom vom Grad n benötigt man n + 1 Konstanten a_i, i = 0, ..., n, wobei man üblicherweise die erste Konstante mit 0 indiziert, da sie zur Potenz x⁰ = 1 gehört. Wird dies später implementiert, werden die Konstanten a_i zu einem Vektor zusammengefasst, bei dem der Index bei 1 beginnt. Oder umgekehrt: Ein Vektor a der Länge n erzeugt ein Polynom vom Grad n - 1.

Abbildung 3: Formeln zur Berechnung von Polynom-Funktionen des Grades n: für eine Variable in Gleichung 1, für 2 Variablen in Gleichung 4. Als Beispiele sind die Formeln für die Taylor-Entwicklung der Exponentialfunktion um x = 0 beziehungsweise um (x, y) = (0, 0) gezeigt.

In Abbildung 3 ist in Gleichung 2 und 3 ein spezielles Beispiel zu sehen: Mit der Wahl der Konstanten aus Gleichung 2 erhält man eine Approximation der Exponentialfunktion; dies ist gerade die Taylor-Entwicklung um x = 0.

Das folgende Skript zeigt eine naheliegende Implementierung der Funktion polynomial(), die zu einem gegebenen x-Wert und Konstanten a ein Polynom erzeugt und es an der Stelle x auswertet:

polynomial <- function(x, a = 1){
  pows <- cumprod( c( 1, rep(x, times = length(a) - 1  ) ) )
  return( sum(a * pows) )
}

polynomial(x = 0.5)
# 1
polynomial(x = 0.5, a = c(1, 1, 1))
# [1] 1.75

Man beachte dabei die Komplikation mit der Länge von a: Hat a die Länge n, wird ein Polynom vom Grad n - 1 erzeugt.

In Zeile 2 wird zuerst der Vektor (1, x, x, ..., x) gebildet, der die selbe Länge hat wie a. Mit dem kumulierten Produkt cumprod() werden anschließend die Potenzen von x erzeugt:

x⁰, x¹, ..., x^n-1.

Zeile 3: Das Polynom wird an der Stelle x ausgewertet, indem man das Skalarprodukt aus a und den Potenzen von x (also pows) bildet.

Zeile 6 bis 9: Zwei simple Anwendungs-Beispiele, einmal für das konstante Polynom p(x) = 1 und einmal für

p(x) = 1 + x + x²,

beide ausgewertet für x = 0.5.

Die Vektorisierung einer Funktion

Die Implementierung der Polynom-Funktion wirkt sehr einfach, hat aber einen Nachteil: so wie polynomial() den Rückgabewert berechnet, kann die Funktion nicht vektorisiert eingesetzt werden.

Das folgende Beispiel demonstriert dies:

polynomial(x = c(0.5, 1), a = c(1, 1, 1))
# [1] 2.5
# Warning message:
#   In a * pows :
#   longer object length is not a multiple of shorter object length

Wenn polynomial() vektorisiert wäre, müsste sich als Ergebnis der Vektor (1.75, 3) ergeben. An der Ausgabe erkennt man aber, dass ein Rückgabewert berechnet wird und eine Warnung ausgegeben wird.

Dazu sollte man sich klarmachen: Wo wird in der Implementierung eine Operation eingesetzt, die die Vektorisierung verhindert?

In der Implementierung von polynomial() wird in Zeile 2 die Funktion rep() verwendet: die Implementierung ist zwar syntaktisch korrekt, berechnet aber für x-Vektoren nicht das gewünschte Polynom. Das folgende Skript zeigt dies:

rep(x = c(0.5, 1), times = 3  )
# [1] 0.5 1.0 0.5 1.0 0.5 1.0

Wird mit diesem Vektor weitergerechnet, erzeugt cumprod() abwechselnd Potenzen von x[1] und x[2] . Um eine vektorisierte Version von polynomial() zu erhalten, müsste man also die Komponenten von x einzeln verwenden, um die Potenzen pows und anschließend das Skalarprodukt zu berechnen. Erst dann werden die Ergebnisse zu einem Vektor zusammengefasst.

Es sollen jetzt mehrere Möglichkeiten vorgestellt werden, wie man die Funktion polynomial() vektorisieren kann; die wird zugleich zu einem besseren Verständnis der apply-Funktionen führen.

1. Möglichkeit: polynomial.iterate()

Über die Komponenten von x wird iteriert; für jede einzelne Komponente werden die oben gezeigten Berechnungen ausgeführt. Dies ist mit einer Schleife verbunden, die meist aufwendiger ist als der Einsatz einer geeigneten apply-Funktion.

2. Möglichkeit: polynomial.sapply()

Die Implementierung verwendet die Funktion sapply(), um über die x-Werte zu iterieren.

3. Möglichkeit: polynomial.vec()

Man verwendet die ursprüngliche Implementierung von polynomial() und verwandelt sie mit Hilfe von Vectorize() in eine vektorisierte Funktion.

1. Möglichkeit: polynomial.iterate()

polynomial.iterate <- function(x, a = 1){
  polynomial <- function(x, a = 1){
    stopifnot(length(x) == 1)
    pows <- cumprod( c( 1, rep(x, times = length(a) - 1  ) ) )
    return( sum(a * pows) )
  }
  y <- vector(mode = "numeric", length = length(x))
  for(i in seq_along(along.with = x)){
    y[i] <- polynomial(x = x[i], a = a)
  }
  return(y)
}

polynomial.iterate(x = c(0.5, 1), a = c(1, 1, 1))
# [1] 1.75 3.00

Zeile 2: Die Funktion polynomial() wird innerhalb von polynomial.iterate() verwendet. Dazu wird der Befehl sicherheitshalber der Befehl stopifnot(length(x) == 1) eingefügt, damit für x kein echter Vektor übergeben werden kann (Zeile 3).

Zeile 7: Es wird der Rückgabewert als Vektor mit der Länge von x vorbereitet.

Zeile 8 bis 10: Die for-Schleife iteriert über die x-Werte und berechnet jeweils den Wert des Polynoms an der Stelle x.

Zeile 11: Der Vektor der y-Werte wird zurückgegeben.

Zeile 14 und 15: Der Test mit den x-Werten (0.5, 1) und dem Polynom

p(x) = 1 + x + x²

liefert jetzt die richtigen Ergebnisse.

2. Möglichkeit: polynomial.sapply()

Die for-Schleife aus der letzten Implementierung lässt sich natürlich vermeiden. Man muss sich nur fragen: welche apply-Funktion ersetzt hier die Schleife?

Iteriert wird über die Werte von x, das heißt für jeden x-Wert soll ein Rückgabewert berechnet werden; die Rückgabewerte werden dann am Besten zu einem Vektor zusammengefasst – eine Liste ist hier nicht nötig, da jeder Rückgabewert eine Zahl ist. Der Vektor a ist dagegen bei jedem Iterationsschritt identisch.

Es ist also nicht nötig, mapply() einzusetzen; und da der Rückgabewert insgesamt ein Vektor sein soll, wird man sapply() wählen. Eine mögliche Realisierung könnte wie folgt aussehen – hier wird wieder die Funktion polynomial() eingesetzt:

polynomial.sapply <- function(x, a = 1){
  polynomial <- function(x, a = 1){
    stopifnot(length(x) == 1)
    pows <- cumprod( c( 1, rep(x, times = length(a) - 1  ) ) )
    return( sum(a * pows) )
  }
  return(sapply(X = x, FUN = polynomial, a = a))
}

polynomial.sapply(x = c(0.5, 1), a = c(1, 1, 1))
# [1] 1.75 3.00

Zeile 7: Hier werden zwei Anweisungen verschachtelt:

Die Funktion sapply() wird mit dem Vektor X und der Funktion polynomial() aufgerufen; an polynomial() muss der Vektor a mit den Koeffizienten des Polynoms weitergereicht werden.
Die Funktion sapply() erzeugt einen Vektor (für jede Komponente von x entsteht ein Zahlenwert), der sofort zurückgegeben werden kann.

Zeile 10: Der Test mit x = c(0.5, 1) liefert wieder die richtigen Ergebnisse.

3. Möglichkeit: polynomial.vec()

Im Paket base gibt es die Funktion Vectorize(), die bisher noch nicht besprochen wurde:

Vectorize(FUN, vectorize.args = arg.names, SIMPLIFY = TRUE, USE.NAMES = TRUE)

Sie bietet die Möglichkeit, eine Funktion FUN in eine vektorisierte Funktion zu verwandeln; man muss dazu angeben, welche Argumente vektorisiert werden sollen – hier ist es lediglich der Eingabewert x. Die weiteren Argumente erklären sich selbst. Rückgabewert ist dann eine Funktion, die man sich unter einem geeigneten Namen abspeichern und die man wie jede andere Funktion einsetzen kann.

Das folgende Skript zeigt wie man polynomial() vektorisieren kann:

polynomial <- function(x, a = 1){
  stopifnot(length(x) == 1)
  pows <- cumprod( c( 1, rep(x, times = length(a) - 1  ) ) )
  return( sum(a * pows) )
}

polynomial.vec <- Vectorize(FUN = polynomial, vectorize.args = "x")

polynomial.vec(x = 0.5, a = c(1, 1, 1))
# polynomial(x = 0.5, a = c(1, 1, 1))

polynomial.vec(x = c(0.5, 1), a = c(1, 1, 1))
# [1] 1.75 3.00

Zeile 1 bis 5: Ausgangspunkt ist wieder die Funktion polynomial(), die mit dem stopifnot()-Befehl definiert wird.

Zeile 7: Mit Hilfe von Vectorize() wird daraus eine neue Funktion definiert, nämlich polynomial.vec().

Zeile 9 bis 13: Die Tests zeigen, dass polynomial.vec() sowohl mit einem x-Wert als auch mit einem Vektor x wie gewünscht arbeitet.

Es stellt sich natürlich die Frage, wie Vectorize() arbeitet – es kann ja nicht sein, dass der Quelltext der Funktion FUN analysiert und dort die passende apply-Funktion eingebaut wird. Ein Blick auf den Quelltext von Vectorize() gibt die Antwort:

function (FUN, vectorize.args = arg.names, SIMPLIFY = TRUE, USE.NAMES = TRUE) 
{
    arg.names <- as.list(formals(FUN))
    arg.names[["..."]] <- NULL
    arg.names <- names(arg.names)
    vectorize.args <- as.character(vectorize.args)
    if (!length(vectorize.args)) 
        return(FUN)
    if (!all(vectorize.args %in% arg.names)) 
        stop("must specify names of formal arguments for 'vectorize'")
    collisions <- arg.names %in% c("FUN", "SIMPLIFY", "USE.NAMES", 
        "vectorize.args")
    if (any(collisions)) 
        stop(sQuote("FUN"), " may not have argument(s) named ", 
            paste(sQuote(arg.names[collisions]), collapse = ", "))
    FUNV <- function() {
        args <- lapply(as.list(match.call())[-1L], eval, parent.frame())
        names <- if (is.null(names(args))) 
            character(length(args))
        else names(args)
        dovec <- names %in% vectorize.args
        do.call("mapply", c(FUN = FUN, args[dovec], MoreArgs = list(args[!dovec]), 
            SIMPLIFY = SIMPLIFY, USE.NAMES = USE.NAMES))
    }
    formals(FUNV) <- formals(FUN)
    FUNV
}

Die Eingabewerte von FUN, die zu vektorisieren sind, werden an mapply() weitergegeben (siehe Zeile 22). Damit geschieht eigentlich nichts anderes als in der Funktion polynomial.sapply(), die einen geeigneten Wrapper um polynomial() gelegt hat. Da bei polynomial() nur ein Eingabewert vektorisiert werden musste, hat die Funktion sapply() ausgereicht, im Allgemeinen muss man mapply() einsetzen.

Anwendung der vektorisierten Funktion: Darstellung von Polynom-Funktionen

Wenn schon der Aufwand betrieben wurde, um die Funktion polynomial() zu vektorisieren, soll sie auch noch in ihrem Einsatz gezeigt werden. In Abbildung 2 war die Exponentialfunktion für einige ausgewählte x-Werte zu sehen. Gemäß Abbildung 3 kann man leicht ein approximierendes Polynom berechnen und dieses mit der Exponentialfunktion vergleichen.

Das nächste Skript führt folgende Berechnungen aus:

Es werden die Koeffizienten berechnet, um die Exponentialfunktion mit einem Polynom 0., 2., ..., 10. Grades zu approximieren.
Die Approximation wird verwendet, um aus den x-Werten wie in Abbildung 2 die Funktionswerte zu berechnen.
Die Exponentialfunktion und die approximierende Funktion werden graphisch dargestellt (genauer: die ausgewählten Punkte auf den entsprechenden Graphen).

Verwendet wird dabei jeweils die Funktion polynomial.sapply().

# Koeffizienten für die Taylor-Entwicklung der Exponentialfunktion
v <- c(1, (1:10))
coeff <- 1  / cumprod(v)

# x-Werte und zugehörige Punkte der Exponentialfunktion
x <- seq(from = -1, to = 2, length.out = 31)
y.exp <- exp(x)

# Liste der y-Werte der Approximation
y <- vector(mode = "list", length = 6)
for(i in (1:6)){
  y[[i]] <- polynomial.sapply(x, a = head(coeff, n = 2*i - 1))
}

# Plots im Gitter vorbereiten
rs <- 2; cls <- 3
par(mfrow = c(rs, cls))
par(bg = "#FFFFF3")

# Plots: in jedem der 6 Plots wird die Exponentialfunktion (rot) und die Approximation (blau) gezeigt
for( i in (1:(rs*cls)) ){
  plot(x = x, y = y.exp, col = "red", type = "l", ylab = "exp(x)")
  points(x = x, y = y[[i]], col = "blue", type = "h")
}

Zur besseren Unterscheidung wird die Approximation mit Histogrammen dargestellt (blau) und die Punkte auf der Exponentialfunktion werden mit Linien verbunden (rot).

Abbildung 4: Darstellung der Exponentialfunktion und ihrer Approximation durch die Taylor-Entwicklung um x = 0 (mit Polynomen vom Grad 0, 2, ..., 10).

Aufgaben:

Überprüfen Sie, ob
- in Zeile 2 und 3 tatsächlich die Koeffizienten aus Gleichung 2 in Abbildung 3 berechnet werden und
- in Zeile 13 tatsächlich Polynome vom Grad 0, 2, ..., 10 berechnet werden.
Berechnen Sie die maximale Differenz zwischen den Werten der Exponentialfunktion und der jeweiligen Approximation. Bei welchem x-Wert ist der Betrag der Differenz maximal?
In den Abbildungen 5 und 6 wird die Gaußsche Glockenkurve f(x) = exp(-x²) sowie ihre Approximation durch Taylor-Polynome vom Grad 0, 2, ..., 16. Entwickelt wird jeweils um x = 0. Geben Sie die Taylor-Polynome an und implementieren Sie die Ausgabe der Graphen wie in Abbildung 6.

Abbildung 5: Darstellung der Gaußschen Glockenkurve (Erklärung oben).

Abbildung 6: Approximation der Gaußschen Glockenkurve durch die Taylor-Entwicklung um x = 0 (mit Polynomen vom Grad 0, 2, ..., 16).

Die Vektorisierung einer Funktion in mehreren Argumenten

Es soll jetzt der Fall untersucht werden, in dem eine Funktion in zwei Argumenten vektorisiert wird. Anbieten würde sich dazu die Funktion aus Gleichung 4 in Abbildung 3, also ein Polynom in den zwei Variablen x und y. Da allgemeine Polynome mit ihrer Vielzahl von Koeffizienten und den Kombinationen der Potenzen sehr schwer zu handhaben sind, werden hier nur Polynome zweiten Grades betrachtet. Das sind Funktionen in x und y, die nur die Terme aus den ersten zwei Zeilen aus Gleichung 4 in Abbildung 3 enthalten. Eine derartige Funktion wird durch sechs Koeffizienten definiert.

Die Implementierung erfolgt in Anlehnung an polynomial(), wobei das Skalarprodukt in zwei Anteile zerlegt wird:

einmal für den konstanten und die beiden linearen Anteile und
einmal für die quadratischen Anteile.

quad.2 <- function(x, y, a){
  # stopifnot(length(x) == 1, length(y) == 1, length(a) == 6)
  lin <- c(1, x, y)
  quad <- c(x*x, x*y, y*y)
  return( sum( head(a, n = 3) * lin ) + sum( tail(a, n = 3) * quad ) )
}

A <- c(1, 1, 1, 0.5, 1, 0.5)
quad.2(x = 0.5, y = 1, a = A)
# [1] 3.625
quad.2(x = 1, y = 1, a = A)
# [1] 5

In Zeile 2 wird dafür gesorgt, dass die Berechnung abbricht, wenn die Argumente mit ungeeigneten Längen gesetzt werden. Man kann sich leicht davon überzeugen, dass die Funktion quad.2() nicht vektorisiert ist – auch wenn man den stopifnot()-Befehl weglässt.

Möchte man quad.2() vektorisieren, kann man dies entweder mit einer geeigneten apply-Funktion selber implementieren oder mit der Funktion Vectorize(). Im Folgenden werden beide Varianten gezeigt.

1. Vektorisierung mit mapply():

Da jetzt sowohl x als auch y mehrere Werte annehmen können, muss man mit der Funktion mapply() arbeiten. Die Implementierung erfolgt nach dem Muster wie bei polynomial.sapply(): Man legt einen Wrapper um die Funktion quad.2(); die Berechnung eines Eingabewertes (x, y) erledigt also quad.2() und der Wrapper sorgt dafür, dass alle Werte (x, y) abgearbeitet werden.

quad.2.mapply <- function(x, y , a) {
  stopifnot(length(x) == length(y))
  
  quad.2 <- function(x, y, a){
    stopifnot(length(x) == 1, length(y) == 1, length(a) == 6)
    lin <- c(1, x, y)
    quad <- c(x*x, x*y, y*y)
    return( sum( head(a, n = 3) * lin ) + sum( tail(a, n = 3) * quad ) )
  }
  
  return( mapply(FUN = quad.2, x, y, MoreArgs = list(a = a)) )
}

quad.2.mapply(x = 0.5, y = 1, a = A)
# [1] 3.625
quad.2.mapply(x = 1, y = 1, a = A)
# [1] 5
quad.2.mapply(x = c(0.5, 1), y = c(1, 1), a = A)
# [1] 3.625 5.000

Die Funktion mapply() (siehe Zeile 11) sorgt dafür, dass der Reihe nach die entsprechenden Komponenten von x und y ausgewählt werden und an quad.2() weitergegeben werden. Mit Hilfe des Argumentes MoreArgs = list(a = a) wird die Liste a, die beim Aufruf von quad.2.mapply() übergeben wurde (etwa in Zeile 14), an quad.2() weitergereicht.

2. Vektorisierung mit Vectorize():

Alternativ wird wieder Vectorize() eingesetzt. Man geht von der – nicht vektorisierten – Funktion quad.2() aus und definiert eine neue Funktion, indem man angibt, dass die beiden Argumente x und y vektorisiert werden sollen:

quad.2 <- function(x, y, a){
  stopifnot(length(x) == 1, length(y) == 1, length(a) == 6)
  lin <- c(1, x, y)
  quad <- c(x*x, x*y, y*y)
  return( sum( head(a, n = 3) * lin ) + sum( tail(a, n = 3) * quad ) )
}

quad.2.vec <- Vectorize(FUN = quad.2, vectorize.args = c("x", "y"))

quad.2.vec(x = c(0.5, 1), y = c(1, 1), a = A)
#[1] 3.625 5.000

Der Test mit den Koeffizienten A von oben liefert identische Ergebnisse wie mit quad.2.mapply().

Anwendung: Darstellung von Polynomen in zwei Variablen

Die Funktion quad.2.mapply() ist somit vektorisiert und mit ihr kann man für beliebige Mengen von (x, y)-Werten schnell die Funktionswerte berechnen und graphisch darstellen; sämtliche Polynome in zwei Variablen bis zum Grad 2 können mit ihr ausgewertet werden. Nach Gleichung 5 in Abbildung 3 kann man zum Beispiel die Koeffizienten festlegen, um die Funktion exp(x) · exp(y) = exp(x + y) durch Polynome zu approximieren.

Das folgende Skript erzeugt drei Graphiken:

Die lineare Näherung für exp(x) · exp(y) in Abbildung 7,
die quadratische Näherung für exp(x) · exp(y) in Abbildung 8,
die Funktion exp(x) · exp(y) in Abbildung 9.

Dabei sind die Näherungen jeweils die Taylor-Entwicklung um (x, y) = (0, 0).

# 1. Lineare Näherung für exp(x) * exp(y)

# Koeffizienten für die lineare Näherung:
L = c(1, 1, 1, 0, 0, 0)

# x- und y-Werte:
x <- seq(from = -2, to = 2, length.out = 21)

# Berechnung der z-Werte:
z <- outer(X = x, Y = x, FUN = quad.2.mapply, a = L)

# Graphik:
persp(x = x, y = x, z, col ="blue", xlab = "x", ylab = "y",
      nticks = 5, axes = TRUE, box = TRUE, ticktype = "detailed")

# 2. Quadratische Näherung für exp(x) * exp(y)

# Koeffizienten für die quadratische Näherung:
Q <- c(1, 1, 1, 0.5, 1, 0.5)

# Berechnung der z-Werte:
z <- outer(X = x, Y = x, FUN = quad.2.mapply, a = Q)

# Graphik:
persp(x = x, y = x, z, col ="blue", xlab = "x", ylab = "y",
      nticks = 5, axes = TRUE, box = TRUE, ticktype = "detailed")

# 3. exp(x) * exp(y)

# Berechnung der z-Werte:
z <- outer(X = x, Y = x, FUN = function(x, y){exp(x)*exp(y)})

# Graphik:
persp(x = x, y = x, z, col ="blue", xlab = "x", ylab = "y",
      nticks = 5, axes = TRUE, box = TRUE, ticktype = "detailed")

Zeile 7: Zum Erzeugen der Graphik wird ein Rechteck in der x-y-Ebene festgelegt, hier das Quadrat mit Seitenlänge 4 um den Ursprung sowie die Dichte der Punkte (hier liegen 21 Punkte auf einer Strecke der Länge 4).

Zeile 10, 22 und 31: Mit Hilfe von outer() werden die z-Werte zu den 21 · 21 Punkten berechnet.

Zeile 13, 25 und 34: Die Funktion persp() zum Plotten wird den x-y- und z-Werten sowie weiteren Argumenten (für die Darstellung) konfiguriert.

Abbildung 7: Die Taylor-Entwicklung für exp(x + y) um (x, y) = (0, 0) mit einer linearen Funktion in x und y.

Abbildung 8: Die Taylor-Entwicklung für exp(x + y) um (x, y) = (0, 0) mit einer quadratischen Funktion in x und y.

Abbildung 9: Die Funktion exp(x + y). Zum Vergleich mit den Taylor-Entwicklungen ist darauf zu achten, das die z-Achsen unterschiedlich skaliert sind.

Die Funktion outer()

Eigenschaften der Funktion outer()

Die Funktion outer() wurde in Matrizen in R: der Datentyp matrix besprochen, dort allerdings unter dem Aspekt, dass man mit ihr Matrizen erzeugen kann. In der typischen Anwendung von outer() werden zwei Vektoren eingegeben, von denen das Kreuzprodukt (oder äußeres Produkt) gebildet wird, also alle möglichen Kombinationen der Komponenten der Vektoren. Auf jede dieser Kombinationen wird eine Funktion FUN angewendet.

Das Paradebeispiel für die Anwendung von outer() ist die Augensumme beim zweimaligen Würfeln zu berechnen:

v <- (1:6)
m <- outer(X = v, Y = v, FUN = "+")
m
#      [,1] [,2] [,3] [,4] [,5] [,6]
# [1,]    2    3    4    5    6    7
# [2,]    3    4    5    6    7    8
# [3,]    4    5    6    7    8    9
# [4,]    5    6    7    8    9   10
# [5,]    6    7    8    9   10   11
# [6,]    7    8    9   10   11   12

Die Funktion outer() soll jetzt mit den apply-Funktionen verglichen werden: Einerseits scheint sie in die Familie zu gehören, da mit ihr zwei Vektoren (oder allgemein Felder) mit einer Funktion FUN verknüpft werden. Andererseits passt sie nicht so recht in die Familie, da mit apply-Funktionen meist Iterationen ausgeführt werden, die die Komplexität der eingesetzten Daten reduzieren (man denke etwa an die Berechnung der gewichteten Mittelwerte oben); dagegen werden mit outer() aus einfachen Eingabewerten komplexere Strukturen aufgebaut.

Die Eingabewerte von outer() sind die beiden Vektoren (oder Felder) X und Y, deren Komponenten mit der Funktion FUN verknüpft werden. Weitere Argumente können mit ... an FUN weitergereicht werden. Der default-Wert für FUN ist die Multiplikation; wird sie verwendet kann die Funktion outer() durch den Operator %o% ersetzt werden:

outer(X, Y, FUN = "*", ...)
X %o% Y

In Abbildung 10 wird versucht, die Arbeitsweise von outer() darzustellen. Es werden zwei Vektoren v und w gezeigt, aus denen eine Matrix aufgebaut wird; im Allgemeinen sind die Eingabewerte Felder beliebiger Dimension, aus denen ein höher-dimensionales Feld aufgebaut wird. Der Dimensionsvektor des Rückgabewertes von outer() setzt sich aus den Dimensionsvektoren der Eingabewerte zusammen:

dim(outer(X, Y, FUN)) = (dim(X), dim(Y)).

Man erkennt an Abbildung 10 wieder die Phasen split-apply-combine:

In der Phase split werden die Eingabe-Vektoren in ihre Komponenten zerlegt.
In der Phase apply werden alle möglichen Kombinationen der Komponenten der Eingabewerte gebildet und darauf wird die Funktion FUN angewendet.
In der Phase combine werden diese Rückgabewerte zu einer Matrix zusammengesetzt (allgemein zu einem Feld).

Abbildung 10: Die Phasen split-apply-combine dargestellt für die Funktion outer() mit zwei Eingabe-Vektoren und einer Matrix als Rückgabewert.

Eine wichtige Anforderung an FUN in outer() ist, dass sie vektorisiert sein muss; im vorigen Abschnitt wurde dazu erklärt, wie man eine Funktion vektorisieren kann. Auch in den folgenden Anwendungen wird darauf eingegangen.

Anwendungen

Additions- und Multiplikationstabellen

In Rechnen mit Restklassen: Teilbarkeitsregeln wurde ausführlich besprochen:

wie man Additions- und Multiplikationstabellen erzeugt,
wie man sie in anderen Zahlensytemen erzeugt,
welche mathematischen Hintergründe (wie Rechnen mit Restklassen, Teilbarkeitsregeln) relevant sind.

Für Additionstabellen muss man FUN = "+" setzen, für Multiplikationstabellen kann man den default-Wert FUN = "*" verwenden. Für andere Zahlensyteme muss man entweder die Funktion FUN anpassen oder die entstehende Matrix nachbearbeiten.

Die entsprechenden Beispiele sollen hier nicht wiederholt werden.

Wertetabellen für zwei-dimensionale Funktionen

Eine Verallgemeinerung der Additions- und Multiplikationstabellen besteht darin, dass man beliebige Funktionen f() einsetzt, um das äußere Produkt outer() zu bilden. Es entstehen dann Wertetabellen für f(). Diese wurden etwa oben verwendet, um die Graphiken zu erzeugen.

An einfachen Beispielen soll die Vorgehensweise nochmals erläutert werden.

1. Beispiel: Wertetabelle für eine zwei-dimensionale Polynom-Funktion

Oben wurde die Funktion quad.2(x, y, a) implementiert, die zu einem gegebenen Wert im Definitionsbereich (x, y) und einem Vektor a für die Koeffizienten den Wert des Polynoms 2. Grades berechnet; dabei wird eine Konvention verwendet, in welcher Reihenfolge die Werte von a die Koeffizienten des Polynoms bestimmen. Man kann diese Konvention leicht an den ersten zwei Zeilen von Gleichung 4 in Abbildung 3 ablesen: der erste Koeffizient steht für den konstanten Term, die nächsten zwei Koeffizienten für die linearen Terme und die letzten drei Koeffizienten für die quadratischen Terme.

Die Implementierung von quad.2() wird nochmals gezeigt:

quad.2 <- function(x, y, a){
  # stopifnot(length(x) == 1, length(y) == 1, length(a) == 6)
  lin <- c(1, x, y)
  quad <- c(x*x, x*y, y*y)
  return( sum( head(a, n = 3) * lin ) + sum( tail(a, n = 3) * quad ) )
}

Versucht man jetzt eine Wertetabelle mit Hilfe von quad.2() für ein spezielles Polynom zu erzeugen, erhält man eine Fehlermeldung:

A <- c(1, 1, 1, 0.5, 1, 0.5)
x <- (-2:2)
z <- outer(X = x, Y = x, FUN = quad.2, a = A)
# Error in dim(robj) <- c(dX, dY) : 
#   dims [product 25] do not match the length of object [1]

Zeile 1: Die Koeffizienten werden so gewählt wie für Abbildung 8.

Zeile 2: Damit die Ausgaben übersichtlicher bleiben, wird nur ein Gitter mit 5 · 5 Punkten gewählt.

Zeile 3: Die Funktion outer() wird aufgerufen. Es wird das Kreuzprodukt der 5 Werte aus x mit sich selbst gebildet, wobei 25 Punkte in der x-y-entstehen. Der Funktionswert wird mit der quadratischen Funktion berechnet, die durch die Koeffizienten A konfiguriert wird.

Zeile 4: Obwohl alles so plausibel erscheint, erhält man eine Fehlermeldung, die nicht leicht zu verstehen ist.

Wie oben schon diskutiert wurde, ist die Funktion quad.2() nicht vektorisiert und die Berechnung der Komponenten der Matrix für den Rückgabewert erfolgt in outer() vektorisiert mit allen 25 Eingabewerten.

Das heißt man muss nur für quad.2() die vektorisierte Version einsetzen und kann damit das Skript ausführen:

# vektorisierte Version von quad.2():
quad.2.mapply <- function(x, y , a) {
  stopifnot(length(x) == length(y))
  
  quad.2 <- function(x, y, a){
    stopifnot(length(x) == 1, length(y) == 1, length(a) == 6)
    lin <- c(1, x, y)
    quad <- c(x*x, x*y, y*y)
    return( sum( head(a, n = 3) * lin ) + sum( tail(a, n = 3) * quad ) )
  }
  
  return( mapply(FUN = quad.2, x, y, MoreArgs = list(a = a)) )
}

A <- c(1, 1, 1, 0.5, 1, 0.5)
x <- (-2:2)
z <- outer(X = x, Y = x, FUN = quad.2.mapply, a = A)
print(z, digits = 3)
#      [,1] [,2] [,3] [,4] [,5]
# [1,]  5.0  2.5  1.0  0.5  1.0
# [2,]  2.5  1.0  0.5  1.0  2.5
# [3,]  1.0  0.5  1.0  2.5  5.0
# [4,]  0.5  1.0  2.5  5.0  8.5
# [5,]  1.0  2.5  5.0  8.5 13.0

Wenn man jetzt die Zuordnung zwischen den Zeilen- und Spaltennummern zu den (x, y)-Werten herstellt, kann man die Zahlenwerte leicht nachvollziehen. Der Funktionswert für den Ursprung (x, y) = (0, 0) ist zum Beispiel im Zentrum der Matrix zu sehen.

Multiplikation von Polynomen

Die Multiplikation von Polynomen ist meist sehr mühsam und fehleranfällig. Abbildung 11 zeigt eine Möglichkeit, wie man sich die Berechnung erleichtern kann: Die Terme werden nach Potenzen sortiert und in eine Multiplikationstabelle eingetragen. Die Terme in den Nebendiagonalen besitzen stets die selbe Potenz; sie müssen nur noch addiert werden.

In der Mitte sind nur noch die Koeffizienten der Polynome eingetragen; wenn man sie wieder richtig den Potenzen zuordnet, kann man das Ergebnis-Polynom rekonstruieren. Rechts sind die Nebendiagonalen der Multiplikationstabellen farbig markiert.

Abbildung 11: Die Multiplikation von Polynomen erleichtert durch das Kreuzprodukt der Summanden (links) beziehungsweise Koeffizienten (mitte). Die Terme auf den Nebendiagonalen gehören jeweils zum selben Grad (rechts).

Von Abbildung 11 ausgehend ist es nicht mehr schwer eine Funktion zu implementieren, die zwei Polynome multipliziert, wobei die Polynome durch ihre Koeffizienten gegeben sind. Rückgabewert ist der Vektor der Koeffizienten des Ergebnis-Polynoms. Dabei wird die Konvention für die Anordnung der Koeffizienten wie in Gleichung 1 in Abbildung 3 verwendet. Das heißt aber, dass der konstante Term zuerst angegeben wird und die höchste Potenz zuletzt – also genau anders herum als in der üblichen Sortierung in der Mathematik (siehe auch Abbildung 11: am Ende werden die Koeffizienten des Ergebnis-Polynoms in der Reihenfolge angegeben, die der Konvention entspricht).

Die Berechnung des Produktes zweier Polynome erfolgt in den Schritten:

Die Multiplikationstabelle der Koeffizienten wird erstellt.
Die Komponenten auf den Nebendiagonalen werden addiert.
Diese Summen liefern – in der richtigen Reihenfolge – die Koeffizienten des Ergebnis-Polynoms.

Am Beispiel der Polynome aus Abbildung 11 werden diese Schritte durchgespielt, anschließend wird eine Funktion definiert, die aus zwei Polynomen das Produkt-Polynom berechnet.

1. Schritt:

Die Multiplikationstabelle wird mit outer() erzeugt, da multipliziert wird, kann man mit dem Operator %o% arbeiten:

# Koeffizienten der Polynome aus Abbildung 11:
f <- c(1, 0, -3, 1)
g <- c(3, -2, 1)

# Multiplikationstabelle:
m <- f %o% g
m
#     [,1] [,2] [,3]
# [1,]    3   -2    1
# [2,]    0    0    0
# [3,]   -9    6   -3
# [4,]    3   -2    1

Man beachte, dass im Vergleich zu Abbildung 11 die Reihenfolgen vertauscht sind.

2. Schritt:

Die eigentlich schwierige Aufgabe besteht darin, die Nebendiagonalen der Matrix m herauszufiltern. Dabei hilft aber folgende Beobachtung: Addiert man für jede Komponente der Matrix m den Zeilen- und Spalten-Index, so stehen auf den Nebendiagonalen immer identische Werte. Die Matrix, die für jede Komponente die Summe aus Zeilen- und Spalten-Index anzeigt, kann man wieder mit outer() erzeugen:

secondaries <- outer(X = seq_along(along.with = f), Y = seq_along(along.with = g), FUN = "+")
secondaries
#    [,1] [,2] [,3]
# [1,]    2    3    4
# [2,]    3    4    5
# [3,]    4    5    6
# [4,]    5    6    7

Man erkennt, dass auf den Nebendiagonalen identische Wert stehen. Wenn man jetzt noch von jedem Wert 1 subtrahiert, ist der Zahlenwert in der Matrix genau der Index der Komponente, für den die Nebendiagonale steht; besser ist also:

secondaries <- outer(X = seq_along(along.with = f) - 1, Y = seq_along(along.with = g), FUN = "+")
secondaries
#     [,1] [,2] [,3]
# [1,]    1    2    3
# [2,]    2    3    4
# [3,]    3    4    5
# [4,]    4    5    6

Man kann die Zahlenwerte in dieser Matrix auch so verstehen:

Die 1 steht für den konstanten Term, der dadurch entsteht, dass die beiden konstanten Terme der Eingabe-Polynome multipliziert werden. (Er steht also für die Potenz 0 im Ergebnis-Polynom.)
Die 6 gehört zur höchsten Potenz des Ergebnis-Polynoms, stimmt aber nicht damit überein: die höchste Potenz ist hier 5 und allgemein ist die höchste Potenz um 1 geringer als die Länge des Koeffizienten-Vektors. Die höchste Potenz im Ergebnis-Polynom entsteht bei der Multiplikation der beiden höchsten Potenzen in den Eingabe-Polynomen. Folglich gilt für den Grad des Polynoms grad(f · g) = grad(f) + grad (g) und somit für die Länge lgth des Produkt-Polynoms f · g: lgth <- length(f) + length(g) - 1 .

Die beiden Matrizen m (die Multiplikationstabelle der Polynome) und secondaries (die Nebendiagonalen) haben identische Dimensionen. Daher ist es leicht, zu einem gegebenen Wert von secondaries zuerst die Indizes zu bestimmen, wo diese Wert angenommen werden und anschließend aus der Multiplikationstabelle diese Indizes anzusprechen:

m[which(secondaries == 1)]
# [1]  3
m[which(secondaries == 2)]
# [1]  0 -2
m[which(secondaries == 3)]
#[1] -9  0  1
m[which(secondaries == 4)]
# [1] 3 6 0
m[which(secondaries == 5)]
# [1] -2 -3
m[which(secondaries == 6)]
# [1] 1

Von diesen Komponenten muss man jetzt nur noch die Summen bilden und hat somit den Vektor der Koeffizienten für das Ergebnis-Polynom; und die Berechnung lässt sich leicht in eine Schleife verpacken:

lgth <- length(f) + length(g) - 1
product <- vector(mode = "numeric", length = lgth )

for( i in (1:lgth) ) {
  product[i] <- sum( m[which(secondaries == i)] )
}

product
# [1] 3 -2 -8  9 -5  1

3. Schritt:

Die Berechnungen werden in eine Funktion mult.polynomial() verpackt. Diese Funktion erhält als Eingabewerte zwei Vektoren für die Koeffizienten der Polynome und sie erzeugt als Rückgabewert einen Vektor mit den Koeffizienten des Ergebnis-Polynoms:

mult.polynomial <- function(a, b){
  secondaries <- outer(X = seq_along(along.with = a) - 1, Y = seq_along(along.with = b), FUN = "+")
  m <- a %o% b
  lgth <- length(a) + length(b) - 1
  product <- vector(mode = "numeric", length = lgth )
  
  for( i in (1:lgth) ) {
    product[i] <- sum( m[which(secondaries == i)] )
  }
  return(product)
}

mult.polynomial(a = f, b = g)
# [1]  3 -2 -8  9 -5  1

Zum Test wurden die beiden Vektoren f und g von oben eingesetzt.

Aufgaben:

1. Welche Polynome und welche Folgen von Koeffizienten werden durch das folgende Skript erzeugt:

n <- 20
a <- c(1, 1)
p <- a

for(i in (1:n)){
  p <- mult.polynomial(a, p)
  print(p)
}

2. Welche Operation wird durch folgende Funktion mit einem Polynom a ausgeführt:

i <- function(a, c = 1){
  lgth.a <- length(a)
  result <- vector(mode = "numeric", length = lgth.a - 1)
  result[1] <- c
  for(i in ((1:lgth.a))){
    result[i + 1] <- a[i] / i
  }
  return(result)
}