そんなものじゃないだろう

> 　記者会見した長沢寛道・同研究科長は「高い値の放射性セシウムが検出された水田で作付けし、データを正確に把握することは、福島県の農業復興に必要だ」と話した。（東京大大学院農学生命科学研究科 長沢寛道研究科長）

作付けして，データを取ってもらえれば科学的知見を得るのに役立つという考えなのだろう。

しかし，農業に携わる人は，金銭的保証をしてくれればよいなどとは思っていないはず。食べてもらえない米を作るなどと言うのは耐えられない苦痛であるはず。科学者は，その苦痛を思いやれないのだろうなあと，やりきれない思いがする。

ダメ出し：アトキンの篩

2011-05-11　R でエラトステネスの篩　の後半に，アトキンの篩のプログラムが掲載されている。

ご本人も「多分条件判定減らすとかしてRに合った実装しないとダメなんだと思う」というとおり，ベクトル化をはかる。ちなみに，元のプログラムでは for を避けるためか while で書いているがこれはほとんど意味がない。

limit = 10000000 で 23 秒かかっていたものが，書き直すと 0.4 秒ほどで計算終了ということになった。

プログラムは以下の通り。

atkin <- function(limit = 1e+06, return = FALSE) {
    sqrt.limit <- sqrt(limit)
    isprime <- logical(limit)
    for (z in c(1, 5)) {
        for (y in seq.int(z, sqrt.limit, by = 6)) {
            y2 <- y * y
            max.x <- floor(sqrt((limit - y2) / 4))
            if (1 <= max.x) {
                n <- 4 * (1:max.x) ^ 2 + y2
                isprime[n] <- !isprime[n]
            }
            x <- y + 1
            max.x <- floor(sqrt((limit + y2) / 3))
            if (x <= max.x) {
                n <- 3 * seq.int(x, max.x, by = 2) ^ 2 - y2
                isprime[n] <- !isprime[n]
            }
        }
    }
    for (z in c(2, 4)) {
        for (y in seq.int(z, sqrt.limit, by = 6)) {
            y2 <- y * y
            max.x <- floor(sqrt((limit - y2) / 3))
            if (1 <= max.x) {
                n <- 3 * seq.int(1, max.x, by = 2) ^ 2 + y2
                isprime[n] <- !isprime[n]
            }
            x <- y + 1
            max.x <- floor(sqrt((limit + y2) / 3))
            if (x <= max.x) {
                n <- 3 * seq.int(x, max.x, by = 2) ^ 2 - y2
                isprime[n] <- !isprime[n]
            }
        }
    }
    for (y in seq.int(3, sqrt.limit, by = 6)) {
        y2 <- y * y
        for (x in 1:2) {
            max.x <- floor(sqrt((limit - y2) / 4))
            if (x <= max.x) {
                n <- 4 * seq.int(x, max.x, by = 3) ^ 2 + y2
                isprime[n] <- !isprime[n]
            }
        }
    }
    for (n in seq.int(5, sqrt.limit, by = 2)) {
        if (isprime[n]) {
            isprime[1:(limit %/% (n * n)) * n * n] <- FALSE
        }
    }
    isprime[2] <- isprime[3] <- TRUE
    if (return) {
        which(isprime)
    }
}

ダメ出し：R で円周率計算（2）

R で円周率計算

mean(runif(points.num)^2+runif(points.num)^2 < 1)*4

だけでよかった。速いし。

乱数をまとめて発生させるという技が裏目に出ることもある

シミュレーションなどで，使用する乱数をまとめて発生させておいて後で使うというようにすれば効率的であるといわれるが，メモリーギリギリでやるような場合には，それが裏目に出るという実例。おまけに，この場合は for ではなく apply を使ったのに...ということなのだ。

> n <- 1000
> gc()
         used (Mb) gc trigger  (Mb)  max used   (Mb)
Ncells 389200 20.8    1307810  69.9   2554318  136.5
Vcells 815860  6.3  101994529 778.2 151544141 1156.2
> set.seed(123)
> system.time({
+     x <- matrix(rnorm(50*n), 50)
+     m <- apply(x, 2, median)
+ })
   ユーザ   システム       経過 
     0.068      0.001      0.069         小規模な場合は，確かにまとめて取る方が効率的
> set.seed(123)
> system.time({
+     m <- numeric(n)
+     for (i in seq_len(n)) {
+         m[i] <- median(rnorm(50))
+     }
+ })
   ユーザ   システム       経過 
     0.086      0.001      0.093

> n <- 1000000
> gc()
           used  (Mb) gc trigger   (Mb)  max used  (Mb)
Ncells   385959  20.7    1259264   67.3   1835812  98.1
Vcells 51544018 393.3  149518151 1140.8 105639609 806.0
> system.time({
+     x <- matrix(rnorm(50*n), 50)
+     m <- apply(x, 2, median)
+ })
   ユーザ   システム       経過 
    73.343      0.863     73.887         大規模な場合は，逆に遅くなってしまう
> system.time({
+     m <- numeric(n)
+     for (i in seq_len(n)) {
+         m[i] <- median(rnorm(50))
+     }
+ })
   ユーザ   システム       経過 
    66.355      0.629     66.733