Programmieren mit Statistischer Software
Eugster
7. Block
SS2012
7 Aufgaben
Aufgabe 1:
Beim Zahlenraten erzeugt der Computer eine Zufallszahl in einem bestimmten Bereich (z.B. [0, 20])
und Sie müssen diese Zahl erraten. Bei jedem Versuch bekommen Sie die Meldung, ob die gesuchte Zahl
größer oder kleiner als die eingegebene Zahl ist; oder ob Sie die Zahl gefunden haben. Ziel bei diesem
Spiel ist es, möglichst wenige Versuche zu brauchen.
Programmieren Sie die Funktion
> zahlenraten <- function(min = 0, max = 20) {
+
...
+ }
welches dieses Spiel umsetzt. Die beiden Argumente min und max definieren den Bereich in dem die
Zufallszahl vom Computer gezogen wird.
Das Einlesen einer Zahl funktioniert mit:
> versuch <- as.numeric(readline("Neuer Versuch: "))
Lösung:
Bei diesem Spiel ist die Anzahl der Schleifendurchläufe nicht bekannt, d.h. eine while-Schleife wird
verwendet. In jedem Schleifendurchlauf wird der Spieler nach einer Zahl gefragt, welche dann entsprechend
ausgewertet wird. Hat der Spieler die Zahl erraten, wird eine logische Variable zahlErraten auf TRUE
gesetzt, was den Abbruch der Schleife bedeutet.
> zahlenraten <- function(min = 0, max = 20) {
+
stopifnot(min < max)
+
+
zahl <- sample(seq(min, max), size=1)
+
zahlErraten <- FALSE
+
i <- 0
+
+
cat("Erraten Sie eine Zahl zwischen", min, "und", max, "...\n")
+
while (!zahlErraten) {
+
i <- i + 1
+
+
versuch <- as.numeric(readline(paste(i, ". Versuch: ", sep="")))
+
+
if (versuch == zahl ) {
+
zahlErraten <- TRUE
1
+
}
+
else {
+
if ( versuch < zahl )
+
cat("Zahl ist größer!\n")
+
else
+
cat("Zahl ist kleiner!\n")
+
}
+
}
+
+
cat("Zahl erraten! Sie haben", i, "Versuche benötigt!\n")
+ }
Eine mögliche Spielfolge ist:
> zahlenraten()
Erraten Sie eine Zahl zwischen 0 und 20 ...
1. Versuch: 10
Zahl ist größer!
2. Versuch: 15
Zahl ist größer!
3. Versuch: 17
Zahl ist größer!
4. Versuch: 18
Zahl ist größer!
5. Versuch: 19
Zahl erraten! Sie haben 5 Versuche benötigt!
Anmerkung: Die optimale Strategie für dieses Spiel ist die binäre Suche.
Aufgabe 2:
Die Kreiszahl π kann mit Hilfe eines Monte-Carlo-Algorithmus näherungsweise bestimmt werden: für
die Berechnung lässt man zufällige Punkte auf ein Quadrat regnen und berechnet, ob sie innerhalb oder
außerhalb eines eingeschriebenen Kreises liegen. Der Anteil der innen liegenden Punkte ist gleich π4 . Mit
steigender Anzahl der Punkten steigt auch die Genauigkeit.
2
Schreiben Sie eine Funktion, welche für eine gegebene Punkteanzahl n π berechnet.
(Hinweis: Für xy-Paare, die innerhalb des Viertelkreises liegen, gilt x2 + y 2 ≤ 1.)
Lösung:
Funktion zur Bestimmung mit Hilfe des Monte-Carlo-Algorithmus:
> mcpi <- function(n) {
+
x <- runif(n)
+
y <- runif(n)
+
+
## xy-Paare die innerhalb des Viertelkreises liegen:
+
inside <- ((x^2 + y^2) <= 1)
+
+
## Pi berechnen:
+
pi <- 4 * sum(inside) / n
+
+
return(pi)
+ }
Aufruf für n = 100:
> set.seed(1234)
> mcpi(100)
[1] 3.32
Simulation um zu sehen, wie sich die Schätzung dem echten π annähert:
> set.seed(1234)
> s <- seq(0, 5000000, by = 10000)[-1]
> mypi <- sapply(s, mcpi)
Visualisierung:
3
>
>
>
>
op <- par(mar = c(4, 4, 0, 0), ps = 8, las = 2)
plot(s, mypi, type = "l", lwd = 2)
abline(h = pi, col = 2, lwd = 2)
par(op)
3.15
mypi
3.14
3.13
3.12
5e+06
4e+06
3e+06
2e+06
1e+06
0e+00
3.11
s
Aufgabe 3:
Bei einer gegebenen Menge von Datenpunkten x1 , x2 , . . . , xn ist der gleitende Mittelwert (moving
average) x̄i der Mittelwert von w (w eine ungerade Zahl größer 0) Datenpunkten, genauer w−1
2 vor xi
w−1
und w−1
nach
x
liegend.
Daraus
folgt,
dass
für
die
ersten
und
letzen
Datenpunkte
kein
gleitender
i
2
2
Mittelwert definiert ist.
Schreiben Sie eine Funktion
> movave <- function(x, w) { ... }
welche einen Vektor x von Datenpunkten und eine Fenstergröße w (ganze ungerade Zahl größer 0) als
Argumente übernimmt. Rückgabewert soll der Vektor der gleitenden Mittelwerte sein, der die selbe Länge
wie x besitzt mit NA bei den ersten und letzten w−1
2 Elementen.
Laden Sie dann den Datensatz LakeHuron,
> data("LakeHuron", package = "datasets")
> x <- as.vector(LakeHuron)
berechnen Sie den gleitenden Mittelwert für w = 3, 5, 21 und zeichnen Sie diese in eine Grafik zusammen
mit den Originaldaten ein.
Lösung:
4
> movave <- function(x, w) {
+
n <- length(x)
+
w2 <- (w - 1) / 2
+
ma <- rep(NA, n)
+
+
for ( i in seq(w2+1, n-w2) )
+
ma[i] <- mean(x[(i-w2):(i+w2)])
+
+
ma
+ }
Lake Huron:
>
>
>
>
>
data("LakeHuron", package = "datasets")
x <- as.vector(LakeHuron)
ma3 <- movave(x, 3)
ma5 <- movave(x, 5)
ma21 <- movave(x, 21)
>
>
>
>
>
>
op <- par(mar = c(4, 4, 0, 0), ps = 8, las = 2)
plot(x, type = "l")
lines(ma3, col = 2)
lines(ma5, col = 3)
lines(ma21, col = 4)
par(op)
582
581
x
580
579
578
577
100
80
60
40
20
0
576
Index
Aufgabe 4:
In der Einheit “Funktionen” haben Sie ein Funktion which maxdev() entwickelt, welche den Index der
Beobachtung mit der maximalen Abweichung vom Median aller Beobachtungen zurückliefert.
Schreiben Sie eine Funktion
5
> which_maxdevs <- function(...) {
+
...
+ }
welche diese Funktion für alle Spalten einer Matrix oder eines data.frames ausführt.
Installieren Sie das Paket benchmark und laden Sie den Datensatz uci621raw. Schauen Sie sich die
Struktur des Objektes an und extahieren Sie die Daten für die das erste Element der dritten und das
erste Element der vierten Dimension. Wenden Sie die programmierte Funktion auf diese Teilmatrix an.
Lösung:
> which_maxdev <- function(x, na.rm = TRUE) {
+
stopifnot(is.numeric(x) & is.vector(x))
+
stopifnot(is.logical(na.rm))
+
+
mdn <- median(x, na.rm = na.rm)
+
devs <- abs(x - mdn)
+
which.max(devs)
+ }
> which_maxdevs <- function(m, na.rm = TRUE) {
+
stopifnot(is.matrix(m) | is.data.frame(m))
+
apply(m, 2, which_maxdev)
+ }
> data("uci621raw", package = "benchmark")
> str(uci621raw)
num [1:250, 1:6, 1:2, 1:21] NA 0.0618 0.0445 0.0413 0.0685 ...
- attr(*, "dimnames")=List of 4
..$ samp: NULL
..$ alg : chr [1:6] "lda" "rf" "knn" "rpart" ...
..$ perf: chr [1:2] "Misclassification" "Time"
..$ ds : chr [1:21] "BreastCancer" "Cards" "chess" "Circle" ...
>
+ dat <- uci621raw[, , 1, 1]
>
> which_maxdevs(dat)
lda
40
rf
151
knn rpart
185
181
svm
216
nnet
106
6
© Copyright 2024 ExpyDoc
