R 2.13.0についてのポストを一部翻訳してみた
via http://dirk.eddelbuettel.com/blog/2011/04/12/#the_new_r_compiler_package
意味が通ってない気がするところは後で直すかもしれません。すみません。
2011年4月11日 新しいR2.13.0のcompilerパッケージについての最初の実験
R 2.13.0が明日リリースされます。新しいバージョンには、RコアメンバのLuke Tierneyによる新しいパッケージcompilerが含まれます。この極めてシンプルかつ直接的な名前のパッケージは、Lukeが何年もの間取り組んできたものです。NEWSファイルには以下のように記載されています。
- compilerパッケージがついに標準パッケージとなりました。使い方を知るには?compiler::compileコマンドを利用してください。このパッケージはRのためのバイトコードコンパイラの実装です。コンパイラはデフォルトでは本リリースには用いられていません。baseパッケージやその他の推奨されるパッケージをコンパイルする方法については、'RInstallation and Administration Manual'を参照してください。
Rのbaseパッケージ、その他推奨パッケージにはcompilerは利用されていないということです。
R/Finance 2011に先立って開催されるRcppワークショップに向けてスライドを作っているときに、compilerの良い例題を思いつきました。昨年の夏と秋に、Radford NealはRのパフォーマンス改善のための詳細な改善案(possible patches)をdevelopment listに送ってくれました。
Some patches were integrated and some more discussion ensued. One strand was on the difference in parsing between normal parens and curly braces. In isolation, these seem too large, but (as I recall) this is due some things 'deep down' in the parser.
しかしながら、この問題については何人かから強い要望として挙げられていました。
And it just doesn't die as a post from last week demonstrated once more.
去年、Christian Robertはこの問題をいくつかの必要な機能群に切り分けてくれました。私は、僅かな高速化のためにコードの10%以上を書きかえるようなことに時間を費やすくらいなら素直にRcppを使う、とpostしました。
ですので、ここでcompilerパッケージがどれだけの仕事をしてくれるのか試してみたいと思います。出発点は昨年と同じです。明示的なforループの中に、1/(1+x)を計算する5つの変数(variant)があります。実際のコードは、これらを効率化するためにベクトル化をしてくれるということはありません。ただし、Christianの議論に従って、パーサによる速度差の発生について言及しておくことは無意味ではありません。
とても便利なrbenchmarkパッケージを使って計測したサマリーは以下です。
> ## cf http://dirk.eddelbuettel.com/blog/2010/09/07#straight_curly_or_compiled
> f <- function(n, x=1) for (i in 1:n) x=1/(1+x)
> g <- function(n, x=1) for (i in 1:n) x=(1/(1+x))
> h <- function(n, x=1) for (i in 1:n) x=(1+x)^(-1)
> j <- function(n, x=1) for (i in 1:n) x={1/{1+x}}
> k <- function(n, x=1) for (i in 1:n) x=1/{1+x}
> ## now load some tools
> library(rbenchmark)
> ## now run the benchmark
> N <- 1e6
> benchmark(f(N,1), g(N,1), h(N,1), j(N,1), k(N,1),
+ columns=c("test", "replications",
+ "elapsed", "relative"),
+ order="relative", replications=10)
test replications elapsed relative
5 k(N, 1) 10 9.764 1.00000
1 f(N, 1) 10 9.998 1.02397
4 j(N, 1) 10 11.019 1.12853
2 g(N, 1) 10 11.822 1.21077
3 h(N, 1) 10 14.560 1.49119
これはChristianの結果を真似たものです。{}を使っているk()関数が最も速いことがわかります。最も遅い関数h()は、相対的に49%遅くなっています(実時間では5秒ほど異なります)。一方、普通の書き方であるf()関数はとてもよく働いてくれることが分かります。
これらの関数をcomplierで処理するとどうなるでしょうか。cmpfun()関数を用いてコンパイル後の変数(variants)を生成してから利用してみます。
> ## R 2.13.0 brings this toy
> library(compiler)
> lf <- cmpfun(f)
> lg <- cmpfun(g)
> lh <- cmpfun(h)
> lj <- cmpfun(j)
> lk <- cmpfun(k)
> # now run the benchmark
> N <- 1e6
> benchmark(f(N,1), g(N,1), h(N,1), j(N,1), k(N,1),
+ lf(N,1), lg(N,1), lh(N,1), lj(N,1), lk(N,1),
+ columns=c("test", "replications",
+ "elapsed", "relative"),
+ order="relative", replications=10)
test replications elapsed relative
9 lj(N, 1) 10 2.971 1.00000
10 lk(N, 1) 10 2.980 1.00303
6 lf(N, 1) 10 2.998 1.00909
7 lg(N, 1) 10 3.007 1.01212
8 lh(N, 1) 10 4.024 1.35443
1 f(N, 1) 10 9.479 3.19051
5 k(N, 1) 10 9.526 3.20633
4 j(N, 1) 10 10.775 3.62673
2 g(N, 1) 10 11.299 3.80310
3 h(N, 1) 10 14.049 4.72871
驚異的に速くなりました！修正はとても容易でした！使い方は簡単です、関数を用意し、コンパイルするだけです。たったそれだけで、限界と思っていたところを超えるスピードが手に入りました。さらに言えば、括弧の書き方の違いによる差が消えていることがわかります。明示的な指数をとっている関数は依然として遅いですが、これは余計な関数呼び出しが含まれてしまうためだろうと思われます。
新しいcompilerパッケージを活用できる場面は極めてたくさんありそうだと言えるでしょう。どういったケースに有効なのか（あるいは有効でないのか）について、さらに多くの人がテストしてくれることを望みます。最後に、Rcppについてもう一度触れておきたいと思います。compiledコードとbyte compiledコードの差を見ることができるでしょう。
> ## now with Rcpp and C++
> library(inline)
> ## and define our version in C++
> src <- 'int n = as<int>(ns);
+ double x = as<double>(xs);
+ for (int i=0; i<n; i++) x=1/(1+x);
+ return wrap(x); '
> l <- cxxfunction(signature(ns="integer",
+ xs="numeric"),
+ body=src, plugin="Rcpp")
> ## now run the benchmark again
> benchmark(f(N,1), g(N,1), h(N,1), j(N,1), k(N,1),
+ l(N,1),
+ lf(N,1), lg(N,1), lh(N,1), lj(N,1), lk(N,1),
+ columns=c("test", "replications",
+ "elapsed", "relative"),
+ order="relative", replications=10)
test replications elapsed relative
6 l(N, 1) 10 0.120 1.0000
11 lk(N, 1) 10 2.961 24.6750
7 lf(N, 1) 10 3.128 26.0667
8 lg(N, 1) 10 3.140 26.1667
10 lj(N, 1) 10 3.161 26.3417
9 lh(N, 1) 10 4.212 35.1000
5 k(N, 1) 10 9.500 79.1667
1 f(N, 1) 10 9.621 80.1750
4 j(N, 1) 10 10.868 90.5667
2 g(N, 1) 10 11.409 95.0750
3 h(N, 1) 10 14.077 117.3083
Rcppは依然として80から最高120倍は高速です。byte compiledコードに対して、速度差は約25-foldです。ただし、これらの例は実際のRコードに比べて非常に単純なものなので、実際の高速化具合はcompilerパッケージ、Rcpp共にもう少し小さいものになるでしょう。
Before I close, two more public service announcements. First, if you use Ubuntu see this post by Michael on r-sig-debian announcing his implementation of a suggestion of mine: we now have R alpha/beta/rc builds via his Launchpad PPA. Last Friday, I had the current R-rc snapshot of R 2.13.0 on my Ubuntu box only about six hours after I (as Debian maintainer for R) uploaded the underlying new R-rc package build to Debian unstable. This will be nice for testing of upcoming releases. Second, as I mentioned, the Rcpp workshop on April 28 preceding R/Finance 2011 on April 29 and 30 still has a few slots available, as has the conference itself.
多変量データの関係を何となく見る方法
多変量のデータがどんな関係を持っているか分からないとき、まずは生データを観察することで何となくデータの傾向が見えてくることがあります。いきなり主成分分析だとかをやってもよく分からないという場合は、統計処理をする前にデータ自体を眺めることでヒントが見えてきたりします。
たとえば、Rの組み込みデータセットの一つfgl
> fgl
RI Na Mg Al Si K Ca Ba Fe type
1 3.01 13.64 4.49 1.10 71.78 0.06 8.75 0.00 0.00 WinF
2 -0.39 13.89 3.60 1.36 72.73 0.48 7.83 0.00 0.00 WinF
3 -1.82 13.53 3.55 1.54 72.99 0.39 7.78 0.00 0.00 WinF
4 -0.34 13.21 3.69 1.29 72.61 0.57 8.22 0.00 0.00 WinF
5 -0.58 13.27 3.62 1.24 73.08 0.55 8.07 0.00 0.00 WinF
6 -2.04 12.79 3.61 1.62 72.97 0.64 8.07 0.00 0.26 WinF
7 -0.57 13.30 3.60 1.14 73.09 0.58 8.17 0.00 0.00 WinF
8 -0.44 13.15 3.61 1.05 73.24 0.57 8.24 0.00 0.00 WinF
9 1.18 14.04 3.58 1.37 72.08 0.56 8.30 0.00 0.00 WinF
10 -0.45 13.00 3.60 1.36 72.99 0.57 8.40 0.00 0.11 WinF
...こんな感じのデータになっています。ガラス質の砂利(fragments)の成分データのよう。2〜9列目が各成分の構成比率となっています。
とりあえず各2成分の組み合わせを見る散布図行列を見てみると
pairs(fgl[,2:9]) #2~9列目の散布図行列
こんな感じになっています。散布図行列で傾向が分かることもありますが、3変数以上の関係を把握するのが難しく、また各散布図間で点がどういう対応関係にあるのかもわかりません。もちろん高次元データを直交座標に映して把握できれば良いわけですが、人間は3次元以上の空間は認識できないので、それは不可能です。それでも何とかしようということで、平行座標プロットというものが考えられています。
平行座標プロットとは、その名の通り、各変数軸(座標軸)を平行に並べたものです。直交は無理なので、平行に並べるわけです。これにより、ある程度多変数間の繋がりを視覚化することができます。
fglデータに対する平行座標プロットはこのようになります。Rでは、MASSライブラリを用いて
library(MASS) parcoord(fgl[,2:9])
とすることで表示できます。簡単です。
これだけでも何となく傾向が分かると思いますが、線が重なってしまい（これは平行座標プロットでは仕方のないことですが）、目で見た傾向（Naが中程度のものはMgが大きくなる傾向があり…）が本当に正しいか不安です。傾向の違う線が途中で混ざっていても分からないからです。
そこで、平行座標プロットに色をつけることでもう少し分かりやすいものとすることを考えます。
# x: 列が変数、行がサンプルのデータ # order.no: 列についてソートした上で色付けしたい場合は列番号を指定 # whiteRow: 表示しない行番号(ソート後) # over: 下から描画したい場合(変数値の大きい側が上書き描画)はT modParcoord <- function(x,order.no=-1,whiteRow=c(),over=F){ #オーダーする変数が指定されている場合は並び変え if(order.no!=-1){ x.order <- x[order(x[,order.no]),] }else{ x.order <- x } #描画しない部分を用意 if(!is.null(whiteRow)){ x.white <- x.order[whiteRow,] x.draw <- x.order[-whiteRow,] }else{ x.white <- NULL x.draw <- x.order } #残った部分を描画 n.draw <- nrow(x)-length(whiteRow) col.alpha <- rep(1.0,n.draw) #半透明度 ここでは指定しない col.draw <- rainbow(n.draw,col.alpha) if(over){ #色順がそのままになるように描画順を逆転 x.draw <- x.draw[order(x.draw[,order.no],decreasing=T),] col.draw <- rev(col.draw) } x.draw <- rbind(x.white,x.draw) #描画しない部分(white)を先に描画 parcoord(x.draw,col=c(rep("white",length(whiteRow)),col.draw)) }
このような関数を用意して、色(グラデーション)つきの平行座標プロットを書いてみたものが次のグラフです。Caの列でソートして値順に色をつけています。
modParcoord(fgl[,2:9],6,over=T)
真ん中あたりのオレンジや黄色、緑といった線の集まりが特に傾向を持っていることが見てとれます。色なしのparcoordに比べ、傾向把握がより信頼できるものになったことが分かるかと思います。
また、modParcoord関数にはoverという引数を持たせています。これは「線を下から書くか上から書くか」を指定する引数です。parcoord関数は行の上から順に線を描いていき、後から書いた線は前の線を塗りつぶす形で描画されます。そのため先に書かれた変数の傾向が見えなくなるという問題を緩和するために、描画の順番を逆転できるようにしています。先ほどのグラフで逆転してみると
となり、かなり見え方が違ってきます。
さて、めでたく生データの傾向が見えたわけですが、実は最初に示した通り、fglデータは右端にカテゴリー変数が付加されています。データ収集者が予めカテゴライズしてくれたものです。この中で、"WinNF"カテゴリだけをプロットしてみると
となります（ここでは色をつける意味はありませんが、modParcoord関数を流用しえいるためこのようになっています）。これは先ほど見つけた「オレンジや黄色、緑」の傾向とほとんど同じものです。目視による分析が正しく傾向を見つけられていたことが裏付けられると言えます。
実際のデータではこのような適切なカテゴリ変数が与えられていないことも多いでしょう。クラスタリングや次元縮約などの方法でカテゴリを見つけていくことになるわけですが、その前に、生データを色々な角度から見ることで、結果の当たりをつけることができる場合もあるということが実感できるのではと思います。
