# # Sitzung des Stochastik-Praktikums, WS 2013/2014, JGU Mainz # 21.10.2013 # # ######################################################### # Grundlegende R-Benutzung # ######################################################### # R als Taschenrechner 1+2 sqrt(7) sin(2.5) cos(pi/2) # Variablenzuweisung: x <- 5 x # alternativ (aber eher nicht gebräuchlich): x=5 # andere Typen von Variablen: y <- "Zeichenkette" z <- TRUE w <- -2+7.5i y; z; w # Vergleiche: 1 < 2 y == "Zeichenkette" x != 4 x <= 3 x <= 3 || y == "Zeichenkette" x <= 3 && y == "Zeichenkette" # Vektor erzeugen: v <- c(1,2,3.1415,-17) # Elementzugriffe: v v[1]; v[3] v[2] <- -2 v # Zugriff auf ein Element ausserhalb des # definierten Bereichs liefert # "NA"="not available" v[5] # R vergrößert Vektoren dynamisch, falls notwendig: v[5] <- 7 v v[10] <- 0.55 v # die meisten mathematischen Operationen auf Vektoren # führt R "elementweise" aus und "rezykliert" gegebenenfalls, # wenn die Operanden verschiedene Länge haben w <- c(1,2,3.1415,-17) w2 <- c(1,2) w+w2 w*w2 w*2 # # Hilfe erhalten: help(Befehl) oder ?Befehl # help(seq) ?seq # Hilfe in einem Browserfenster öffnen: # (je nach System ggf. einen anderen Browser eingeben bzw. # einfach help.start() ) help.start(browser="firefox") # # Mehr über Vektoren # # "fortgeschrittenere" Arten, Vektoren zu erzeugen: seq(from=1, to=10, by=3) rep(1,5) 2:12 # Vektor-Indizierung mit Vektoren, Teilauswahl v <- 2*(1:10) v v[1:4] v<6 v[v<6] # # Listen # l <- list(1.2, "text", -5, FALSE) l l[[1]] l[[4]] # Indizierung ausserhalb der Grenzen führt bei Listen # zu Fehlermeldung: l[[5]] # Listenelemente können benannt sein, dann ist Zugriff # via Name möglich li <- list(Zahl=1.2, Wort="text", Dings=-5, Wert=FALSE) li li$Zahl li$Wert <- TRUE li # # R kennt auch Matrizen # z <- 1:8 m <- matrix(z, nrow=2, ncol=4) m m[1,3] m[,4] m[2,] # # Rs Demos # demo() # # Uniform auf [0,1] verteilte (Pseudo-)Zufallszahlen # runif(1) runif(12) # # Skript einlesen # Übrigens: das Arbeitsverzeichnis abfragen/wechseln mit # getwd(), setwd(...) # source("beispielskript.R") # # Ausgabe in Datei umleiten # sink(file="sink-bsp.dat") sin(2*pi*(0:11)/12) sink() # Ausgabe zurück auf Konsole sin(2*pi*(0:11)/12) # # Grafikausgaben # a <- 1:10 plot(a) plot(a, main="Die Zahlen von 1 bis 10", col="blue") plot(a, main="Die Zahlen von 1 bis 10", col="blue", ylab="Zahl", xlim=c(0,20), ylim=c(-2,12)) x11() # neues Grafikfenster öffnen absz <- cos(2*pi*(0:11)/12) ordi <- sin(2*pi*(0:11)/12) plot(absz, ordi) plot(absz, ordi, type="l") plot(absz, ordi, type="l", lwd=3) plot(absz, ordi, type="l", lwd=3, lty=2) points(absz, ordi, col="tomato2") # Punkte zu einem Plot hinzufügen # Grafiken als Postscript-Datei erzeugen: postscript("Beispielbild1.ps") plot(a, main="Die Zahlen von 1 bis 10", col="blue") dev.off() # oder ein neues Grafikgerät öffnen, um das alte zu beenden # Analog: Grafiken als PDF-Datei erzeugen: pdf("Beispielbild1.pdf") plot(a, main="Die Zahlen von 1 bis 10", col="blue") dev.off() # oder ein neues Grafikgerät öffnen, um das alte zu beenden # Grafikparameter setzen x11() par(mfrow=c(2,2)) # Hier beispielsweise: Mehrere Plots auf einmal plot(absz, ordi) plot(absz, ordi, type="l") plot(absz, ordi, type="l", lwd=3) plot(absz, ordi, type="l", lwd=3, lty=2) points(absz, ordi, col="tomato2") # Punkte zu einem Plot hinzufügen # aktuelle Grafik in eine Datei kopieren: dev.copy2pdf(file="Beispielbild2.pdf") # (analog: dev.copy2eps) dev.off() # ######################################################### # Daten einlesen, einfache grafische Analysen # ######################################################### # Elementar: scan q <- scan("quadratzahlen.txt") q # Für "richtige" Daten häufig sinnvoller: # data.frames mit read.table einlesen #weinsortenRLP <- read.table("Wein_RLP_2012.txt", row.names=1, header=TRUE) weinsortenRLP <- read.table("Wein_RLP_2012.txt", header=TRUE) # Auf Spalten eines data.frame zugreifen weinsortenRLP$Ahr weinsortenRLP[,3] # alternativ # Spalten als Variable verfügbar machen attach(weinsortenRLP) Ahr # Histogramm zeichnen hist(RLP_ges, main="Weinsorten in RLP", xlab="Anbaufläche (ha)", ylab="Anzahl", col="green") # Man kann einen sog. "Rug" hinzufügen, um auch die Einzelwerte # darzustellen rug(RLP_ges) # Klassenzahl/größen manipulieren hist(RLP_ges, main="Weinsorten in RLP", xlab="Anbaufläche (ha)", ylab="Anzahl", col="green", breaks=25) # Tortendiagramm rote <- 25:40 pie(RLP_ges[rote], labels=Rebsorten[rote], main="Rotweinanbau in RLP") # Balkendiagramm barplot(Rheinhessen[25:30], names.arg=Rebsorten[25:30], col="darkred", cex.names=0.7, main="Anbauflächen einiger Rotweinsorten in Rheinhessen", ylab="Anbauflächen (ha)") detach(weinsortenRLP) # Ein anderer Datensatz: ertraege <- read.table("Wein_Ertraege.txt", header=TRUE, dec=",") attach(ertraege) # Numerische Kenngrößen: # empirische Mittelwerte, Mediane, empir. Streuung der # Jahreserträge einzelner Rebsorten mean(Riesling) mean(Bacchus) mean(Portugieser_Blauer) median(Riesling) sd(Riesling) # Beobachtete Jahreserträge als "Stripchart" # visualisieren: stripchart(Riesling) # Mehrere Stripcharts auf einmal: stripchart(list("Riesling"=Riesling, "Bacchus"=Bacchus, "Blauer Portugieser"=Portugieser_Blauer), xlab="Ertrag in hl pro ha") # "jitter" verschiebt (hier) die y-Koordinaten der Werte # zufällig, macht "Bindungen" (=gleiche oder nahezu gleiche Werte) # (besser) sichtbar: stripchart(list("Riesling"=Riesling, "Bacchus"=Bacchus, "Blauer Portugieser"=Portugieser_Blauer), method="jitter", xlab="Ertrag in hl pro ha") # Jahreserträge einzelner Rebsorten als Boxplots boxplot(Riesling) boxplot(list(Riesling, Bacchus, Portugieser_Blauer)) boxplot(list(Riesling, Bacchus, Portugieser_Blauer), names=c("Riesling", "Bacchus", "Blauer Portugieser")) # Ertrag als Funktion des Jahrs: plot(Jahr, Bacchus) plot(Jahr, Bacchus, type="l") # ######################################################### # Beispiele einfacher Schleifen und Funktionen # ######################################################### # # Summe der ersten n Quadratzahlen, für n=1,..,N: # N <- 100 # Einmal "von Hand": summen <- numeric(N) # Erzeuge "leeren" Vektor z <- 0 for (i in 1:N) { z <- z+i^2 summen[i] <- z } plot(summen, main="Summe der ersten n Quadratzahlen", ylab="", xlab="n") # und einmal mit Rs Vektor-Operationen: summen.alternativ <- cumsum((1:N)^2) # Überzeugen wir uns, dass dasselbe ausgerechnet wurde: summen-summen.alternativ max(abs(summen-summen.alternativ)) # Übrigens: 1+...+n=n(n+1)(2n+1)/6 # Füge den Funktionsgraphen von x(x+1)(2x+1)/6 der Grafik hinzu curve(x*(x+1)*(2*x+1)/6, col="blue", add=TRUE) # # Funktionsdefinitionen # f <- function(a,b) { a-b } f(3,5) # liefert 3-5=-2 # # Approximative Quadratwurzelberechnung nach Heron von Alexandria # genauigkeit <- 1e-7 heronapprox <- function(a) { b <- 1 while(abs(b^2-a)>genauigkeit) { b<-(b+a/b)/2 } b } # Beispiel: heronapprox(7) # Zum Vergleich: sqrt(7) # Heron-Approximation als Funktionsiteration auffassen # (und diese Auffassung illustrieren): a <- 10 # Setze "leeres" Zeichenfenster mit Koordinatensystem auf # (type="n") plot(0, xlim=c(0,2*sqrt(a)), ylim=c(0,2*sqrt(a)), type="n", xlab="x", ylab="y", main=paste("Heron-Approx. von sqrt(",a,") als Funktionsiteration")) # Zeichne Funktion ein f <- function(x) {(x+a/x)/2} curve(f, col="blue", lwd=3, from=0.01, to=2*sqrt(a), add=TRUE) # Füge Gerade (mit Achsenabschnitt 0 und Steigung 1, d.h. Diagonale) ein abline(0,1) schritte <- 4 b <- 1 lines(c(b,b), c(-0.2,f(b)), col="darkgray") for (i in 1:schritte) { bneu <- f(b) lines(c(b,bneu,bneu), c(bneu,bneu,f(bneu)), col="darkgray") b<-bneu }