Jakob Richter

Workflow für Simulationsstudien

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.
Read the rest of this entry »

Wiener Prozess

Drei Methoden um mit R den Wiener Prozess zu simulieren:

Normalverteilte Zuwächse

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
set.seed(44147)
T <- 1 #Zeitintervall
Ns <- c(10,100,1000)
 
plotres <- lapply(Ns, function(N){
	Delta <- T/N # Zeitschritte
	W <- c(0, cumsum( sqrt(Delta) * rnorm(N)))
	t <- seq(0,T, length.out=N+1)
	res <- cbind.data.frame(W=W,t=t,N=N)
	return(res)
})
 
plotres <- do.call("rbind",plotres)
plotres <- as.data.frame(plotres)
plotres$N <- as.factor(plotres$N)
 
library(ggplot2)
g <- ggplot(plotres,aes(x=t,y=W,color=N))
g + geom_line(size=0.7)

Random Walk

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
set.seed(44137)
T <- 1
ns <- c(10,100,1000)
 
randres <- lapply(ns,function(n){
	t <- seq(0,T,length.out=n+1)
	S <- cumsum(sample(x=c(-1,1),size=n,replace=T))
	W <- c(0,S)
	W <- W * sqrt(T/n) #hier liegt ein Fehler im Buch vor.
	cbind.data.frame(W=W,t=t,N=n)
})
 
randres <- do.call("rbind",randres)
randres <- as.data.frame(randres)
randres$N <- as.factor(randres$N)
 
library(ggplot2)
h <- ggplot(randres,aes(x=t,y=W,color=N))
h <- h + geom_line(size=1)
 
#Punkte an den Sprungstellen fuer kleine N
datpoints <- randres[as.vector(randres$N)<=100,]
h + geom_point(data=datpoints,aes(x=t,y=W,color=N),size=3)

Karhunen-Loève Approximation

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
set.seed(44137)
 
phi <- function(i,t,T){
	(2 * sqrt(2*T)) / ((2*i+1)*pi) * sin(((2*i+1)*pi*t)/(2*T))
}
 
T <- 1
ns <- c(10,50,100) #sollten vielfache voneinander sein.
Zs <- rnorm(max(ns))
 
karres <- lapply(ns, function(n){
	Gen <- n*10
	t <- seq(0,T,length.out=Gen+1)
	#Waehlen immer die mittlere Zufallszahl aus.
	mitte <- ceiling(seq(from=max(ns)/(2*n), to=max(ns), by=max(ns)/n))
	Z <- Zs[mitte]
	W <- numeric(Gen+1)
	for(i in 2:(Gen+1)){
		W[i] <- sum(Z*sapply(1:n, function(x) phi(x,t[i],T)))
	}
	cbind.data.frame(W=W,t=t,N=n)
})
 
karres <- do.call("rbind",karres)
karres <- as.data.frame(karres)
karres$N <- as.factor(karres$N)
 
library(ggplot2)
i <- ggplot(karres,aes(x=t,y=W,color=N))
i + geom_line(size=1)

Für weiter interessierte hier die Folien des Vortrags:
01_Brownsche_Bewegung_und_drei_Methoden_zur_Simulation.pdf