# # Statistik fuer Informatiker, SS 2018, JGU Mainz # # ############################################### # 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() # ######################################################### # 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+4+...+(n-1)^2+n^2=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 (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) # 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 } # ############################################### # Knappe Anmerkungen zu Histogrammen # (ein Mini-Ausflug in die deskriptive Statistik) # ############################################### # # n gegebene Zahlen x_1,...,x_n ("Beobachtungswerte") # # zerlege in k Klassen gemäß Zerlegungsintervallen # (z_0,z_1], (z_1,z_2],...,(z_{k-1},z_k], # # sei h_k := # { i <= n : z_{k-1} < x_i <= z_k } # (Häufigkeit der k-ten Klasse), # wir nehmen an, dass z_0 kleiner als der # kleinste Wert der x_i und z_k größer oder gleich # dem größten Wert der x_i ist. # # Ein Histogramm (in seiner einfachsten Form) # besteht dann aus k Balken, # wobei Balken j auf der x-Achse zwischen # z_{j-1} und z_j aufsitzt und Höhe h_j hat # (eine gute Wahl der Zerlegungsintervalle ist # eine gewisse Kunst für sich, R verwendet # standardmäßig "plausible Heuristiken", # siehe ?hist und ?pretty ) # z.B. x <- scan("beispieldaten.txt",sep=",") # man sieht ziemlich wenig, wenn man nur die # "Rohdaten" als Zahlenliste anschaut: x # hist(x) # erzwingen wir die Klassengrenzen von "von Hand" hist(x, breaks=seq(1.0,4.0,by=0.25)) # es kann (gerade zu Vergleichszwecken) # manchmal günstiger sein, die Balkenhöhen so zu # skalieren, dass die Gesamtfläche der Balken = 1 # ergibt # (man spricht dann auch von einem "Dichtehistogramm") hist(x, breaks=seq(1.0,4.0,by=0.25), prob=T) # Titel und Beschriftungen hinzufügen: hist(x, breaks=seq(1.0,4.0,by=0.25), prob=T, main="Histogramm der Beispieldaten", xlab="x-Werte", ylab="Dichte") # ############################################### # Bonusmaterial: # Making of fuer die Startbeispiel-Folien # ############################################### # Parametrisierung der Randkurve von B radiusfkt <- function(w) { w <- w+pi/7 (sin(w)*sqrt(abs(cos(w)))/(sin(w)+7/5)+2-2*sin(w))/7 } t <- 2*pi*seq(0,1,by=0.005) r2 <- sapply(t, radiusfkt) plot(r2*cos(t)+0.5,r2*sin(t)+0.6, xlim=c(0,1),ylim=c(0,1),type='b') rect(0,0,1,1) drin <- function(x,y) { # entscheide, ob Punkt (x,y) innerhalb der Figur xx <- x-0.5; yy <- y-0.6 ra <- sqrt(xx^2+yy^2) if (ra==0) return(TRUE) wi <- sign(xx)*asin(yy/ra)+pi*sign(yy)*ifelse(xx<0,1,0) if (ra <= radiusfkt(wi)) TRUE else FALSE } for (i in 1:1000) { x <- runif(1); y <- runif(1) points(x,y, col=ifelse(drin(x,y),'red','green')) } repl <- 10000 erf <- 0 for (i in 1:repl) { if (drin(runif(1),runif(1))) { erf <- erf+1 } } erf/repl # zeichne Einheitsquadrat und die # Menge B radiusfkt <- function(w) { w <- w+pi/7 (sin(w)*sqrt(abs(cos(w)))/(sin(w)+7/5)+2-2*sin(w))/7 } t <- 2*pi*seq(0,1,by=0.001) r2 <- sapply(t, radiusfkt) # oeffne ein "leeres" Grafikfenster: plot(0, xlim=c(0,1),ylim=c(0,1), xlab='', ylab='',type='n') # zeichne Einheitsquadrat ein: rect(0,0,1,1) # zeichne B ein: polygon(r2*cos(t)+0.5,r2*sin(t)+0.6,col='lightblue') points(r2*cos(t)+0.5,r2*sin(t)+0.6, xlim=c(0,1),ylim=c(0,1),type='l',lwd=2) # waehle zufaellig Punkt (x,y) # im Einheitsquadrat x <- runif(1); y <- runif(1) # und zeichne ein: points(x,y, pch=20, col=ifelse(drin(x,y),'darkblue','grey')) # oeffne ein "leeres" Grafikfenster: plot(0, xlim=c(0,1),ylim=c(0,1), xlab='', ylab='',type='n') # zeichne Einheitsquadrat ein: rect(0,0,1,1) dateiname <- "quadrat_ohneB.pdf" dev.copy2pdf(file=dateiname) # zeichne B ein: polygon(r2*cos(t)+0.5,r2*sin(t)+0.6,col='lightblue') # (zeichne Randlinie fett nach) points(r2*cos(t)+0.5,r2*sin(t)+0.6, xlim=c(0,1),ylim=c(0,1),type='l',lwd=2) dateiname <- "quadrat_B.pdf" dev.copy2pdf(file=dateiname) # erzeuge 100 zufaellige Punkte: erfolge <- 0 for (i in 1:100) { x <- runif(1); y <- runif(1) if (drin(x,y)) { erfolge <- erfolge+1 points(x,y,col='darkblue',pch=20) # auch pch='.' ? } else points(x,y, col='gray',pch=20) } # Anteil Punkte in B: erfolge/100 title(paste('100 zufällige Punkte, Anteil in B :',erfolge/100)) dateiname <- "quadrat_B_100Pkte-3.pdf" dev.copy2pdf(file=dateiname) ## Verteilung des Anteils zeichnen pw <- 0.25557220603 plot(0, xlim=c(0,1),ylim=c(0,0.15), xlab='Wert von M', ylab='Wahrscheinlichkeit',type='n') for (k in 0:100) { lines(c(k/100,k/100),c(0,dbinom(k,size=100,prob=pw))) } dateiname <- "binomialgewichte100.pdf" dev.copy2pdf(file=dateiname) ## Verteilung des Anteils zeichnen pw <- 0.25557220603 plot(0, xlim=c(0,1),ylim=c(0,0.05), xlab='Wert von M', ylab='Wahrscheinlichkeit',type='n') for (k in 0:500) { lines(c(k/500,k/500),c(0,dbinom(k,size=500,prob=pw))) } dateiname <- "binomialgewichte500.pdf" dev.copy2pdf(file=dateiname) ## Illustration des Gesetzes der grossen Zahlen: N <- 10000 erf.prot <- numeric(N) for (i in 1:N) { x <- runif(1); y <- runif(1) erf.prot[i] <- ifelse(drin(x,y),1,0) } anteil.cum <- cumsum(erf.prot)*(1/(1:N)) ausw <- seq(from=100, to=10000, by=100) plot(ausw,anteil.cum[ausw],type='l', xlab='Anzahl Punkte',ylab='Beobachteter Anteil in B') abline(h=pw, lty=2) dateiname <- "illustrationGGZ.pdf" dev.copy2pdf(file=dateiname)