# Stochastik-Praktikum, WS 2010/2011 # Sitzung vom 8.11.2010 # Themen: # Grafische Datenanalyse, kurze Rekapitulation # (empirische) Verteilungsfunktion, Quantile # Schleifen und Funktionen in R # Simulation einer Irrfahrt, Rand-Auftreffproblem # Simulation von (Pseudo-)Zufallszahlen # ######################################################### # Nochmal kurz: Die Wein-Daten vom letzten Mal: wein <- read.table("whe.txt", header=TRUE, dec=",") wein attach(wein) Riesling stripchart(Riesling) stripchart(Riesling, method="jitter") boxplot(Riesling) # Empirische Verteilungsfunktion (der Rieslingerträge): plot.ecdf(Riesling, xlab="Ertrag", main="Empirische Vert. der Riesling-Erträge") # von Hand: x11() plot(sort(Riesling), (1:21)/21, type="S") q<- 0.5 quantile(Riesling, q) abline(h=q, col="darkgray",lty=2) abline(v=quantile(Riesling, q), col="darkgray",lty=2) # Ertrag als Funktion des Jahrs: plot(Jahr, Bacchus, main="Jahreserträge in hl/ha", type="l", ylab="") points(Jahr, Bacchus) points(Jahr, Riesling, col="red") points(Jahr, Riesling, col="red", type="l") legend("topright", legend=c("Bacchus", "Riesling"), col=c("black","red"), lty=1) # Es gibt offenbar bessere und schlechtere Weinjahre: # Bacchus gegen Riesling als Scatterplot plot(Bacchus, Riesling, main="Jahreserträge in hl/ha, 1989-2009") text(Bacchus, Riesling-1, labels=Jahr, cex=0.7, col="red") mean(Bacchus) abline(v=mean(Bacchus)) mean(Riesling) abline(h=mean(Riesling)) detach(wein) # ######################################################### # Schleifen und Funktionen in R: # 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^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 # 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="red") for (i in 1:schritte) { bneu <- f(b) lines(c(b,bneu,bneu), c(bneu,bneu,f(bneu)), col="red") b<-bneu } # Was Sie schon immer über Schleifen wissen wollten, # aber nie zu fragen wagten ... for(i in 1:10)print(i) # Schleife zaehlt von 1 bis 10 farben=c("rot", "gruen", "blau") for(i in farben){ # for schleife muss nicht numerisch sein cat("Farbe ist",i,"\n") # cat druckt, "\n" = Zeilenumbruch } wf=c(sin,cos,tan) # Vektor von winkelfunktionen for(w in wf){ # der Laufindex nimmt als Werte sin, cos, und tan an print(w(0.2)) # die jeweilige Winkelfunktion wird berechnet und ausgegeben } n <- 10 while(n>5){ # while schleife wird durchlaufen solange Bedingung wahr ist print(n<-n-1) # n bekommt neuen Wert n-1. Wert dieses Ausdrucks ist n-1. } n <- 10 while(n>1){ n <- n-1 if(n<5) next # next springt zum Ende der while Schleife und damit wieder zum Anfang print(n) } print(n) n <- 10 while(n>1){ n <- n-1 if(n<5) break # break springt aus der while Schleife raus print(n) } print(n) # break und next gehen auch in for schleife n <- 10 repeat{ # repeat auquivalent zu while(TRUE) n <- n-1 print(n) if(n>-5) next if(n<=1) # wird nur ausgefuehrt, falls n <= -5 und beendet dann die Schleife break } print(n) x <- c(-4,7, 8, 11) y <- c("rot","gruen", "gelb","blau") bed <- c(TRUE, FALSE, TRUE, FALSE) ifelse(bed, x, y) # waehlt x[i], falls bed[i]=TRUE, sonst y[i], vektorwertig ifelse(x>0, x, y) # Für sehr viele Simulationsprobleme benötigt man # sog. "Pseudozufallszahlen": # Beobachtungen oder Werte x[1], x[2], x[3], ... die keine # erkennbare Regelmäßigkeit haben und für die Zwecke der # Berechnung als unabhängig und zufällig generiert angenommen werden # (dürfen), d.h. wir tun beispielsweise so, als ob die x[1], x[2], ... # durch wiederholtes, unabhängiges Drehen eines Glücksrad gewonnen # worden wären -- obwohl sie aus Praktikabilitätsgründen durch einen # (deterministischen) Algorithmus (einen sog. (Pseudo-)Zufallsgenerator) # im Computer erzeugt worden sind. # # Rs eingebauter Zufallsgenerator: # 10 uniform auf [0,1] verteilte Pseudozufallszahlen runif(10) runif(10) # Reproduzierbarkeit erreicht man durch explizites Setzen des # random seed: set.seed(143) runif(10) set.seed(143) runif(10) werte2 <- runif(100) hist(werte2, prob=TRUE) plot(werte2, werte2[c(2:N,1)]) # Münzwürfe (und allgemeiner binomialverteilte ZVn) # erzeugt rbinom(1, prob=0.5) # "wirft eine faire Münze" rbinom(15, prob=0.3) # "wirft 10 verzerrte Münzen (Erfolgsw'keit=0.3) # erzeuge eine Binomial(20, 0.5)-verteilte (Pseudo-)Zufallsgröße: rbinom(1, size=20, prob=0.5) # ########################################################## # Das "Problem der Punkte" von Fra Luca Paccioli (1494), # (gelöst von Pascal und von Fermat im 17. Jh.): # Spieler A und B werfen eine faire Münze # A gewinnt, wenn insgesamt 4 Köpfe gefallen sind, bevor # 4 Zahlen gefallen sind, sonst gewinnt B. # # Die von L. Paccioli berichtete Streitfrage war: # Wenn nach 3 Runden 2x Kopf und 1x Zahl gefallen sind und # das Spiel nun unterbrochen werden muss, wie soll man den Einsatz # aufteilen? # Wir gehen die Frage zunächst empirisch ann: simuliere.Spiel <- function(koepfe, zahlen) { # Werfe solange faire (Pseudozufalls-)Münze, bis # kumulierte Anzahl beobachtete Köpfe # oder kumulierte Anzahl beobachtete Zahlen >= 4 while(koepfe < 4 && zahlen < 4) { if (runif(1)<0.5) koepfe <- koepfe+1 else zahlen <- zahlen+1 } # Gebe 1 zurück, wenn "Kopf-Spieler" gewinnt, sonst 0 if (koepfe==4) 1 else 0 } # Einmal spielen () simuliere.Spiel(0,0) # Zur Kontrolle: Das ursprüngliche Spiel ist (theoretisch) fair: versuche <- 1000 gewinne <- sum(replicate(versuche, simuliere.Spiel(0,0))) gewinne/versuche # Gegeben, dass nach 3 Runden 2x Kopf und 1x Zahl gefallen sind, # wie wahrscheinlich ist es, dass der "Kopf-Spieler" A gewinnt? versuche <- 1000 gewinne <- sum(replicate(versuche, simuliere.Spiel(2,1))) # geschätzte Gewinnw'keit von A: gewinne/versuche v <- c(1000, 5000, 10000, 20000, 50000) ews <- numeric(5) for (i in 1:5) { ews[i] <- sum(replicate(v[i], simuliere.Spiel(2,1)))/v[i] } plot(log(v), ews, xlab="log(Anz. Wdh.)", ylab="geschätzte Erfolgsw'keit für A") # (theoretische Antwort: 11/16) abline(11/16,0, lty=2) # ########################################################## # Simulation einer Irrfahrt, Rand-Auftreffproblem # Simuliere einen "Irrfahrtspfad"=Gewinnprozess im Münzwurfspiel: # In jeder Runde wird Münze mit Erfolgsw'keit p geworfen # bei Ergebnis "1" erhält Spieler A einen Euro von Spieler B, # bei Ergebnis "0" zahlt Spieler A einen Euro an Spieler B. n <- 20 # Anz. Schritte p <- 0.5 # Erfolgsw'keit der Münze k <- 0 # Startwert wuerfe <- (2*rbinom(n, size=1, prob=p)-1) plot(wuerfe) pfad <- k+cumsum(wuerfe) # pfad[i]=Gewinn von Spieler A nach i Runden plot(pfad, xlab="Schritte", main="Ein Irrfahrtspfad") points(pfad, type="l", add=TRUE) # füge zur besseren Sichtbarkeit Linien hinzu # Nehmen wir an, A und B spielen, bis einer 10 Euro verloren hat # (und der aktuelle Gewinn von A sei k) MW.Spiel <- function(k) { gewinn <- k while(gewinn > -10 && gewinn <10) gewinn <- gewinn+(2*rbinom(1,1,p)-1) if (gewinn==10) 1 else 0 } MW.Spiel(0) replicate(10, MW.Spiel(0)) sum(replicate(1000, MW.Spiel(0)))/1000 # Wie sieht die Gewinnwahrscheinlichkeit von A als Funktion des # Startwerts (="Vorsprung") aus? startwerte <- (-9):9 ewsA <- numeric(19) versuche <- 1000 for (i in 1:19) { ewsA[i]<-sum(replicate(versuche, MW.Spiel(startwerte[i])))/versuche cat("Startwert",startwerte[i],"geschätzte Gewinnw'keit für A:",ewsA[i],"\n") } plot(startwerte, ewsA, xlab="Startwert", ylab="geschätzte Gewinnw'keit für A") abline(0.5,1/20,lty=2) # Gerade 0.5+(1/20)*x einzeichnen # Analog kann man die Frage für p != 1/2 angehen # (dann sieht die exakte Antwort allerdings anders aus). # In ähnlicher Weise kann man z.B. die Frage studieren, wie # lange das Spiel dauert (und die Verteilung der Anz. Würfe bis zur # Entscheidung empirisch untersuchen). # ######################################################### # Simulation von (Pseudo-)Zufallszahlen # Für sehr viele Simulationsprobleme benötigt man # sog. "Pseudozufallszahlen": # Beobachtungen oder Werte x[1], x[2], x[3], ... die keine # erkennbare Regelmäßigkeit haben und für die Zwecke der # Berechnung als unabhängig und zufällig generiert angenommen werden # (dürfen), d.h. wir tun beispielsweise so, als ob die x[1], x[2], ... # durch wiederholtes, unabhängiges Drehen eines Glücksrad gewonnen # worden wären -- obwohl sie aus Praktikabilitätsgründen durch einen # (deterministischen) Algorithmus (einen sog. (Pseudo-)Zufallsgenerator) # im Computer erzeugt worden sind. # Ein einfaches Beispiel sind die # linearen Kongruenzengeneratoren: # Beginne mit einem (ganzzahligen) Startwert ("Zufallssame", "random seed") # naechster Wert ist # (a*x+c) mod M # z.B. M <- 2048; a <- 65; c <- 1 # "Selbstbau"-linearer Kongruenzengenerator: linKonGen <- function() { neu <- (a*x + c) %% M # folgendes setzt den Wert der (globalen) Variable x auf den von neu assign("x", neu, .GlobalEnv) # gebe "uniform" auf [0,1] verteilten Wert aus: neu/M } # Startwert: x <- 5 N <- 100 werte <- numeric(N) for (i in 1:N) werte[i]<-linKonGen() # Übrigens: Dasselbe leistet (wesentlich schneller) # werte <- replicate(N, linKonGen()) # Anschauen: plot(werte) lines(werte,add=TRUE) hist(werte, prob=TRUE) # Zeichne Wert gegen Nachfolger-Wert: plot(werte, werte[c(2:N,1)]) # Perioden können ein Problem sein (speziell bei kleinem M): M <- 16 a <- 5; c<-1 x <- 5 replicate(30, linKonGen()) # (in diesem Beispiel: Wert 4 = Wert 17, etc.)