C 曲線

前に作ったはずなのに見つからなかった
ドラゴンカーブのプログラムとほんのちょっとの違いで，結果はかなり違うものになる

drawC <- function(a.x, a.y, b.x, b.y, n)
{
    x <- b.x-a.x
    y <- a.y-b.y
    c.x <- a.x+(x-y)/2
    c.y <- b.y-(x-y)/2
    if (n == 0) {
        lines(c(a.x, c.x, b.x), c(a.y, c.y, b.y), col="red")
    }
    else {
        drawC(a.x, a.y, c.x, c.y, n-1)
        drawC(c.x, c.y, b.x, b.y, n-1)
    }
}
old <- par(mar=c(0, 0, 0, 0))
plot(c(-6, 4), c(-8, 8), type="n", axes=FALSE, xlab="", ylab="", asp=1)
drawC(0, 4, 0, -4, 15)
par(old)

描画結果

n=15

n=12

n=10

πを求めるシミュレーション

Taglibro de H
モンテカルロ法
http://ito-hi.blog.so-net.ne.jp/2006-09-05

短くすればよいと言うものではないが，二次元配列を利用すると良い
配列の二乗，列和に colSums，列和が 1 以下になる平均値（！）を求め（sum(p)/n って，mean(p) のこと）4倍する

set.seed(124)
n <- 100000
mean(colSums(matrix(runif(2*n), 2)^2) < 1)*4

答え　3.14484

フラクタルの例として木

これも以前書いたもの

drawTree <- function(a.x, a.y, b.x, b.y, n, col=1)
{
    STEM.RATIO <- 0.25
    BRANCH.RATIO <- 0.6
    xx <- b.x-a.x
    yy <- a.y-b.y
    angle1 <- atan(yy/xx)+pi/6
    angle2 <- atan(yy/xx)-pi/6*1.5
    center.length <- sqrt(xx^2+yy^2)*(1-STEM.RATIO)
    branch.length <- BRANCH.RATIO*center.length
    sig <- ifelse(xx >= 0, 1, -1)
    c.x <- a.x+STEM.RATIO*xx
    c.y <- a.y-STEM.RATIO*yy
    d.x <- c.x+sig*(branch.length*cos(angle1))
    d.y <- c.y-sig*(branch.length*sin(angle1))
    e.x <- c.x+sig*(branch.length*cos(angle2))
    e.y <- c.y-sig*(branch.length*sin(angle2))
    segments(a.x, a.y, c.x, c.y, col=col)
    if (n <= 0) {
        segments(c.x, c.y, c(b.x, d.x, e.x), c(b.y, d.y, e.y), col=col)
    }
    else {
        drawTree(c.x, c.y, b.x, b.y, n-1, col=col)
        drawTree(c.x, c.y, d.x, d.y, n-1, col=col)
        drawTree(c.x, c.y, e.x, e.y, n-1, col=col)
    }
}
plot(c(00,500), c(0, 500), type="n", axes=FALSE, xlab="", ylab="", asp=1)
for (i in 1:10) {
    x <- sample(500, 1)
    y1 <- sample(500, 1)
    y2 <- y1-runif(1, min=100, max=400)
    drawTree(x, y2, x, y1, sample(4, 1)+2, col=i)
}

描画結果

コッホ曲線

これも，以前書いたもの

drawKoch <- function(a.x, a.y, b.x, b.y, n)
{
    c.x <- (2*a.x+b.x)/3
    c.y <- (2*a.y+b.y)/3
    d.x <- (a.x+2*b.x)/3
    d.y <- (a.y+2*b.y)/3
    x <- b.x-a.x
    y <- a.y-b.y
    d <- sqrt(x^2+y^2)/sqrt(3)
    if (x >= 0) {
        a1 <- atan(y/x)+pi/6
        e.x <- a.x+d*cos(a1)
        e.y <- a.y-d*sin(a1)
    }
    else {
        a2 <- atan(y/x)-pi/6
        e.x <- b.x+d*cos(a2)
        e.y <- b.y-d*sin(a2)
    }
    if (n <= 0) {
        lines(c(a.x, c.x, e.x, d.x, b.x), c(a.y, c.y, e.y, d.y, b.y), col="red")
    }
    else {
        drawKoch(a.x, a.y, c.x, c.y, n-1)
        drawKoch(c.x, c.y, e.x, e.y, n-1)
        drawKoch(e.x, e.y, d.x, d.y, n-1)
        drawKoch(d.x, d.y, b.x, b.y, n-1)
    }
}
plot(c(50, 450), c(70, 428), type="n", axes=FALSE, xlab="", ylab="", asp=1)
r <- 180
ox <- oy <- 250
p.x <- r*cos(pi*(3/6))+ox
p.y <- r*sin(pi*(3/6))+oy
q.x <- r*cos(pi*(7/6))+ox
q.y <- r*sin(pi*(7/6))+oy
r.x <- r*cos(pi*(11/6))+ox
r.y <- r*sin(pi*(11/6))+oy
n <- 7
drawKoch(p.x, p.y, q.x, q.y, n)
drawKoch(q.x, q.y, r.x, r.y, n)
drawKoch(r.x, r.y, p.x, p.y, n)
grid(40)
points(c(p.x, q.x, r.x), c(p.y, q.y, r.y))

描画結果

ドラゴン曲線

Taglibro de H
ドラゴン曲線
http://ito-hi.blog.so-net.ne.jp/2008-09-06

前に作っていたものだけど

drawDragon <- function(a.x, a.y, b.x, b.y, n)
{
    x <- b.x-a.x
    y <- a.y-b.y
    c.x <- a.x+(x+y)/2
    c.y <- b.y+(x+y)/2
    if (n <= 0) {
        lines(c(a.x, c.x, b.x), c(a.y, c.y, b.y), col="red")
    }
    else {
        drawDragon(a.x, a.y, c.x, c.y, n-1)
        drawDragon(b.x, b.y, c.x, c.y, n-1)
    }
}
plot(c(50, 450), c(70, 428), type="n", axes=FALSE, xlab="", ylab="", asp=1)
p.x <- 170
p.y <- 140
q.x <- 400
q.y <- 350
drawDragon(p.x, p.y, q.x, q.y, 10)
grid(40)

描画結果

スピログラフ（ついでに）

Taglibro de H
http://ito-hi.blog.so-net.ne.jp/2009-11-05

ついでにスピログラフのプログラムを while などを使わずに書いてみる。
ずいぶん短くなるし，描画もあっという間であるけれど。
while で描く方が，とろとろ描いていておもしろいのだけど（笑）。

spirograph <- function(
    gr,                # 外円の半径
    sr,                # 内円の半径
    pr,                # 内円内のペンの位置
    accuracy=100,        # 描画精度
    col = "black",        # 描画色
    lwd = 1)            # 描画線幅
{
    gcd <- function(a, b) {
        r <- a %% b
        return(ifelse(r == 0, b, gcd(b, r)))
    }

    gr <- trunc(gr)
    sr <- trunc(sr)
    pr <- trunc(pr)
    if (gr <= 0 | sr <= 0 | pr <= 0) {
        stop("Parameters must be postive.")
    }
    theta0 <- seq(0, 2*pi*sr/gcd(gr, sr), length=accuracy*sr/gcd(gr, sr))
    theta1 <- theta0*(gr/sr)
    x <- (gr-sr)*sin(theta0)-pr*sin(theta1)
    y <- (sr-gr)*cos(theta0)-pr*cos(theta1)
    old <- par(mar=c(0, 0, 0, 0))
    plot(x, y, type="n", axes=FALSE, xlab="", ylab="", asp=1)
    lines(x, y, lwd=lwd, col=col)
    par(old)
}

layout(matrix(1:4, 2))
spirograph(150, 108, 60)
spirograph(72, 552, 600, lwd = 3)
spirograph(11, 30, 45, col="red")
spirograph(15*29, 15*23, 50)
layout(1)

描画結果

ユークリッドの互除法

良くあるパターンであるが，引数を入れ替える必要はない

Taglibro de H
http://ito-hi.blog.so-net.ne.jp/2009-11-05

  ## 最大公約数
  gcd <- function(a, b) {
    if (a < b) {
      t <- a
      a <- b
      b <- t
    }
    r <- a %% b
    return(ifelse(r == 0, b, gcd(b, r)))
  }

a が b より小さい場合でも，一回余分に関数呼び出しがあるだけで，ちゃんと解が求まる。

  gcd2 <- function(a, b) { # これでいいのだ
    r <- a %% b
    return(ifelse(r == 0, b, gcd2(b, r)))
  }

> gcd(168, 48)
[1] 24
> gcd(48, 168)
[1] 24

> gcd2(168, 48)
[1] 24
> gcd2(48, 168)
[1] 24

matrix ではなく data.frame を使うべき場合

Taglibro de H

http://ito-hi.blog.so-net.ne.jp/2008-12-02?comment_success=2010-09-01T15:26:42&time=1283322402
R: order()で並べ替え　[統計] [編集]

こういうデータがあったとする。

m <- matrix(c(
    261, "コナラ","Quercus serrata",
     65, "ウリハダカエデ","Acer rufinerve",
    288, "アケビ","Akebia quinata",
    177, "モチツツジ","Rhododendron macrosepalum",
    288, "ミツバアケビ","Akebia trifoliata",
     51, "タラノキ","Aralia elata",
    261, "アラカシ","Quercus glauca",
    261, "アベマキ","Quercus variabilis",
    177, "コバノミツバツツジ","Rhododendron reticulatum",
    171, "エゴノキ","Styrax japonica"),
  ncol = 3, byrow = TRUE)

途中省略

ただし、日本語ではうまく並び替えできない。
> m[order(as.numeric(m[,1]), m[,2]),]
      [,1]  [,2]                 [,3]                      
 [1,] "51"  "タラノキ"           "Aralia elata"            
 [2,] "65"  "ウリハダカエデ"     "Acer rufinerve"          
 [3,] "171" "エゴノキ"           "Styrax japonica"         
 [4,] "177" "モチツツジ"         "Rhododendron macrosepalum"
 [5,] "177" "コバノミツバツツジ" "Rhododendron reticulatum"
 [6,] "261" "コナラ"             "Quercus serrata"         
 [7,] "261" "アベマキ"           "Quercus variabilis"      
 [8,] "261" "アラカシ"           "Quercus glauca"          
 [9,] "288" "アケビ"             "Akebia quinata"          
[10,] "288" "ミツバアケビ"       "Akebia trifoliata"  

matrix では 1 列目も character になってしまっている。
また，2 列目は character だけど，エンコーディングのバイトが辞書順になっていない。
解決法は，元のデータを data.frame として定義する。そうすれば，うまくいく。

> m <- read.csv(stdin(), header=FALSE)
0:     261, "コナラ","Quercus serrata"
1:      65, "ウリハダカエデ","Acer rufinerve"
2:     288, "アケビ","Akebia quinata"
3:     177, "モチツツジ","Rhododendron macrosepalum"
4:     288, "ミツバアケビ","Akebia trifoliata"
6:      51, "タラノキ","Aralia elata"
7:     261, "アラカシ","Quercus glauca"
8:     261, "アベマキ","Quercus variabilis"
9:     177, "コバノミツバツツジ","Rhododendron reticulatum"
10:     171, "エゴノキ","Styrax japonica"
11:
>

> m
    V1                  V2                        V3
1  261              コナラ           Quercus serrata
2   65      ウリハダカエデ            Acer rufinerve
3  288              アケビ            Akebia quinata
4  177          モチツツジ Rhododendron macrosepalum
5  288        ミツバアケビ         Akebia trifoliata
6   51            タラノキ              Aralia elata
7  261            アラカシ            Quercus glauca
8  261            アベマキ        Quercus variabilis
9  177  コバノミツバツツジ  Rhododendron reticulatum
10 171            エゴノキ           Styrax japonica

> m[order(m[,1], m[,2]),]
    V1                  V2                        V3
6   51            タラノキ              Aralia elata
2   65      ウリハダカエデ            Acer rufinerve
10 171            エゴノキ           Styrax japonica
9  177  コバノミツバツツジ  Rhododendron reticulatum
4  177          モチツツジ Rhododendron macrosepalum
8  261            アベマキ        Quercus variabilis
7  261            アラカシ            Quercus glauca
1  261              コナラ           Quercus serrata
3  288              アケビ            Akebia quinata
5  288        ミツバアケビ         Akebia trifoliata