Jakob Richter

Eine häufige Aufgabe ist es, für verschiedene Parameter ein und die selbe Simulation mit R durchzuführen und dann die Ergebnisse auszuwerten. Persönlich bevorzuge ich es, erst alles zu simulieren und dann auszuwerten. Dieses Vorgehen hat mehrere Vorteile: Bei langandauernden Simulationen müssen diese nur einmalig ausgeführt werden. Außerdem: Wenn ich mit meiner Auswertung unzufrieden bin, kann ich einfach die Auswertungsprozedur modifizieren ohne noch einmal alles durchlaufen lassen zu müssen und der Quelltext ist auch nachvollziehbarer. Insbesondere ist das Vorgehen aber auch angepasst auf die grafische Darstellung und Auswertung mit ggplot2.

Darum ist es immer mein Ziel alle Ergebnisse der Simulation in einem data.frame zu speichern, wobei in jeder Zeile des data.frames der simulierte Wert und seine Parameter mit denen er erzeugt wurde stehen sollen.

Beispiel

1
2
3
4

parameterA <- c(10,50,100)
parameterB <- c(0.3,0.9)
parameterC <- round(c(0.1,0.2)*pi,3)
parameter <- expand.grid(A=parameterA,B=parameterB,C=parameterC)

Hier treffen wir schon auf die erste tolle Funktion expand.grid(). Sie erzeugt einen data.frame der alle möglichen Kombinationen von den Werten in den Eingabevektoren erstellt. Natürlich geht das auch mit mehr als dreien und auch mit Zeichenketten anstelle von Zahlen.

Kommen wir nun zur Simulation:

6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

set.seed(234)
simerg <- apply(parameter,1,function(pars){
	A <- pars[1]
	B <- pars[2]
	C <- pars[3]
	xycord <- cbind.data.frame(x=0,y=0)
	richtung <- 0
	for(i in 2:A){
		neuerichtung <- richtung + C * rnorm(1,1,i/A*2)
		richtung <- neuerichtung
		schrittweite <- rnorm(1,mean=B,sd=B)
		xneu <- xycord[i-1,"x"] + schrittweite*sin(richtung)
		yneu <- xycord[i-1,"y"] + schrittweite*cos(richtung)
		xycord[i,] <- cbind.data.frame(x=xneu,y=yneu)
	}
	#Ausgabe als data.frame sehr nützlich
	#Netterweise werden A und B über alle Zeilen wiederholt
	cbind.data.frame(Steps=A,Weite=B,Drift=C,xycord) #letzte Zeiel wird automatisch zurückgegeben
})
#Warnings können ignoriert werden

Die nächste Funktion, die kein Fremdbegriff für R-Nutzer sein darf: apply(). Angewendet auf einen data.frame kann sie über die Zeilen apply(...,MARGIN=1,...) oder über die Spalten apply(...,MARGIN=2,...) laufen. Wenn bspw. MARGIN=1 wird jeweils jede Zeile als Vektor (also in dem Fall transponiert) an die in apply() angegebene Funktion übermittelt. Die Funktion können wir unter Verwendung von {...} übersichtlich über mehrere Zeilen schreiben, wie man es gewohnt ist.
Damit wir am Ende alle Simulationsergebnisse gut zusammenfassen können ist es von Vorteil diese gleich in einem data.frame zusammenzupacken. Mit cbind.data.frame() können wir die Faktoren und die Ergebnisse zusammenfassen und den Spalten gleichzeitig schon jetzt bedeutsame Namen geben (hier A,B und Z).
Schauen wir uns das Ergebnis mal an:

28
29
30
31
32
33
34
35
36
37
38

simerg
# [[1]]
#    Steps Weite Drift          x           y
# 1     10   0.3 0.314  0.0000000  0.00000000
# 2     10   0.3 0.314 -0.1221397 -0.29132710
# ...
# [[2]]
#    Steps Weite Drift          x           y
# 1     50   0.3 0.314 0.00000000  0.00000000
# 2     50   0.3 0.314 0.09795731  0.28858860
# ...

Da die apply() hier einen data.frame zurückgegeben hat, war sich R unsicher und hat jedes Ergebnis in eine Liste gesteckt (Das macht R nicht standardmäßig, aber wenn es “denkt”, dass es damit den sicheren Weg geht, denn in ein einzelnes List-Item kann alles rein – unabhängig davon, was in den anderen Listeneinträgen enthalten ist.). Wir sehen natürlich, dass das blöd ist, denn alle Listeneinträge haben die gleiche Struktur (gleiche Spaltenanzahl und gleiche Spaltennamen) und passen gut untereinander in den heiß geliebten data.frame.
Dafür gibt es einen kurzen und genialen Trick:

39

simulationsergebnis <- do.call("rbind",simerg) #Funktioniert natürlich nur, wenn alle data.frames in der Liste die gleichen Spaltennamen haben. Die Anzahl der Zeilen kann glücklicherweise unterschiedlich sein!

Jetzt ist simulationsergebnis ein data.frame und kinderleicht auszuwerten.

41
42
43
44
45
46
47
48
49
50
51
52

#folgendes ist für ggplot sinnvoll:
simulationsergebnis$Steps <- as.factor(simulationsergebnis$Steps)
simulationsergebnis$Weite <- as.factor(simulationsergebnis$Weite)
simulationsergebnis$Drift <- as.factor(simulationsergebnis$Drift)
 
#Auswertung der Simulation
library("ggplot2")
g <- ggplot(data=simulationsergebnis, aes(x=x,y=y,color=Steps))
g + geom_path() + facet_wrap(Weite~Drift,scales="free")
 
#Zahlenauswertung
aggregate(x~Steps+Weite+Drift,data=simulationsergebnis,FUN=range)

Zum Verstehen kann ich nur raten alles einmal auszuführen und sich den Inhalt der Variablen nach jedem Schritt anzuschauen und alles mal mit anderen Werten und eigenen Ideen anzureichern.