3D Surface Plot with lattice/rgl

slashdot を眺めていると hirax さんによる Mathematica を利用した曲面プロットのチュートリアルが紹介されていた。 Excel や Maxima には移植されているのに、残念なことに R には移植されていない模様。 ということで R と lattice/rgl(OpenGL) で描いてみた。 — lattice Version lattice は Bell Labs で開発された S/S-PLUS 向け可視化フレームワーク Trellis Graphics を R で実装したもの。 2004年には John Chambers Award を受賞している。 以前書いたコードをそのまま流用したのが以下 expand.grid {base} は 直積を作ってくれる関数。 — rgl Version rgl は R で OpenGL を呼び出すパッケージ。 2003年には John Chambers Award

R fortunes

R 2.14 をインストールしたおりにふと Rprofile.site ファイルを見ると、これまでよりファイルサイズが大きくなっていた。 中身を確認すると、 CRAN ミラー設定の後ろに次のようなコメントがあった。 R 環境でも Unix/Linux の fortune が提供されるようになったのかと、コメント通りに install.packages(‘fortunes’) で fortunes パッケージをインストールし、コメントアウトを解除し、R を再起動すると予想通り fortune メッセージが表示された。 現時点(1.4-2)では 292 個のメッセージがある。イースターエッグを見つけたみたいですごく嬉しい。

Set Environment Variables in R

R で環境変数を設定 R で環境変数を設定するには環境設定ファイル内でKEY=VALUE形式で指定するとよい。 環境変数は次の3つのレイヤーで指定できる。 user specific site specific system specific Windows では system レベルの指定が出来なかったり、OS によりパスが異なる。 また、オプション(–vanilla/–no-environ)を与えて特定の環境変数ファイルを読み込まないようにすることもできる。 Windows/Unix 環境を対象にこれらをまとめたのが次の表。(click to enlarge) パッケージのサーチパスを環境変数で定義 R でマイナーアップグレードする際(たとえば 2.13.1  から 2.13.2)、マイナーバージョン間でパッケージのサーチパスを共有しておけば、インストール済みのパッケージを一から再インストールしなくてすむ。 環境変数 R_LIBS でパッケージのサーチパスを設定できるので、この値をRの環境変数設定ファイルにKEY=VALUE形式で記載すればよい。 D:\R\R-2.13.2>type etc\Renviron.site R_LIBS=”D:/R/lib-2.13″ 関連ソースコード src/main/Renviron.c::process_user_Renviron src/main/Renviron.c::process_site_Renviron src/main/Renviron.c::process_system_Renviron src/unix/system.c::Rf_initialize_R src/gnuwin32/system.c::R_SetWin32/cmdlineoptions

Bland–Altman Plot in ggplot2

Peter Dalgaard 著 “Introductory Statistics with R” の Ch. 5.6 “Paired t test” を読んでいたところ「対応のある2群の差を可視化して確認するには、対応同士の平均と差の分散図を描くとヨイ」といったようなことが書かれていた。このプロットを Bland–Altman Plot と呼ぶらしい。 せっかくなので ggplot2 で描いてみた。 bland_altman_plot <- function(dat){   dat$diff <- dat$pre – dat$post   ggplot(dat, aes(pre, diff)) +   geom_hline(yintercept = mean(dat$diff) + c(-2, 0, 2) * sd(dat$diff), linetype=2, color=’brown’) +   geom_point() +  theme_bw() +   xlab(‘mean of

Move To Front Rule In Action(amazon.com)

服部哲弥先生の「Amazonランキングの謎を解く―確率的な順位付けが教える売上の構造 (DOJIN選書)」(2011年。化学同人社)という本を読んでいたところ、Amazonの売り上げランキングや2chのスレッドランキングが move to front rule というモデルで解説されていた。注文が発生すると、最大瞬間風速的にランキングは1位になる。その後は、他の本の注文が発生するたびにずるずると順番を下げていき、再び注文があったタイミングで1位に返り咲く。 書籍はアマゾンランキングの一般論がテーマなので、具体的に書籍を特定したうえでその本のランキングをプロットというのは行わず、move to front rule に基づく一般的なランキングの遷移がプロットされているだけであった。 やはり、実際の本の動きが気になるので、amazon.com/amazon.co.jp でランキングが上位の本と数十万台をさまよっている本を選び、1時間ごとのランキングを約1.5ヶ月間(2011年6月上旬から)にわたりプロットしたのが以下 (画像をクリックすると amazon サイトに飛びます) amazon.com Best Sellers Average Sellers amazon.co.jp Best Sellers Average Sellers MEMO ランキングが10万台の本は確率過程(stochastic process)のテキストがおてごろだったため、利用させていただく。 cronで1時間ごとに実行していたが、LANの問題、プログラムの不具合、その他により、一部データが欠落 今はランキングを plot しているだけなので、もう少し何とかしたいところ。 眠いので、今日はここまで。

rworldmap:Mapping Worldwide Suicide Rates

R Journal Volume 3/1, June 2011 で位置データを可視化する “rworldmap” という R パッケージが紹介されていた。 Andy South : rworldmap: A New R package for Mapping Global Data. The R Journal Vol. 3/1, June 2011. ISSN 2073-4859 rworldmap は R プログラムが本職でない人でもお手軽に使えるよう設計されており、同解説記事では rworldmap のビジョンが次のように宣言されている。 The vision for rworldmap is to produce a package to facilitate the visualisation and mapping of global

Convergence of E

和達　三樹，十河　清著の「キーポイント確率統計 (理工系数学のキーポイント)」という本にでていた自然対数の底 e の近似計算（P.51) (4.29)、 (4.30)の2つの近似が紹介されているが、e への収束の速さは大きく異なる。 収束の速さを R でアニメーションにしたのが以下のコード(plot 関数バージョンと ggplot2 バージョン) plot version e1(green) : (4.29) e2(red) : (4.30) base.plot<- function() { plot(1:3, 1:3, type=”n”, asp=1, ann=FALSE, axes=’FALSE’, frame.plot=TRUE) title(main=”A peek inside the \n estimation of e”) axis(side=2, at=c(seq(1, 3,by=0.25)), labels=c(as.character(seq(1, 3,by=0.25))), par(cex=0.8,las=1)) } e.est <- function(iter) { sum.e1<- sum.e2 <- stor.rat<- numeric(0)

Move To Front Rule In Action(2ch.net)

服部哲弥先生の「Amazonランキングの謎を解く―確率的な順位付けが教える売上の構造 (DOJIN選書)」(2011年。化学同人社)という本を読んでいたところ、Amazonの売り上げランキングや2chのスレッドランキングが move to front rule というモデルで解説されていた。注文が発生すると、最大瞬間風速的にランキングは1位になる。その後は、他の本の注文が発生するたびにずるずると順番を下げていき、再び注文があったタイミングで1位に返り咲く。 move to front rule を満たし、アマゾンのランキング以上に直観的に理解できる身近な例として２チャンネルのスレッドランキングがある。本の品格を保つためか(?)後ろの方(Ch7.3 P.174-)でちょろっと紹介されているだけだった。 数学的な話はかなり高度で難解なのだが、ランキングの時系列遷移はデータをとってプロットすればだれでもできる。 ２ちゃんねるの板の中からアクティブレベルに応じて4つを選び、各板での1時間ごとのランキングを約1ヶ月間(2011年6月)にわたりプロットしたのが以下 ランク変動が活発な板 ニュース速報＋ GIF アニメーション   ふつうの板 プロ野球 GIF アニメーション のんびりした板 プログレ GIFアニメーション 宝塚、四季 GIF アニメーション MEMO R の animations パッケージをつかい、Windows 環境固有の凶悪な仕様とバトルしつつ swf 形式で作成したものの、WordPress にアップロードできないことに気づいたり、mp4 変換してYouTubeにアップロードしたら、アップロードが完了した瞬間に “This video has been removed because it is too long.” と削除されたり(アップロード後に動画が40分近くあることに気づく)と紆余曲折を経て、レガシーな gif アニメーションに落ち着いた。 とはいえ、何メガもあるgifファイルを何個ものせると重いので、ffmpeg で動画変換して

Café Wall Illusion In gplot2

どこかのサイトを巡回中に Café Wall Illusion という錯視を見つけた。 Bristol にあるカフェの外壁に錯視が利用されているのを Bristol 大学の研究者が 1973年に気づき、1979 年に論文にまとめられた([1])。 Figure 1. The original of the Cafe Wall, St Michael’s Hill, Bristol. ([1] から) みるからに R でさくっと描けそうだったので R & ggplot2 で描いてみた — Program 処理の流れとしては 格子をつくり 黒のマス目を決定し 白地に黒を塗り 横線を引けば完了 l1 <- expand.grid(x=c(seq(0, 4) * 4 , seq(0, 4) * 4+1), y=seq(0, 2) * 4) # first

Setting The Default CRAN Mirror

パッケージをインストールしたりアップグレードするときに install.packages などとやると CRAN のミラーの選択画面が表示される。いつも同じミラーを選ぶのならデフォルトで設定して、選択の手間を省きたい。 CRANミラーの管理 R パッケージのミラー情報はオプションの ‘repos’ で管理されている。 デフォルトでは CRAN ミラーが “@CRAN@” になっているおり、パッケージ管理時にデフォルトのままだとミラー選択に誘導されるため、これをお好みのミラーサイトに設定すればよい。 Option の getter/setter R オプションの getter/setter は次のようになっている。 getter : getOption(‘key’) setter : options(key=value) 今回の CRAN ミラーの場合、次のようになる getter : getOption(‘repos’) setter : options(repos=foo) Windows 環境でのミラー設定 オプションは次のファイルで設定する R_HOME/etc/Rprofile.site R_HOME は環境変数で、実際のパスは “Sys.getenv(‘R_HOME’)” で確認可能。 上記ファイルには次のようなひな形があるので、”http://my.local.cra&#8221; の箇所をお好みのミラーに変更し、コメントアウトすればよい。 # set a CRAN mirror # local({r <-