ダメ出し：Excel にすり寄る必要はないし，それ自体間違っている

「ggplot2 bar_plot with label」で，

「エクセル脳の皆さんに贈る」などといっているが，あなたこそ「ggplot脳」なのでは？

提示されたグラフは統計学的には不適切きわまりない。不適切なグラフに，どんなラベルを付けようが，不適切なグラフが適切なグラフになるわけはない。

適切なグラフは box-plot である。それが，「エクセル脳の皆さん」にも，「ggplot脳の皆さん」にも受け入れがたくても，仕方ない。

左のグラフを描くには，以下のようにする。
実に単純だ。なお，言うまでもないが，左の図の太い水平線は，平均値ではなく中央値だ。

boxplot(Sepal.Length~Species, data=iris)
median <- by(iris$Sepal.Length, iris$Species, median)
text(1:3+0.35, median, median, pos=4, xpd=TRUE)

ダメ出し：コンパイラなんかに頼らない

「TokyoR25」にて，library(compiler) の威力を見せつけようと？？しているようだが，compile したときとしないときの比較ではなくて，単に標準誤差はサンプルサイズの平方根に反比例するという統計学上の当たり前の話をなぞっているだけのような。

ここでも，もう何回となく取り上げてきて，うんざりということだが，プログラムを書く一人一人についてはうんざりでもなんでもなく，え～そんなことがあるの？ということなんだろうねえ。

元のプログラムは，シミュレーションごとに標本を生成している。それは，ダメだっっていってるでしょ。まとめて生成するの。必然的に for も使わなくて済むから。

> library(ggplot2)
> library(compiler)
>
> age <- data.frame(age = c(23, 25, 27, 27, 29, 26, 31, 31, 28, 29, 23, 26, 29,
+     30, 30, 27, 26))
>
> getMean <- cmpfun(function(times, s = 10) {
+     sample_mean <- numeric(times)
+     for (i in 1:times) {
+         res <- sample(age$age, size = s, replace = TRUE)
+         sample_mean[i] <- mean(res)
+     }
+     return(sample_mean)
+ })

だからね，13秒もかかっちゃうわけよ。

> system.time(res1 <- getMean(1e+06))
   ユーザ   システム       経過  
    13.384      0.143     13.650

あいかわらず，ggplot の図も「きたない」し。

> ggplot(data.frame(res1), aes(x = res1)) + geom_histogram(binwidth = 0.1, colour = "white") +
+     geom_vline(x = mean(res1), color = "red") + annotate(x = 26, y = 40000,
+     geom = "text", label = paste(sep = ":", "mean", round(mean(res1), 1))) +
+     annotate(x = 26, y = 30000, geom = "text", label = paste(sep = ":", "sd",
+         round(sd(res1), 3))) + xlim(24, 32)

まとめて標本を生成するとはこういうこと。

> getMean2 <- cmpfun(function(times, s = 10) {
+     res <- matrix(sample(age$age, size = s*times, replace = TRUE), times) # はい，ここですよ！
+     return(rowMeans(res)) # シミュレーション結果を返すのも一発で
+ })

で，こうなるわけさ。

> system.time(res2 <- getMean2(1e+06))
   ユーザ   システム       経過  
     0.231      0.057      0.287

> ggplot(data.frame(res2), aes(x = res2)) + geom_histogram(binwidth = 0.1, colour = "white") +
+     geom_vline(x = mean(res2), color = "red") + annotate(x = 26, y = 40000,
+     geom = "text", label = paste(sep = ":", "mean", round(mean(res2), 1))) +
+     annotate(x = 26, y = 30000, geom = "text", label = paste(sep = ":", "sd",
+         round(sd(res2), 3))) + xlim(24, 32)

定石なんだから，ああだこうだいわずに，いつでも採用すべき方策。

ダメ出し：方針が間違えている その2

そもそも，この種の問題を解くのに，ブルート・フォースを使うのは良くない。

問題をよく考えれば，理論的な解がちゃんと存在するし，そこまで行かなくても，コンピュータ・パワーを利用すると瞬殺の答えが得られる。

以下のようなプログラムを書くと，n=10 なら 0.028 秒！！，n=398 でも 0.065 秒！！！

n=399 以上は，numeric が扱える上限を超える。gmp ライブラリあたりを使えば良い。

dice.n <- function(n=10) {
    d <- c(1:6, 5:1)
    for (n in 3:n) {
        ld <- length(d)
        d2 <- numeric(ld+5)
        for (i in 1:6) {
            d2[i:(i+ld-1)] <- d2[i:(i+ld-1)]+d
        }
        d <- d2
    }
    return(d)
}

> system.time({
+ z <- dice.n(10)
+ barplot(z)
+ })
   ユーザ   システム       経過  
     0.028      0.002      0.031

> system.time({
+ z <- dice.n(398)
+ barplot(z)
+ })
   ユーザ   システム       経過  
     0.065      0.004      0.068

もう正規分布だね。

ダメ出し：方針が間違えている

「dice」にて，

> サイコロ３個ふって出目をカウント
> ６進法を使ってみた。

最後に

> サイコロ10個振ると？

とあるが，質問しているのか，挑戦しているのか，

少なくとも，掲載されているプログラムでは時間が掛かってやる気が起きない。

n <- 8
dice <- function(dec = 2, size = 2) (sum(mapply(as.integer, strsplit(dec2base(dec,
    6, size), "")) + 1))
a <- rep(0, 6*n)
names(a) <- 1:(6*n)

for (i in 0:(6^n-1)) {
    b <- dice(i, size = n)
    a[b] <- a[b] + 1
}
print(a[n:(6*n)])

8個のサイコロを振ったら，126秒ほどかかる。10個振ると4500秒くらいもかかると思う。

＝＝＝＝＝＝＝

別のアプローチとして，expand.grid というのを使う。

x expand.grid(x, x, x) としてできるものは，各行は３つのサイコロの出目である。行数は 6^3 行である。

10個のサイコロを振るなら，expand.grid(x,x,x,x,x,x,x,x,x,x) とすれば，6^10行，10列の行列ができる。rowSums で出目の和をとり，table 関数で集計し，barplot で図を描く。

実行時間は 50 秒ほど。

> system.time({
+ n <- 10
+ x <- 1:6
+ y <- eval(parse(text=paste("expand.grid(", paste(rep("x", n), collapse=","), ")", collapse="")))
+ barplot(table(rowSums(y)), las=1)
+ })
   ユーザ   システム       経過  
    50.398      7.074     57.070

 

ダメ出し：初心者の手本になるようなプログラムを描こう

「mandelbrot set」に，以下のようなプログラムが掲載されている。

mandelbrot <- function(step = 0.1, pch = 19, col = 30, n = 2) {
    a <- outer(seq(-2, 1, step), seq(-1, 1, step), function(x, y) complex(r = x,
        i = y))
    a <<- a[Mod(a) <= 2]  #絶対値２以下のものだけにする（あまり意味はないが）
    d <- c()
    dc <- c()
    for (c in a) {
        z <- 0
        for (i in 1:100) {
            z <- z^n + c
            if (Mod(z) >= 2) {
                d <- c(d, c)
                dc <- c(dc, i)  #何回目で発散したか
                break
            }
        }
    }
    # 回数で色分け
    plot(d, pch = pch, xlab = paste("step:", step, "n = ", n), xlim = c(-2,
        1), ylim = c(-1, 1), col = colors()[dc + col])
    dc <<- dc
}

c() とは何事かと思うだろうが，NULL と同じである。名は体を表すのだから，NULL って書こうよ。

> identical(NULL, c())
[1] TRUE

それに，d <- NULL にしておき，後で d <- c(d, foo) などとドンドン要素を追加していくのは，プログラム作法に反している。計算時間が指数関数的に増えていく。

必要なだけのメモリを確保して，要素を代入していくのが定石だ。未だにこの辺りのことを体得していないプログラマがいるのは嘆かわしい。

それに，a, b, c, d, dc のように，意味のない変数名を使うのも困ったことだ。

<<- を使っているのも問題。必要ならば，list にして return 文で返すべきだ。

同でもよいことだけど，注釈記号の後には少なくとも空白１個置きましょう。気持ち悪いから。

ということで，一番目の指摘点だけを書き換えるだけで，プログラムの速度はかなり速くなる。

step により指数関数的に増えるのだが，step=0.01 のときは，元のプログラムの３倍ほどの速さになる。

mandelbrot <- function(step = 0.1, pch = 0, col = 30, n = 2) {
    a <- outer(seq(-2, 1, step), seq(-1, 1, step), function(x, y) complex(real = x,
        imaginary = y))
    a <<- a[Mod(a) <= 2]  # 絶対値２以下のものだけにする（あまり意味はないが）
    len <- length(a)          # 追加
    d <- complex(len)     # 書き直し
    dc <- integer(len)       # 書き直し
    m <- 0                           # 追加
    for (c in a) {
        z <- 0
        for (i in 1:100) {
            z <- z^n + c
            if (Mod(z) >= 2) {
                m <- m+1      # 追加
                d[m] <- c     # 書き直し
                dc[m] <- i  # 何回目で発散したか  # 書き直し
                break
            }
        }
    }
    # 回数で色分け
      d <- d[1:m]   # 書き直し
    plot(d, pch = pch, xlab = paste("step:", step, "n = ", n), xlim = c(-2,
        1), ylim = c(-1, 1), col = colors()[dc + col], asp=1, cex=10*step)
    dc <<- dc[1:m]
}

ダメ出し：簡明直截なプログラミングをすれば，誤りもなくなる

「円を描く」にて

t <- seq(0, 2, pi/1000) * pi
x <- cos(t)
y <- sin(t)
plot(x, y, type = "l", asp = 1)

などと書いているが，

なんで，t <- seq(0, 2, pi/1000) * pi なのだろうか？意味不明である。
引数か何か誤解があるのでは？

それに，描かれた円は，閉じていない!!!

座標(1,0) のあたりを拡大すると，以下のようになっている！！

こんなことがないように，

t <- seq(0, 2*pi, length=1000)

でよいだろ

t <- seq(0, 2*pi, by=0.001)

は端点の処理という点で（t に 2*pi が含まれない），元のプログラムと同じで，今回の目的には合わない。

> t <- seq(0, 2, pi/1000) * pi
> sin(t[length(t)])
[1] -0.00611687              # ０に近くないとマズイ

> t <- seq(0, 2*pi, length=1000)
> sin(t[length(t)])
[1] -2.449294e-16            # ちゃんとやればこうなる。こうならないとだめ。

ダメ出し：ちゃんとしたグラフを描こう

「ggplot2:scale_shape_manual　ggplot２メモ：水準を形で分ける」において，以下のようなグラフを描いている。

水準が９個あり，色や線種を変えて折れ線グラフを描くと見にくいということで，記号も付けて描いたということだが，本当に見やすいのか？例えば，この折れ線の中から日本がどれなのか，探すのは簡単か？

簡単でない理由は，折れ線の凡例が本体と分離されているからだ。凡例が枠外にあるからではない。

以下のような，標準の matplot を使って描いた図と比べれば，差は一目瞭然。

見た目が派手（きれい？）なグラフを描くことに腐心するより，見やすいグラフを描くようにしたいものだ。パワポで見せびらかすにはよいかも知れないが，学会誌に ggplot2 で描いたような図は載らないだろう。

上の図を描いた元のプログラム

library(rdatamarket)
library(reshape2)
library(ggplot2)
getDM <- function(url) {
    data <- dmseries(url)
    data <- data.frame(year = index(data), as.data.frame(data))
    data <- melt(data, id.vars = "year")
    colnames(data) <- c("year", "country", "value")
    invisible(data)
}
url <- "http://datamarket.com/data/set/15n6/health-expenditure-public-of-gdp#display=line&ds=15n6|hnt"
res <- getDM(url)
smp <- subset(res, country %in% c("Sudan", "Mexico", "Kuwait", "Germany",
    "Finland", "France", "Canada", "Japan", "United.States"))
ggplot(smp, aes(x = year, y = value, group = country)) + geom_line(aes(color = country)) +
    geom_point(aes(color = country, shape = country)) + scale_shape_manual(values = 0:length(unique(smp$country)))

下の図を描いたプログラム
x <- matrix(smp$value, 15)
old <- par(mar=c(3, 3, 1, 6), mgp=c(1.8, 0.8, 0), las=1, cex.axis=0.8)
matplot(x, type="o", xaxt="n", xlab="year", ylab="value", lty=1, col=c(1:4, 6, 8), pch=1:9, cex=0.5)
axis(1, at=1:15, labels=1995:2009)
name <- c("Canada", "Finland", "France", "Germany", "Japan", "Kuwait", "Mexico", "Sudan", "United States")
axis(4, x[15,], name)
par(old)

ダメ出し：文字変換

「Jリーグ；J1チーム総年俸と順位の関係」に

全角数字列を半角数字列に変換するプログラム片がある。しかしまあ，これだけ簡単なものを複雑に書くものだなあというのが感想。無駄も目立つ。

# 全角数字→半角数字へ変換
nkf <- function(num.str = "３０００") {
    chars <- strsplit(num.str, "")[[1]]
    ret.num <- paste(apply(matrix(chars), 1, twobyte2onebyte), sep = "", collapse = "")
    return(ifelse(ret.num == "", 0, ret.num))
}

# 2byte-->1byte
twobyte2onebyte <- function(x = "３") {
    switch(x, ０ = return("0"), １ = return("1"), ２ = return("2"), ３ = return("3"),
        ４ = return("4"), ５ = return("5"), ６ = return("6"), ７ = return("7"),
        ８ = return("8"), ９ = return("9"), return(""))
}

chartr を使えば一発。この関数ならば，どのような文字変換でも，簡単に行う事ができる。カタカナ・ひらがな変換とか，暗号とかね。

nkf2 <- function(num.str = "３０００") {
    paste(chartr("０１２３４５６７８９", "0123456789", unlist(strsplit(num.str, ""))), collapse="")
}

全角文字ベクトルを作って，速さを検証してみた。

> zen <- chartr("0123456789", "０１２３４５６７８９", as.character(sample(0:9, 10000000, replace=TRUE)))
> system.time(nkf(zen))
   ユーザ   システム       経過  
    11.007      0.113     11.062
> system.time(nkf2(zen))
   ユーザ   システム       経過  
     7.475      0.040      7.483

まあ，大きな違いはない。