SSブログ

R: Mandelbrot [統計]

WikiNews: Mathematician Benoît Mandelbrot dies aged 85.

Andrew Gelmanのブログ記事: Mandelbrot and Akaike: from taxonomy to smooth runways (pioneering work in fractals and self-similarity).

ということで、マンデルブロー集合をかいてみる。
Mandelbrot.png
Mandelbrot2.png
コード
library(snow)

px <- 900
py <- 600
min.x <- -2
max.x <- 1
min.y <- -1
max.y <- 1
g <- expand.grid(x = seq(min.x, max.x, length = px),
                 y = seq(min.y, max.y, length = py))

man <- function(p, lim = 255, th = 2.0) {
	x <- p[, 1]
	y <- p[, 2]
	u <- 0
	z <- 0
	while (u < lim) {
		z <- z^2 + (x + y * 1i)
		if (abs(z) > th) break
		u <- u + 1
	}
	return(u)
}

##
cl <- makeCluster(2, type = "MPI")
clusterExport(cl, list("man", "g"))
z <- parSapply(cl, 1:nrow(g), function(i) man(g[i, ]))
stopCluster(cl)

matz <- matrix(z, nrow = px)
image(seq(min.x, max.x, length = px),
      seq(min.y, max.y, length = py), matz,
      xlab = "x", ylab = "y")


library(snow)

px <- 512
py <- 512
min.x <- -0.58
max.x <- -0.55
min.y <- -0.57
max.y <- -0.54
g <- expand.grid(x = seq(min.x, max.x, length = px),
                 y = seq(min.y, max.y, length = py))

man <- function(p, lim = 255, th = 2.0) {
	x <- p[, 1]
	y <- p[, 2]
	u <- 0
	z <- 0
	while (u < lim) {
		z <- z^2 + (x + y * 1i)
		if (abs(z) > th) break
		u <- u + 1
	}
	return(u)
}

##
cl <- makeCluster(2, type = "MPI")
clusterExport(cl, list("man", "g"))
z <- parSapply(cl, 1:nrow(g), function(i) man(g[i, ]))
stopCluster(cl)

matz <- matrix(z, nrow = px)
png(width = 900, height = 900)
image(seq(min.x, max.x, length = px),
      seq(min.y, max.y, length = py), matz,
      col = heat.colors(256),
      xlab = "x", ylab = "y")
dev.off()


並列化してみたけど、インタプリタでは やはりそれなりに時間がかかる。まあ、8ビット機の昔はアセンブラでも1晩以上かかっていたのだが。
タグ:R
nice!(0)  コメント(0)  トラックバック(1) 
共通テーマ:学問

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 1