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

「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

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

おっぱい関数

「おっぱい関数」を材料として3次元表示についてまとめる。

「おっぱい関数」は元々誰が作ったか不詳だが，検索すると
Julia + Excel　http://www.hirax.net/diaryweb/2012/03/01.html
Maxima　https://plus.google.com/u/0/117731699005945352808/posts/Yswj5tswtNd
などが引っかかる。変形したものもある。
Python　http://kemeconoajito.blog88.fc2.com/blog-entry-175.html
R によるものは，http://siguniang.wordpress.com/2011/11/09/3d-surface-plot-with-latticergl/
にある。

library(lattice)
 
g <- function(x, y) {
 1/8 * (6 * exp(-((2/3 * abs(x) - 1)^2 + (2/3 * y)^2) - 1/3 * (2/3 * y + 1/2)^3) +
        2/3 * exp(-2.818^11 * ((abs(2/3 * x) - 1)^2 + (2/3 * y)^2)^2) +
        2/3 * y - (2/3 * x)^4)
}
m <- 50
x <- seq(-3, 3, length.out=m)
y <- seq(-3, 3, length.out=m)
grid   <- expand.grid(x=x, y=y)
grid$z <- g(grid$x, grid$y)
wireframe(z ~ x * y, grid, shade=TRUE)
library(rgl)
grid$col <- ifelse(sqrt((abs(grid$x) - 3/2) ^ 2 + grid$y ^2 ) < 0.45, 'hotpink', 'lightpink')# nipple
open3d()
rgl.surface(x, y, grid$z, color=grid$col)

この関数部分を小生が書き直した

n <- 100
x0 <- y0 <- seq(-3, 3, length=n+1)
xy <- expand.grid(x=x0, y=y0)
x <-  xy[,1]
y <-  xy[,2]
t0 <- (2 * abs(x) / 3 - 1) ^ 2
t1 <- 2 * y / 3
t2 <- t1 ^ 2
t4 <- -t0 - t2 - (t1 + 0.5) ^ 3 / 3
z <- 0.125 * (6 * exp(t4) +
     2 * exp(-exp(11) * (t0 + t2) ^ 2) / 3 +
     2.5 * t1 -16 * x ^ 4 / 81)
library(rgl)
plot3d(x, y, z, col = "burlywood1", aspect = c(1, 1, 0.6), axes=FALSE, xlab="", ylab="", zlab="")

# 以下の，ふたつの方法もやってみる
open3d()
rgl.surface(x0, y0, z, color="burlywood1")
library(lattice)
wireframe(z ~ x * y, xy, shade=TRUE, aspect=c(1, 0.5), xlab="", ylab="", zlab="")

# 以下で出力される csv ファイルを Excel で 3d 表示
write(matrix(z, n+1), file="oppai2.csv", ncolumns=n+1, sep=",")

変形したものも書いておく。
http://kemeconoajito.blog88.fc2.com/blog-entry-175.html

n <- 100
xy <- expand.grid(seq(-3, 3, length=n+1), seq(-3, 3, length=n+1))
x <-  xy[,1]
y <-  xy[,2]
t0 <- (2 * abs(x) / 3 - 1) ^ 2
t1 <- 2 * y / 3
t2 <- t1 ^ 2
t3 <- (t0 + t2) ^ 2
t4 <- -t0 - t2 -  (t1 + 0.5) ^ 3 / 3
z <- 0.125 * (7 * exp(t4) + 1.5 * exp(-exp(10) * t3) +
     exp(-exp(5) * t3) + 2.5 * t1 - 16 * x ^ 4 / 81)
library(rgl)
plot3d(x, y, z, col = "burlywood1", aspect = c(1, 1, 0.6), axes=FALSE, xlab="", ylab="", zlab="")

ついでに見つかった，2次元の「おっぱい関数」も R で書いたものを載せておく。

# https://twitter.com/math_neko/status/227672102487068672/photo/1
y <- seq(.Machine$double.eps, 1.2, length=300)
x <- 3*y*log(y)-1/36*exp(-(36*y-36/exp(1))^4)
plot(1, 1, type="n", xlim=c(-1.2, 1), ylim=c(-0.1, 1))
lines(x, y)

# http://goo.gl/rtZS7
n <- 5000
curve(sqrt(1-(x-1)^2), n=n, xlim=c(-2.5, 2.5), ylim=c(-0.1,  1.25), asp=1)
curve(sqrt(1-(x+1)^2), n=n, add=TRUE)
curve(sqrt(0.01-(x-1)^2)+1, n=n, add=TRUE)
curve(sqrt(0.01-(x+1)^2)+1, n=n, add=TRUE)

# http://i.imgur.com/IHAhR.jpg
t <- seq(-1.4, 1.5, length=200)
y <- 2*t^3-0.5*t-1
x <- 0.9*(exp(-t^4+t^2+t)+0.2*(exp(-100*((y-0.2)^2)))-0.5*exp(-(y-3.5)^2)-0.15*exp(-(y-5)^2))
plot(x, y, type="l", asp=1)

ダメ出し：ちゃんと関数があるんだから，モ～

久しぶりの大物

「多次元度数分布を取る」より

プログラム自体は，やっと <- の前後に空白をおくなどで少しは読みやすくなったけど，まだ，){ や if( や for( があったり，カンマの後に空白を置かないなど相変わらず，読みにくい。

3次元クロス集計をする自前の関数を書いたというところだけど，リストや for 文を惜しげもなく使っているので，処理時間がべらぼうにかかる。

> set.seed(123456789)
> X <- matrix(rnorm(100000*3),ncol=3) + matrix(rnorm(100000*3,mean=3,sd=2),ncol=3)
>
> ns <- c(5,4,3)
> system.time(mult.hist.out <- multi.hist(X,ns = ns))
   ユーザ   システム       経過  
    21.665      0.346     21.797

これでは困る。R にはちゃんとした関数 xtabs があるのだから，それ使おうヨ。

ただ，xtabs は，元のデータがカテゴリー化されているというのが前提なので，カテゴリー化は cut 関数を使って自分で行う。

もとのプログラムも cut 関数で使う breaks と同じ働きをする変数の設定はしているが，ちょっと小細工をしている。しかし，「*0.0001」はいいかげんなものだ。意味不明のオマジナイノジュモンは使わない。

    if(is.null(bins)){
        range.X <- apply(X,2,range)
        range.X[2,] <- range.X[2,] + min(range.X[2,]-range.X[1,])*0.0001
        bins <- list()
        for(i in 1:length(X[1,])){
            bins[[i]] <- seq(from=range.X[1,i],to=range.X[2,i],length=ns[i]+1)
        }
    }

正統的なプログラムにしよう。
div 関数が元のデータのカテゴリー化を行う。cut 関数で使う breaks は，一番最初の区間の左端は含まれないので（マニュアル嫁） seq(range.x[1], range.x[2], length=ns+1) だけだと最小値の行き場所がなくなる（NA になる）ので，.Machine$double.eps を引いてやる。

mx <- function(x) {
    x <- data.frame(x)
    div <- function(x, ns) {
        range.x <- range(x)
        breaks <- seq(range.x[1], range.x[2], length=ns+1)
        breaks[1] <- breaks[1]-.Machine$double.eps
        return(as.integer(cut(x, breaks=breaks)))
    }    
    return(xtabs(~., mapply(div, x, ns)))
}

このプログラムの実行時間は

> system.time(mx(X))
   ユーザ   システム       経過  
     0.810      0.034      0.842

となり，26倍ほど速い。

自分でプログラムするにしても，もう少し速いプログラムを書ける（ようになってネ）。

ちなみに，元記事のプログラム名は multi.hist となっているが，hist はヒストグラム（度数分布図）。
度数分布表は frequency table，クロス集計表は contingency table とかいうから，名前を変えたほうがよい。

 

xtable の hline.after と add.to.row

booktabs=FALSE（デフォルト）では hline.after で，どこでも，何本でも \hline を引くことができる。
booktabs=TRUE にして，hline.after では，複数の水平線（toprule, midrule, bottomrule）を引こうとしてもうまくいかない。

正解は，add.to.row を使うこと（たぶん）。

print(xtable(x), add.to.row=list(list(0:5, 4), c("[-1.5mm]\n ", "\\midrule\n")) )

これにより以下の LaTeX ソースが生成される。

\begin{table}[!Hhtbp]
\begin{center}
\begin{tabular}{rrrrrl}
  \toprule
 & Sepal.Length & Sepal.Width & Petal.Length & Petal.Width & Species \\
  [-1.5mm]
  \midrule
1 & 5.10 & 3.50 & 1.40 & 0.20 & setosa \\
   [-1.5mm]
 2 & 4.90 & 3.00 & 1.40 & 0.20 & setosa \\
   [-1.5mm]
 3 & 4.70 & 3.20 & 1.30 & 0.20 & setosa \\
   [-1.5mm]
 4 & 4.60 & 3.10 & 1.50 & 0.20 & setosa \\
   [-1.5mm]
  \midrule
5 & 5.00 & 3.60 & 1.40 & 0.20 & setosa \\
   [-1.5mm]
  \bottomrule
\end{tabular}
\end{center}
\end{table}

Mountain Lion 用の graphviz その2

というのは，必要ないようである。

そもそも，Mountain Lion にしてから graphviz（dot コマンド） が動かなかったのだが，エラーメッセージを落ち着いて読めば，

> pathDiagram(sem.dhp, "example1")
Running  dot -Tpdf -o example1.pdf  example1.dot
dyld: Library not loaded: /usr/lib/libltdl.7.dylib
  Referenced from: /usr/local/bin/dot
  Reason: image not found

ということで，libltdl.7.dylib がどこにあるか調べたら，いろんな所にあるけど，もともとは（?）Xcode がインストールするファイルのようで，ほかにも /opt/local/lib/ の中にも libltdl.7.dylib (2011/01/17) があるので，それを /usr/lib/ にコピーしてみたら，

sudo cp /opt/local/lib/libltdl.7.dylib /usr/lib/

ちゃんと動いた。しかし，Xcode の中にある最新 (2012/07/26) の

/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.7.sdk/usr/lib/libltdl.7.dylib

は，graphviz とマッチングが取れていないようで，以下のようなエラーメッセージを吐く

> pathDiagram(sem.dhp, "example1")
Running  dot -Tpdf -o example1.pdf  example1.dot
dyld: Library not loaded: /usr/lib/libltdl.7.dylib
  Referenced from: /usr/local/bin/dot
  Reason: no suitable image found.  Did find:
    /usr/lib/libltdl.7.dylib: mach-o, but wrong filetype
    /usr/lib/libltdl.7.dylib: mach-o, but wrong filetype