COS论坛 | 统计之都

统计学论坛

1. 1楼

   

   cloud_wei
   

   资深会员
   注册于: 2007/04/21
   发帖数: 929

   用R来玩几个经典的分形图，很好玩。感兴趣的运行一下，跟帖发个图片上来，多谢支持，谢谢:)
   
   第一个：
   

   plot.tri <- function(n = 1000, col ="blue", ani=FALSE, cex=1.2){
     p <- runif(n); 
     X <- rbind(rep(0, n), rep(0, n))
     B <- cbind(c(0,0),c(0.25,0.433),c(0.5,0))
     if(ani) plot(0,0,xlim=c(0,1),ylim=c(0,0.85),type="n",xlab="",ylab="")
     for(i in 2:n){ 
           pp <- p[i]; 
           ind <- rank(c(c(1/3,2/3,1), pp), ties.method="min")[4]
           X[,i] <- 0.5*X[,i-1] + B[,ind]
           if(ani) points(X[1,i], X[2,i],pch = ".", cex = 1, col = col)
       }
     if(!ani) plot(X[1,], X[2,],pch = ".", cex = cex, col = col, xlab="", ylab="") 
   }
   
   ### example
   plot.tri(100000, ani=TRUE)

   
   
   第二个：
   

   plot.koch <- function(k,col="blue"){ 
     plot(0,0,xlim=c(0,1), ylim=c(-sqrt(3)/6,sqrt(3)/2), asp = 1,type="n",xlab="", ylab="")
     plotkoch <- function(x1,y1,x2,y2,n){
       if (n > 1){
         plotkoch(x1,y1,(2*x1+x2)/3,(2*y1+y2)/3,n-1); 
         plotkoch((2*x1+x2)/3,(2*y1+y2)/3,(x1+x2)/2-(y1-y2)*sqrt(3)/6,(y1+y2)/2-(x2-x1) *sqrt(3)/6,n-1);
         plotkoch((x1+x2)/2-(y1-y2)*sqrt(3)/6,(y1+y2)/2-(x2-x1)*sqrt(3)/6,(2*x2+x1)/3,(2 *y2+y1)/3,n-1); 
         plotkoch((2*x2+x1)/3,(2*y2+y1)/3,x2,y2,n-1)
       }    else { 
         x=c(x1,(2*x1+x2)/3,(x1+x2)/2-(y1-y2)*sqrt(3)/6,(2*x2+x1)/3,x2); 
         y=c(y1,(2*y1+y2)/3,(y1+y2)/2-(x2-x1)*sqrt(3)/6,(2*y2+y1)/3,y2); 
         lines(x,y,type="l",col=col) 
       }
     }
     plotkoch(0,0,1,0,k) 
     plotkoch(0.5,sqrt(3)/2,0,0,k) 
     plotkoch(1,0,0.5,sqrt(3)/2,k)
   }
   
   ## example
   for(i in 1:5){
     plot.koch(i,col=i)  
     Sys.sleep(1)
   }

   
   
   第三个：
   

   plot.leaf <- function(n=100000, col="green",cex=2){
     x <- c(.5, .5);
     plot(x[1], x[2], xlim=c(-3, 3), ylim =c(0, 10),type="n")
     p <- c( .85,  .92,  .99,  1.00);
     A <- rbind(c(.85, .04), c(-.04,.85), c(.20,-.26), c(.23,.22),
                c(-.15,.28), c(.26,.24),  c(0, 0), c(0, .16))
     B <- cbind(c(0, 1.6), c(0, 1.6), c(0,.44), c(0,0))
   
     for (i in 1:n){
        ran  <-  runif(1);
        ind <- rank(c(p, ran), ties.method="min")[5]
        x <- A[(2*ind-1):(2*ind),]%*%x + B[,ind]
       points(x[1],x[2], pch=".", cex=cex, col=col)
     }
   } 
   plot.leaf(cex=1.6)

   
   
   第四个：
   

   plot.tree <- function(x1,y1,x2,y2,n,xlim=c(-1,1), ylim=c(0,2), col="blue", add=FALSE){
     if(!add)
         plot(0,0,xlim=xlim, ylim=ylim, type="n",xlab="", ylab="",asp=1)
     tree <- function(x1,y1,x2,y2,n){
       flag <- 0;
       theta <- pi/6;
       if (x2<x1) flag <- 1;
       if (n>1){
         tree(x1,y1,(2*x1+x2)/3,(2*y1+y2)/3,n-1);
         tree((2*x1+x2)/3,(2*y1+y2)/3,(2*x2+x1)/3,(2*y2+y1)/3,n-1);
         tree((2*x2+x1)/3,(2*y2+y1)/3,x2,y2,n-1);
         tree((2*x1+x2)/3,(2*y1+y2)/3,(2*x1+x2)/3+sin(pi/2-atan((y2-y1)/(x2-x1))-theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3),(2*y1+y2)/3+cos(pi/2-atan((y2-y1)/(x2-x1))-theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3),n-1);
         tree((2*x2+x1)/3,(2*y2+y1)/3,(2*x2+x1)/3+sin(pi/2-atan((y2-y1)/(x2-x1))+theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3),(2*y2+y1)/3+cos(pi/2-atan((y2-y1)/(x2-x1))+theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3),n-1);
       } else {
         x <- c(x1,x2);
         y <- c(y1,y2);
         xx <- c((2*x1+x2)/3,(2*x1+x2)/3+sin(pi/2-atan((y2-y1)/(x2-x1))-theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3));
         yy <- c((2*y1+y2)/3,(2*y1+y2)/3+cos(pi/2-atan((y2-y1)/(x2-x1))-theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3));
         xxx <- c((2*x2+x1)/3,(2*x2+x1)/3+sin(pi/2-atan((y2-y1)/(x2-x1))+theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3));
         yyy <- c((2*y2+y1)/3,(2*y2+y1)/3+cos(pi/2-atan((y2-y1)/(x2-x1))+theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3));
         lines(x, y, type="l",col=col)
         lines(xx, yy, type="l",col=col)
         lines(xxx, yyy, type="l",col=col)
       }
     }
     tree(x1,y1,x2,y2,n)
   }
   
   ## example
   for(i in 1:5){
      plot.tree(0,0,0,1.5,i,col=i)
      Sys.sleep(1)
   }

   8 月 前回复 # 回复
2. 2楼

   

   lyxmoo
   

   版主
   注册于: 2009/02/01
   发帖数: 156

   我来传一个

   8 月 前回复 # 回复
3. 3楼

   

   lovelyday
   

   版主
   注册于: 2008/12/17
   发帖数: 127

   Re
   分形好玩

   8 月 前回复 # 回复
4. 4楼

   

   lyxmoo
   

   版主
   注册于: 2009/02/01
   发帖数: 156

   这个是example_2

   8 月 前回复 # 回复
5. 5楼

   

   微微
   

   初级会员
   注册于: 2007/11/23
   发帖数: 60

   穿上衣服膜拜一下~~
   怎么不发到R呢？

   8 月 前回复 # 回复
6. 6楼

   

   yanlinlin82
   

   无头衔
   注册于: 2008/11/27
   发帖数: 510

   经典漂亮的分形图形，可惜限于R的效率，画到数据很多时有些慢。

   8 月 前回复 # 回复
7. 7楼

   

   谢益辉
   

   个人主页
   注册于: 2006/05/19
   发帖数: 7,262

   是的，用R画会非常慢……
   
   这几个例子都收入animation包的demo()里了，将出现在animation 1.1-0中。

   8 月 前回复 # 回复
8. 8楼

   

   cloud_wei
   

   资深会员
   注册于: 2007/04/21
   发帖数: 929

   本身就是递归算法，慢是正常的，matlab也一样的~

   8 月 前回复 # 回复

RSS 订阅本帖

回复

您必须登录才能回复。