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.
Noch keine Stimmen abgegeben
Noch keine Kommentare

Einordnung des Artikels

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:

  1. Die Funktion mapply(), die zahlreiche Eingabewerte besitzt und somit mehrere Funktionalit├Ąten bereitstellt.
  2. 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.
  3. 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.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.

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

  1. Ist es ├╝berhaupt m├Âglich im Argument ... von mapply() Dataframes einzugeben?
  2. 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.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 (a0, a1, ..., an) 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 ai, i = 0, ..., n, wobei man ├╝blicherweise die erste Konstante mit 0 indiziert, da sie zur Potenz x0 = 1 geh├Ârt. Wird dies sp├Ąter implementiert, werden die Konstanten ai 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.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:

x0, x1, ..., xn-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 + x2,

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 + x2

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:

  1. 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.
  2. 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).Abbildung 4: Darstellung der Exponentialfunktion und ihrer Approximation durch die Taylor-Entwicklung um x = 0 (mit Polynomen vom Grad 0, 2, ..., 10).

Aufgaben:

  1. ├ť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.
  2. Berechnen Sie die maximale Differenz zwischen den Werten der Exponentialfunktion und der jeweiligen Approximation. Bei welchem x-Wert ist der Betrag der Differenz maximal?
  3. In den Abbildungen 5 und 6 wird die Gau├čsche Glockenkurve f(x) = exp(-x2) 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 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).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:

  1. Die lineare N├Ąherung f├╝r exp(x) ┬Ě exp(y) in Abbildung 7,
  2. die quadratische N├Ąherung f├╝r exp(x) ┬Ě exp(y) in Abbildung 8,
  3. 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 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 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.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.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).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:

  1. Die Multiplikationstabelle der Koeffizienten wird erstellt.
  2. Die Komponenten auf den Nebendiagonalen werden addiert.
  3. 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)
}