# # 1. Sitzung des Stochastik-Praktikums, WS 2020/2021, JGU Mainz # 2.11.2020 # # ######################################################### # 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() # z.B. speziell demo(graphics) # # 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 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() # ein neues Grafikfenster oeffnen # (je nach System unterschiedlich, abrufbar mit options("device") # zeigt an, welche Grafik"geraete" das laufende R kennt 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 q0 <- scan("quadratzahlen.txt") q0 # Für "richtige" Daten häufig sinnvoller: # data.frames mit read.table einlesen weinsortenRLP <- read.table("Wein_RLP.csv", header=TRUE, sep=",") # schauen wir die Daten an weinsortenRLP View(weinsortenRLP) # Auf Spalten eines data.frame zugreifen weinsortenRLP$Ahr weinsortenRLP[,3] # alternativ # Spalten als Variable verfügbar machen attach(weinsortenRLP) Ahr # Histogramm zeichnen hist(Rheinland.Pfalz.gesamt, main="Weinsorten in RLP", xlab="Anbaufläche (ha)", ylab="Anzahl", col="green") # Man kann einen sog. "Rug" hinzufuegen, um auch die Einzelwerte # darzustellen rug(Rheinland.Pfalz.gesamt) # Klassenzahl/-groessen manipulieren hist(Rheinland.Pfalz.gesamt, main="Weinsorten in RLP", xlab="Anbaufläche (ha)", ylab="Anzahl", col="green", breaks=25) rug(Rheinland.Pfalz.gesamt) # Tortendiagramm rote <- 28:44 pie(Rheinland.Pfalz.gesamt[rote], labels=Rebsorte[rote], main="Rotweinanbau in RLP") # Balkendiagramm barplot(Rheinhessen[28:32], names.arg=Rebsorte[28:32], col="darkred", cex.names=1.1, 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=",") # (beachte: "deutsches" Dezimaltrennzeichen [=Komma] in der Original-Datei) ertraege View(ertraege) attach(ertraege) # Numerische Kenngroessen: # empirische Mittelwerte, Mediane, empir. Streuung der # Jahreserträge einzelner Rebsorten mean(Riesling) mean(Bacchus) mean(Portugieser_Blauer) median(Riesling) sd(Riesling) # Beobachtete Jahresertraege 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") # Jahresertraege 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") # (durchgezogene Striche statt einzelner Punkte) # vergleichen wir zwei verschiedene Sorten: plot(Jahr, Elbling, type="l", col="blue") points(Jahr, Riesling, type="l", col="darkgreen") # fuegen wir noch eine Legende hinzu: legend("topright", legend=c("Elbling","Riesling"),col=c("blue","darkgreen"),lty=1) detach(ertraege) # ######################################################### # 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) # Ueberzeugen wir uns, dass dasselbe ausgerechnet wurde: summen-summen.alternativ max(abs(summen-summen.alternativ)) # Uebrigens: 1+4+...+(n-1)^2+n^2=n(n+1)(2n+1)/6 # Fuege 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 (ca. 10-75 n.Chr (?)) # [auch "babylonisches Wurzelziehen" genannt] # 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) # Fuege 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 }