#Problem 5a: chaos game - Sierpinski Dreieck #wir starten in einem festen Punkt (x0,y0), ziehen zufällig eine der Funktionen f1,f2,f3 und wenn die #gezogene Funktion auf (x0,y0) an. Mit dem so erhaltenen Punkt verfahren wir analog #Frage: wie viel Prozent der Punkte liegen im unteren linken Dreieck? #Schritt 0: Definition des IFSP f1 <- function(x){r <- x/2;return(r)} #die drei Funktionen f2 <- function(x){r <- x/2 + c(1/2,0);return(r)} f3 <- function(x){r <- x/2 + c(1/4,sqrt(3)/4);return(r)} IFS <- list(f1=f1,f2=f2,f3=f3) #Schreiben Funktionen in Liste p <- c(1/3,1/3,1/3) #Auswahlwahrscheinlichkeiten fnum <- 1:length(IFS) #Nummern der Funktionen #Schritt 1: generieren einen Pfad des chaos games n <- 10000 #Länge des Pfads z <- c(0,0.5) #Startpunkt A <- data.frame(x=rep(0,n),y=rep(0,n)) index <- sample(fnum,size=n,replace=TRUE,prob=p) #ziehen gleich n Funktionen mit Zurücklegen for(i in 1:n){ z <- IFS[[index[i]]](z) #iteratives Anwenden der gezogenen Funktionen A[i,] <- z #Abspeichern der Koordinaten (analog random walk) } plot(A,cex=0.25) #plot der Punkte #Modifikation Schritt 1: flottere Version des Codes #blockweise Berechnung mehrerer runs des chaos games n <- 50000 p <- c(1/6,1/3,1/2) block <- 2000 beg <- 100 R<-trunc(n/block) indices <- c(0,block*(1:(R-1)),n) AA <- vector("list",length=R) for(r in 1:R){ z <- c(runif(1),runif(1)) nr <- indices[r+1]-indices[r] A <- data.frame(x=rep(0,nr),y=rep(0,nr)) index <- sample(fnum,size=nr,replace=TRUE,prob=p) for(i in 1:nr){ z <- IFS[[index[i]]](z) A[i,] <- z } AA[[r]] <- A[beg:nr,] } A <- do.call("rbind",AA) plot(A,cex=0.25) #nicer plot including histogram #plot points and two-dim histogram library(ggplot2) library(gridExtra) p <- ggplot(A, aes(x=x, y=y)) p <- p + geom_point(size=0.1,colour="black") p <- p + labs(title="Sierpinski Triangle: Orbit Chaos Game") p <- p + scale_x_continuous(name="",breaks=seq(0,1,length=4),labels=c("0","1/3","2/3","1")) p <- p + scale_y_continuous(name="",breaks=seq(0,1,length=4),labels=c("0","1/3","2/3","1")) p <- p + theme_bw() p01a<-p farbe<-rainbow(100,start=.40,end=.17) p <- ggplot(A, aes(x=x, y=y)) p <- p + stat_binhex(binwidth = c(0.01,0.01),colour=NA) p <- p + theme_bw() p <- p + scale_fill_gradientn(colours=farbe,name="count") p <- p + labs(title="Sierpinski Triangle: Histogram") p <- p + scale_x_continuous(name="",breaks=seq(0,1,length=4),labels=c("0","1/3","2/3","1")) p <- p + scale_y_continuous(name="",breaks=seq(0,1,length=4),labels=c("0","1/3","2/3","1")) p <- p + theme_bw() p01b<-p grid.arrange(p01a,p01b,ncol=2,widths=c(0.45, 0.55)) #-------------------------------------------------------------------------------- #-------------------------------------------------------------------------------- # star library(ggplot2) library(gridExtra) #define functions f1 <-function(x){x/3} f2 <-function(x){x/3+c(2/3,0)} f3 <-function(x){x/3+c(1/3,1/3)} f4 <-function(x){x/3+c(0,2/3)} f5 <-function(x){x/3+c(2/3,2/3)} IFS<-list(f1=f1,f2=f2,f3=f3,f4=f4,f5=f5) fnum <- 1:length(IFS) #choose starting point and run chaos game: #start chaos game: n <- 100000 block <- 2000 beg <- 100 R<-trunc(n/block) indices <- c(0,block*(1:(R-1)),n) AA <- vector("list",length=R) for(r in 1:R){ z <- c(runif(1),runif(1)) nr <- indices[r+1]-indices[r] A <- data.frame(x=rep(0,nr),y=rep(0,nr)) index <- sample(fnum,size=nr,replace=TRUE) for(i in 1:nr){ z <- IFS[[index[i]]](z) A[i,] <- z } AA[[r]] <- A[beg:nr,] } A <- do.call("rbind",AA) #plot points and two-dim histogram p <- ggplot(A, aes(x=x, y=y)) p <- p + geom_point(size=0.1,colour="black") p <- p + labs(title="Star: Orbit Chaos Game") p <- p + scale_x_continuous(name="",breaks=seq(0,1,length=4),labels=c("0","1/3","2/3","1")) p <- p + scale_y_continuous(name="",breaks=seq(0,1,length=4),labels=c("0","1/3","2/3","1")) p <- p + theme_bw() p01a<-p farbe<-rainbow(100,start=.40,end=.17) p <- ggplot(A, aes(x=x, y=y)) p <- p + stat_binhex(binwidth = c(0.01,0.01),colour=NA) p <- p + theme_bw() p <- p + scale_fill_gradientn(colours=farbe,name="count") p <- p + labs(title="Star: Histogram") p <- p + scale_x_continuous(name="",breaks=seq(0,1,length=4),labels=c("0","1/3","2/3","1")) p <- p + scale_y_continuous(name="",breaks=seq(0,1,length=4),labels=c("0","1/3","2/3","1")) p <- p + theme_bw() p01b<-p grid.arrange(p01a,p01b,ncol=2,widths=c(0.45, 0.55))