wrong, rogue and log:Stat

カテゴリ：Stat

• 一様乱数の積は一様乱数にはならない
  [ 2011-06-14 23:07 ]
• Rにおけるgrepはgreprが便利
  [ 2010-12-16 23:58 ]
• Incanter (+ Clojure)で遊び中…
  [ 2010-10-30 15:39 ]
• Incanter
  [ 2010-10-28 21:40 ]
• Rcppのチュートリアルビデオ
  [ 2010-10-28 21:20 ]
• Rは好きだ、しかし…
  [ 2010-09-14 11:10 ]
• RクリエイターRoss IhakaのRへの諦念
  [ 2010-09-14 00:10 ]
• Andrew GelmanのRへの不満
  [ 2010-09-09 19:28 ]
• Rcppはムッチャ速い
  [ 2010-09-08 18:48 ]
• パーレンよりカーリーブラケットの方が速い
  [ 2010-09-06 22:12 ]

一様乱数の積は一様乱数にはならない

Twitterのタイムラインに最近も出ていたのだが、1年位前にRコミュニティで流行ったこれは面白い。

Understanding “randomness”
http://stackoverflow.com/questions/3956478/understanding-randomness

I can't get my head around this, which is more random?
rand()
OR
rand() * rand()
I´m finding it a real brain teaser, could you help me out?

これは質問者は"more random"と言っているように、この質問の時点では確率変数の概念を全く理解していなくて、その意味で教科書的というか額に入れて飾りたいようなナイスな間違えの一つである。一般的に普通のプログラム言語のrand()は擬似一様乱数を生成する関数であるのだが、一様乱数を二つ掛けたらといって、もっとスゴイ一様乱数になるわけではない。

この質問のすぐ下にそれなりの模範回答がある。ただ、その回答は回答者のオリジナルではなく、Eric W. Weissteinの"Wolfram Math Worldの"Uniform Product Distribution"に添付されているMathematicaのノートブックを孫引きしているのだ。

http://mathworld.wolfram.com/UniformProductDistribution.html

しかし、結構この説明は明快なので、同じような説明をRでやってみることにする。まずは0 < x < 1の一様乱数をプロットしてヒストグラム描くことにする。ここで擬似一様乱数はrunif()によって与えられる。あと蛇足だが、Rのrunif()は最初からMersenne-Twisterで乱数を生成するので、それなりに信頼のおける擬似乱数生成をデフォルトで行う。

この擬似一様乱数で生成された確率変数のヒストグラムをみると、当たり前だが「一様乱数」の名前通り目で見て一様な分布になっている。

＞ u1 <- runif(10000)
＞ hist(u1, prob = TRUE)

つぎに、一様乱数の積をとったもの、つまりStackOverflowの質問の rand() * rand()について擬似乱数から生成した分布を見てみる。

＞ u2 <- runif(10000) * runif(10000)
＞ hist(u2, prob = TRUE)

このヒストグラムを見ても、一様な分布ではないことがわかる。実はこの分布こそが、さきほどのWolfram MathWorldにある"Uniform Product Distribution"なのである。まずは、n=2のときのUniform Product Distributionとして、- log(y)という関数を上の一様乱数の積のヒストグラムに重ねて書いてみる。

＞ y <- seq(0, 1, .01)
＞ lines(y, -log(y))

これをみると明らかなのだが、一様乱数の積の密度分布は- log(x)になるのだ。一様分布ではなく。一般的には一様な乱数をとる確率変数X_1X_2...X_nの積の密度分布の関数型は以下のようになる。



n = 2にすると-log(x)になる。この一様確率変数の積がどうしてこのような形になるのかというのは、Wolfram MathWorldにあるようにデルタ関数の積分を考えればいいのだが、ちゃんと議論するには結構大変なので、ここでは割愛する。少し検索してみると、この一様確率変数の積の密度分布について、最初に定式化を与えた記述は以下の本にあるそうだ。

Springer, M.D., 1979. The Algebra of Random Variables. Wiley, New York

また、この密度分布の関数型をより一般的な形にしたものが、産能大の石原辰雄氏によって導かれていて、もしかしたらこちらの数学的帰納法で示した証明のほうが簡単かもしれない。

「独立で非同一な分布に従う一様確率変数の和および積の分布」
http://ci.nii.ac.jp/naid/110001878196

更には、以下のペーパーは別の証明を与えている。

Dettmann & Georgiou "Product of n independent Uniform Random Variables"
http://www.maths.bris.ac.uk/~macpd/georgiou/ProductRVs%20revised.pdf

上記の例は確率変数の「積」だけれど、「和」はどうかというと、独立な一様乱数をもつ確率変数の和をたくさん重ね合わせると中心極限定理によりその重ねあわせた分布は正規分布になる。というか、別に一様乱数の確率変数でなくても、重ね合わせる確率変数の数が十分大きく、独立かつほぼ同一の分布なら正規分布になる、というのが中心極限定理の直感的理解だ。以下は一様乱数を確率密度分布とする5つの独立な確率変数の和の例。

＞ u3 <- runif(n) + runif(n) + runif(n) + runif(n) + runif(n)
＞ hist(u3, prob = TRUE)
＞ curve(dnorm(x, mean=mean(u3), sd=sd(u3)), 0, 5, add = TRUE)

実線の曲線が正規分布であるが、この5つ程度の確率変数の和であっても中心極限定理の雰囲気は味わうことができると思う。

さてこの議論の着地点であるが、実は一様乱数から任意の確率変数を生成する「逆変換法」と非常に密接に関わっているのだ。今の時代は数学ソフトウェアや統計計算環境が充実しているので、主要な確率密度分布については予め組み込みであり、ほとんどの場合はあまり余計な心配をしなくていい。しかし、その組み込み関数に存在しない確率密度分布の確率変数を生成したいというときに、「逆変換法」を用いると一様乱数をネタにして必要な確率変数を生成することができるのだ。

この「一様乱数からの確率変数の生成」ということから考えると、StackOverFlowの質問にある"rand() * rand()"というのはある確率密度（この場合-log(x))を持ったあらたな確率変数を生成しているのであって、「もっとすごい一様乱数」を生成しているのではないということなのだ。

長くなったので、このあとの解説は次回に回す。

Rにおけるgrepはgreprが便利

Rを使っていてムチャクチャにストレスがたまることの一つにgrepがある。Rのgrepはベクトルをつかうことを前提で返り値がベクトル位置だから、普通のPOSIXのgrepのようにマッチした文字列だけを欲しい場合、その手続があまりにタコすぎて嫌になる。例えば、こんな風。

> astr <- strsplit("thisistest123456789", "")[[1]]
> pos <- grep("[0-5]", astr)
> astr[pos]
[1] "1" "2" "3" "4" "5"

文字列をstrsplitでキャラクタのベクトルにスプリットしてから、先頭要素をとりだし、それをgrepして位置を割り出し、最後のその位置の文字列を求めるという結構面倒な手順。

こんなときには、stringrライブラリを使うと一発で解決する。

> library(stringr)
> str_extract_all("thisistest123456789", '[0-5]')
[[1]]
[1] "1" "2" "3" "4" "5"

これはStackOverflowから取得したネタ。ホントStackOverflowは少し時間が空いたときに眺めるには、最高のサイトだと思う。

via:
regex - How can I use grep with parameters in R ?
- Stack Overflow
http://bit.ly/f46GWo

Incanter (+ Clojure)で遊び中…

前の記事で紹介したIncanterだが、未実装部分がかなりあるものの、遊んでみるとやはり素晴らしい。

例えば、行列の実装だけれど、以下のwikiの記事をみればわかるけれど、Parallel ColtのDense Matrixで実装されていて、Rの行列よりも効率が良い。

liebke / incanter: Matrices
http://github.com/liebke/incanter/wiki/matrices

やっぱりParallel Colt等のJavaの定評あるライブラリ資源を直接使えるのはスゴイよな。Clojure自身がコンカレンシーのマシーナリーを持っているし、HadoopなどのJavaの数ある並列インフラを「ライブラリ」レベルで利用できる。また、ClojureはGoogle AppEngineで動作するようだから、Incanterも動作するかも。

更には、JVM上で動くからCほどじゃないにしても、JITの最適化により、Python, Ruby, Rなどでは敵わないほどの実行速度なのは素晴らしいと思う。

豊富で実践配備された性能抜群のライブラリとJVMというJavaの資産が利用可能で、動的型付け言語＋インタプリタだから足の早いインクリメンタルな開発が可能な上にコンカレンシーも可能になっている（ここまではClojureの機能）インフラに、統計解析ライブラリを実装するというスジの良さにはメロメロだ。

Incanterに対する懸念点は、Rほどにライブラリが増えるかどうかがわからない、というところかな。それは、Lispのカッコと前置記法に世の中の研究者や実務者が慣れるかどうか未知数ということと同義だと思う。まあ、Emacs Lisp使っている人は結構いるので、普及の可能性はかなりあるとは思うけれど。

統計解析関係のIncanterでの取り扱いは、以下のドキュメントあたりから、ページを手繰ってみると雰囲気が掴める。非常にRから影響を受けていることがわかる。

Introduction to Incanter
http://data-sorcery.org/2009/05/31/introduction-to-incanter/

いずれにしろ、もうしばらく遊んでみるつもり。

Incanter

いや、もうこれヤバすぎ。ClojureベースのRクローン。

Incanter
http://incanter.org/

巨大データもコンカレンシーもなんでもござれ。なにせJVMでClojureだから。Javaの膨大なライブラリがシームレスに使える。もうこれはコンセプトだけで萌えまくる。ちょっと動かしてチュートリアルをやったけれど、ホントに動いて感動した。

まだRクローンというには機能的な実装が足りないようだけれど、そこはlisp方言。一旦火がついたらとんでもない速度で開発が進むと思う。ライブラリリポジトリはCIANとか名づけられるのだろうか。

スゴイ楽しみだ！

Rcppのチュートリアルビデオ

Rcppの作者Dirk Eddelbuettel(とRomain Francois)によるRcpp, RInside、そしてRProtoBufのチュートリアルビデオ。Google Tech Talksでの講演。もちろんinlineパッケージについての言及がある。

Integrating R with C++: Rcpp, RInside, and RProtobuf
http://www.youtube.com/watch?v=UZkaZhsOfT4



RProtoBufはGoogleの有名なProtocol BuffersにRインタフェースをつけたもの。これはRからのシリアライザとして是非とも使うべきだな。

RInsideはC++のコードにRのコードを埋め込むというもの。まだこちらはあまり洗練されているわけでなく、コッソリ裏でRでEvalした結果を返すという実装のよう。

それにしてもDirk Eddelbuettelの作るR関連のコードはハッカー視点のセンスが素晴らしいと思う。

Rは好きだ、しかし…

最近、Rの非効率性やパフォーマンスの悪さをここで良く取り上げていて、それじゃあ柏野はRが嫌いなのか、と誤解されることもあると思うので、その弁明を。

僕にとってRはかなり好きなコンピュータ言語だ。使っていてワクワクする。言語を使いながらマイナーな分野の統計学の勉強ができるフリーソフトウェア / オープンソースソフトウェアの言語なんて他には存在しない。最近ではCRANにあるevirやismev, evdを利用して極値統計学を勉強した。

クイックなプロトタイピングができるのも本当に素敵だ。よくわからないデータを探索的に調べていくためには最高のブッシュナイフだと思う。だからこそ、@ITで市井の知識労働者を対象とした連載をしている。

実践！ Rで学ぶ統計解析の基礎
http://www.atmarkit.co.jp/fcoding/index/stat.html

しかし、エンタープライズレベルに必要な最低数百GB、普通に数TBのデータを目の前にしたときに、これをどうするかということになると、Rではどうしょうもないと思う。竹槍で核爆弾に立ち向かう気分になる。全データをメモリに一度読み込んでから処理する仕様、変数割り当てや参照のたびにコピーが発生するという非効率性、名前空間やスコーピングの不備による変数の安全性問題と並列化への非対応、カッコの種類や数によって速度が変化する絶望的なパーシングの遅さを誇るインタプリタ、長いコードを書いていくと可読性が悪くなる奇妙なシンタックス、救いようがないエクセプションハンドリング、どれをとってもお金を稼ぐ現場で使うモダンな言語としての資格を備えていない。

GoogleのHal Varianも言うように、これからは"Data matters"なのだ。

Why Data Matters
http://googleblog.blogspot.com/2008/03/why-data-matters.html

ただし、その場合のデータとは、今までのような数百MB程度のデータであることは少なくなっていくだろう。Revolution Analyticsのように、Hadoopのような並列環境に強引にRを載せて、数TBのデータクランチングを行わせようという試みもあるが、筋の悪い言語設計をシステム実装でカバーするのは、将来に多大なツケを産むことになるだろう。僕自身も当初はRを並列化プラットフォームに載せて処理をするということを考えたけれど、いくつかの試みの後でそれはもうあきらめた。現在すぐに利用するならば、必要な統計関数は自前で実装するしかないと思っている。ただし、現在ではどのプラットフォームもどの言語も数学ライブラリが充実しているので、それを組み合わせる作業になるだけだからそれほどドローバックはない、という素敵な事実が背景にあるから、そう決断できるのだが。

Rはこれからも教育用の言語として生き残るだろう。また、少人数の研究者が自ら集めた比較的小さなデータに対してデータ解析・統計解析をするときには、切れ味鋭い道具であり続けるだろう。しかし、日々生まれる巨大データとこの道具で格闘して世界を変えるのは難しいというのが厳しい現実である。Rを創始した一人であるRoss Ihakaが自らRの末期ガンの宣告をしたように、下手な延命はせず、役に立つ分野ではそのまま役に立たせ、時代の要請に沿わなくなったら、そのまま逝かせてやるほうがいいと思う。そして、その屍から新たな生命が生まれるのだ。

---
コメントありがとうございます。

私も基本的にはデータを前処理するところと、解析するところ、そしてプレゼンテーションをするところは分離されていくと思っています。そしてだからこそRのような中途半端に全てを網にかけている立ち位置は、エンタープライズの世界では使われにくくなっていくのだと思います。そしてdata squashingのような技術が発達すれば、そこである程度の解析までできてしまいますから、ますます統計解析言語/環境としてRの役割は薄れていくと思います。もちろん、Rの機能としては、ますますプレゼン機能を充実させるという進化のさせかたもあります。現在別アプリであるGGobiを取り込むのはかなり良いアイディアかもしれません。

それにしても、RのS3/S4メソッド混在という文法上の不整合や、スコーピングやネームスペースの不在、そして非効率な実行速度をどうにかしないと、破綻するということは変わりないと思います。

RクリエイターRoss IhakaのRへの諦念

なんともまあ、RのクリエイターであるRoss Ihakaも、Rのあまりの非効率的でノロマな実装に

I've come to the conclusion that rather than "fixing" R, it would be much more productive to simply start over and build something better.

と言い放っちゃった…。結局、Sからの負の遺産や最初の言語設計の不十分さが積もり積もって、もうチョットやそこらの修正では、Rの一連のパフォーマンス問題や非効率性の解決はできないみたいだ。

http://xianblog.wordpress.com/2010/09/06/insane/comment-page-4/#comment-5187

IhakaのRに対する不満は大きく二つあって、一つは、Rはスカラーの型をもっていないのとRのインタプリタ実装が激遅なので、スカラー演算がアホみたいに遅いこと、もう一つは、データフレーム等では値参照をしているのでそのたびにコピーが発生して死ぬほど遅くなるということ。その他にもスコーピングが実装されていないので、あるコーディングでは変数がグローバル変数になったりローカル変数になったりというような狂ったような仕様も不満みたいだ。まあ当然ですね。特にスコーピングがダメだということは、言語としてはマトモな並列化はできないということだから、進化の袋小路に入ってしまっていると思う。

Ihaka自身は、Rを今以上に効率化したり修正したりする作業について、労多くて得るものが少ないと見切っているみたい。それよりもいっそのこと全て新しく作り直すほうがいいと言っている。そしてその場合、その新しい統計言語は、Rとはセマンティックスは変わるだろうし、マルチコアや並列計算に対応しなければならないとも言っている。

いやー、一々がごもっともだ。

Andrew GelmanのRへの不満

カッコによる速度差の一連の議論から派生して、Andrew GelmanがRに不満を述べている。

The future of R
http://bit.ly/dpmNIB

I don't know that R's library is so brilliant as all that--if necessary, I don't think it would be hard to reprogram the important packages in a new language.

I would say, though, that the problems with R are not just in the technical details of the language. I think the culture of R has some problems too. As I've written before, R functions used to be lean and mean, and now they're full of exception-handling and calls to other packages. R functions are spaghetti-like messes of connections in which I keep expecting to run into syntax like "GOTO 120."

ほんとだよね。プロのエンジニアから言わせてもらうと、RのSyntaxはタコそのもの。そしてなによりエクセプションハンドリングがヒドすぎ。これで複数人で少し巨大なコード書いたら絶対にオカシクなると思う。今はメインのユーザが研究者で、研究者は自分の研究のところしか興味ないからあまり感じないのだろうけれど。

前エントリにも示したように、カッコの種類や数でパフォーマンスが変わるなんて話をしたら、聞いていた人間がため息ついていたよ。普通の言語では信じられない事態だけに、チョットこれはいただけないと思う。

まあ、このように言語そのものには不満があるのだが、やっぱり組込みの統計関数やCRANのパッケージ群はムチャクチャ魅力であるのは間違いない。だから、個人が探索用の「ブッシュナイフ」として使うのがRを利用する一番良い使い方で、パフォーマンスが必要な場合やミッションクリティカルな場合は別実装をするべきなんだと思う。Googleはまさにそうしているよね。

Rcppはムッチャ速い

Debian Projectに所属しているRcppの開発者Dirk Eddelbuettelの超高速コードである。先日のRadford Nealのエントリを受けたChristian Robertのブログエントリをみて、Rcppで実装してみたとのこと。

Straight, curly, or compiled?
http://dirk.eddelbuettel.com/blog/2010/09/07/#straight_curly_or_compiled

こちらもrbenchmarkをインストールして、そのまま試してみよう。

> 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}
> library(Rcpp)
> library(inline)
> l <- cxxfunction(signature(ns="integer", xs="numeric"),
+ 'int n = as(ns); double x=as(xs);
+ for (int i=0; i+ return wrap(x); ',
+ plugin="Rcpp")
> library(rbenchmark)
> N <- 1e6
> benchmark(f(N, 1), g(N, 1), h(N, 1), j(N, 1), k(N, 1), l(N, 1),
+ columns=c("test", "replications", "elapsed", "relative"),
+ order="relative", replications=10)
test replications elapsed relative
6 l(N, 1) 10 0.092 1.0000
5 k(N, 1) 10 13.469 146.4022
1 f(N, 1) 10 13.929 151.4022
4 j(N, 1) 10 15.671 170.3370
2 g(N, 1) 10 16.319 177.3804
3 h(N, 1) 10 20.774 225.8043

lがRcpp実装だが、ハヤっ！このくらい速いと、カーリーブラケットがどうの、パーレンがどうのという問題なんてなくなっちゃうな。なるほどね、inlineパッケージのcxxfunctionを使うと、メチャクチャ簡単にRcppが利用できるんだね。これは素晴らしい！これからRcppを使うことに決めましたよ。

パーレンよりカーリーブラケットの方が速い

Tronto大のRadford Nealのブログエントリである。

Two Surpising Things about R
http://bit.ly/bQrOTk

... But it turns out that in R you can use both parentheses and curly brackets! The curly brackets are normally used to group statements, but an expression is one type of statement, and the last (or only) in the group provides the value.
...
Here’s the even more surprising thing. It occurred to me that before rushing out and using this feature, I should check that it doesn’t introduce some horrible inefficiency, as might be the case if curly brackets were optimized for their more common use in grouping statements.

Rの数式表現にはカーリーブラケット { }もパーレン( )もどちらも使えるのだが、カーリーブラケットの方が計算が速くなるそうだ。

早速僕のiMacでも試してみよう。

> a <- 5; b <- 1; c <- 4
> f <- function (n) for (i in 1:n) d <- 1/{a*{b+c}}
> g <- function (n) for (i in 1:n) d <- 1/(a*(b+c))
> system.time(f(1000000))
ユーザ システム 経過
2.073 0.023 2.137
> system.time(g(1000000))
ユーザ システム 経過
2.268 0.017 2.342

ほんとだ。0.2sec違う。Radford Nealの結果と同じくらいだ。へえ。今度からカーリーブラケットでコーディングするか。

あと、これはR 2.11.1のr52822という特定バージョンのモンキーパッチになってしまうが、同じくRadford Nealが作成した、高速化のためのパッチが以下にある。

http://www.cs.toronto.edu/~radford/R-mods.html

もちろんパッチにあてたあとに./configure, makeが必要だ。