# Franziska Wandtner, 2012, wandtner@math.uni-frankfurt.de ########################################################################## # Wir geben uns eine Fläche im Einheitsquadrat vor und fragen nach deren Flächeninhalt. # Diesen wollen wir mit Hilfe einer Monte-Carlo-Methode annähern: # Es werden 10 rein zufällige Punkte in das Einheitsquadrat "geworfen" und jeweils bestimmt, ob sie in die markierten Fläche gefallen sind oder nicht. # Der relative Anteil der Punkte, die in die markierte Fläche gefallen sind, ist ein Schätzwert für den Flächeninhalt. # Im zweiten Teil soll untersucht werden, wie gut dieser Schätzwert ist, wenn man n rein zufällige Punkte wählt. # Hierzu wird das Verfahren mit n Punkten k mal wiederholt und ein Balkendiagramm der Schätzwerte erstellt. n und k können von Ihnen oben im Code festgelegt werden. # Beachte: Für k=1000 benötigt man schon eine spürbare Rechenzeit. Falls es zu lange dauert, kann man unten die Anzahl der Wiederholungen verringern. # Die Schätzwerte werden schließlich mit dem tatsächlichen Flächeninhalt verglichen. ########################################################################## rm(list=ls()) ############################ # Anzahl der Punkte und der Wiederholungen bei Teilen 2 und 3 können HIER festgelegt werden. Ansonsten brauchen Sie nichts ändern. n<-100 k<-5000 ########################################################################### #Graphikparameter und Fenstergröße für Graphik festlegen # Für Mac: quartz(width=6, height=6, pointsize = 12, family = "Helvetica", antialias = TRUE) # Für Windows: # windows(width=6, height=6) ########################## Unterprogramme ################################## # Graphik x<-c(0.1, 0.5, 0.6, 0.8) y<-c(0.9, 0.1, 0.6, 0.7) graph<-function() { par(mar=c(0,0,1,0),oma=c(0,0,1,0)) plot(c(0,1),c(0,1),t="n",axes=F,xlab="",ylab="",main=main) polygon(c(0,1,1,0),c(0,0,1,1)) # Figur polygon(x,y,col="turquoise") } ######################### Vorbereitungen ################################### # p ausrechnen p<-(x[1]*y[2]-x[2]*y[1])+(x[2]*y[3]-x[3]*y[2])+(x[3]*y[4]-x[4]*y[3])+(x[4]*y[1]-x[1]*y[4]) p<-abs(p/2) # Geraden zwischen Punkten f1<-function(x){1.1-2*x} f2<-function(x){-2.4+5*x} f3<-function(x){0.3+0.5*x} f4<-function(x){(13/14)-(2/7)*x} ############################### Hauptprogramm ############################### # Teil 1: Hier sollte man nichts ändern. # EINZELNE PUNKTE: # Plot vorbereiten main<-"Wie groß ist der Anteil der blauen Fläche?" graph() # Zufällige Punkte einzeln (10 Stück) drinnen<-0 readline("Einzelne Punkte") xr<-runif(10) yr<-runif(10) for(i in 1:10) { points(xr[i],yr[i],col="gray2",cex=1,pch=19) # Analyse: Ist Punkt drinnen? if(xr[i]==0.1&yr[i]==0.9){drinnen<-drinnen+1} if(xr[i]>0.1&xr[i]<0.5&yr[i]>=f1(xr[i])&yr[i]<=f4(xr[i])){drinnen<-drinnen+1} if(xr[i]==0.5&yr[i]==0.1){drinnen<-drinnen+1} if(xr[i]>0.5&xr[i]<0.6&yr[i]>=f2(xr[i])&yr[i]<=f4(xr[i])){drinnen<-drinnen+1} if(xr[i]==0.6&yr[i]==0.6){drinnen<-drinnen+1} if(xr[i]>0.6&xr[i]<0.8&yr[i]>=f3(xr[i])&yr[i]<=f4(xr[i])){drinnen<-drinnen+1} if(xr[i]==0.8&yr[i]==0.7){drinnen<-drinnen+1} if(i<=9){readline("Nächster Punkt")} } readline("Prozent in Menge") main<-"Schätzung des Anteils" graph() points(xr,yr,col="azure4",cex=1,pch=19) text(0.7,0.85,paste("Anteil Treffer:",drinnen*10,"%")) readline("Und jetzt genauer!") text(0.45,0.5,"Und jetzt genauer!",col="red",cex=1.3,font=2) readline("Mehr Punkte!") text(0.45,0.45,"Mehr Punkte!",col="red",cex=1.3,font=2) ############################################################################# # Teil 2 readline("Simuliere mehr Punkte") # VIELE PUNKTE # Plot vorbereiten main<-"Wie groß ist der Anteil der blauen Fläche?" graph() # n rein zufällige Punkte xr<-runif(n) yr<-runif(n) points(xr,yr,col="gray2",cex=0.3,pch=19) #Beschriftung text(0.8,0.1,paste("n =",n)) # Analyse: Ist Punkt drinnen? drinnen<-0 for(i in 1:n) { if(xr[i]==0.1&yr[i]==0.9){drinnen<-drinnen+1} if(xr[i]>0.1&xr[i]<0.5&yr[i]>=f1(xr[i])&yr[i]<=f4(xr[i])){drinnen<-drinnen+1} if(xr[i]==0.5&yr[i]==0.1){drinnen<-drinnen+1} if(xr[i]>0.5&xr[i]<0.6&yr[i]>=f2(xr[i])&yr[i]<=f4(xr[i])){drinnen<-drinnen+1} if(xr[i]==0.6&yr[i]==0.6){drinnen<-drinnen+1} if(xr[i]>0.6&xr[i]<0.8&yr[i]>=f3(xr[i])&yr[i]<=f4(xr[i])){drinnen<-drinnen+1} if(xr[i]==0.8&yr[i]==0.7){drinnen<-drinnen+1} } drinnen #Schätzer für relativen Anteil pdach1<-drinnen/n readline("Anteil Treffer") main<- "Schätzung des Anteils" graph() text(0.8,0.1,paste("n =",n)) points(xr,yr,col="azure4",cex=0.3,pch=19) text(0.7,0.85,paste("Anteil Treffer:",round(pdach1*100,2),"%")) readline("Wie typisch war dieser Schätzwert?") text(0.5, 0.5,"Wie typisch war dieser Schätzwert?",col="red",cex=1.3,font=2) readline(paste("Wiederhole", k, "Mal!")) text(0.5, 0.45,paste("Wiederhole", k, "Mal!"),col="red",cex=1.3,font=2) ################## readline(paste(k, "Schätzungen - dauert eine Weile!!")) #VIELE SCHÄTZUNGEN: dauert eine Weile!! # k Schätzungen mit jeweils n Punkten (Festlegungen siehe oben) drinnen<-rep(0,k) for(j in 1:k) { xr<-runif(n) yr<-runif(n) # Analyse: Ist Punkt drinnen? for(i in 1:n) { if(xr[i]==0.1&yr[i]==0.9){drinnen[j]<-drinnen[j]+1} if(xr[i]>0.1&xr[i]<0.5&yr[i]>=f1(xr[i])&yr[i]<=f4(xr[i])){drinnen[j]<-drinnen[j]+1} if(xr[i]==0.5&yr[i]==0.1){drinnen[j]<-drinnen[j]+1} if(xr[i]>0.5&xr[i]<0.6&yr[i]>=f2(xr[i])&yr[i]<=f4(xr[i])){drinnen[j]<-drinnen[j]+1} if(xr[i]==0.6&yr[i]==0.6){drinnen[j]<-drinnen[j]+1} if(xr[i]>0.6&xr[i]<0.8&yr[i]>=f3(xr[i])&yr[i]<=f4(xr[i])){drinnen[j]<-drinnen[j]+1} if(xr[i]==0.8&yr[i]==0.7){drinnen[j]<-drinnen[j]+1} } } pdach<-drinnen/n breaks<-seq(max(min(pdach)-0.1,0),max(pdach)+0.1,by=1/(n*10)) quartz(width=8, height=6, pointsize = 12, family = "Helvetica", antialias = TRUE) hist<-hist(pdach,breaks=breaks,xlab=paste("Schätzwert des relativen Anteils der Fläche (",n," Punkte)"),freq=T,main=paste("Verteilung der Schätzwerte (aus",k, "Wiederholungen)"),col="black") # Balkendiagramm zeichnen hist(pdach,breaks=breaks,xlab=paste("Schätzwert des relativen Anteils der Fläche (",n," Punkte)"),freq=T,main=paste("Verteilung der Schätzwerte (aus",k, "Wiederholungen)"),col="black") readline("wahrer Anteil") points(0.195,0,col="red",cex=1.3,pch=19) text(max(min(pdach)-0.1,0),max(hist$counts),"Tatsächlicher Anteil der blauen Fläche: 0.195",col="red",cex=0.9,adj=c(0,0)) readline("erste Schätzung") arrows(max(min(pdach)-0.1,0),0.5*max(hist$counts),pdach1,0,col="blue") text(max(min(pdach)-0.1,0),0.6*max(hist$counts), paste("Unser erster Schätzwert war:",pdach1),cex=0.9,col="blue",adj=c(0,0)) ############################################################################# readline("Ende") graphics.off()