【統計】R

Rの環境操作とクロージャ

折を見て、Common Lisp関連の書籍を読んでいます。なかなか進まないのですが。。。
何を調べていた時か忘れてしまいましたが、こちらで次のようなコードを見つけました。

Javascriptとかでもよく見かけるクロージャのサンプルですね。
で、はたと気になりました。Lispの影響を強く受けているRではどのように書けるのかと。
ちょっと書いてみました。

Rでは、CLのletにあたる関数がありません(多分)。なのですが、そもそもRには環境自体を作成、操作できる関数があるので、そちらを使うこととします。

…期待どおりの動作をしてくれません。環境とクロージャ環境中のオブジェクトを表示してみます。


関数が呼ばれる度に新しい環境が作成されていて、counterは新しい環境内の変数となっているようです。

ということで、Rでは関数内で定義された変数の扱いがCLと違うようで、CLでは親環境にまで派生しますが、Rでは親環境にまで派生しないようです。
Rのスコーピング規則は、Ross先生がSとの互換性を捨ててまでこだわったところらしいので、このような仕様にしているのは何か理由があるのでしょう。破壊的な代入が親環境にまで影響するのが嫌だったとかなのかなー
詳しくはわからないんですが、Pythonのジェネレータっぽいことしたい時はどうすれば良いんだろう。クラス作れとかってことなんでしょうかね。

ちなみに、Rのレキシカルスコープについては、以下の書籍が詳しいです。

Rの基礎とプログラミング技法

著者/訳者:U.リゲス

出版社:シュプリンガー・ジャパン(株)( 2006-10-22 )

定価:

単行本 ( 258 ページ )

ISBN-10 : 4431712186

ISBN-13 : 9784431712183


EmacsでRのコード補完を快適にしてくれるauto-complete-acr.elがオムニ補完に対応しました。

Emacs Wiki更新情報を眺めていると、ac-R.elなるものを発見。

auto-complete.elのRインターフェイスが出てる!! “EmacsWiki: ac-R.el” – http://j.mp/cbrPMEMon Aug 23 22:48:03 via KeySnail

案の定、auto-complete.elのR用拡張でした。
ざーっと見た感じ、オムニ補完ができそうだし、古いバージョンのRでも補完ができるようになっていたりと、とても高機能っぽいです。
ちなみにオムニ補完とは、もともとVimの補完方法につけられている名称で(多分)、VSのインテリセンスのようなものだと思うとわかりやすいかもしれません。
これを使うと例えば、JavaやPythonのメソッドをさくっと補完していくことができるのですごく便利なのです。

で、早速使ってみたのですが、これがなかなか思うように動いてくれません。
まず、思いもよらないところでオムニ補完が働いたりするくせにオムニ補完してほしいところでしてくれません。これはオムニ補完のトリガーとなっている正規表現が変てこなことに起因するのですが、あまり肌に合いませんでした。
それと、オムニ補完時は、補完候補が読みこまれているオブジェクトに限定されるようになっているのですが、変なところでオムニ補完が働くせいで入力したいものが補完候補になかったりします。

というわけで結局, 肌にあったauto-complete-acr.elやねーということになったわけですが、せっかくなのでauto-complete-acr.elもオムニ補完ができるようにしてみました。
んですが書いた後に、そもそもCLOS風味なオブジェクトシステムのRでオムニ補完って必要なくね?という事に気がついたのですが、とりあえずアップデートしてGitHubに置いときます。

auto-complete-acr.el at master from myuhe’s auto-complete-acr.el – GitHub

インストールなぞは以前書いたこちらを参考にどうぞ。

ESS on Emacsで快適に補完できるauto-complete-acr.el0.2リリースしました

肝心のオムニ補完ですが、今のところデータフレームの列にアクセスする$演算子とS4クラスのスロットにアクセスする@演算子に対応しています。
動作している動画を見れば一目瞭然かと思います。動画ではデータフレーム”dat”で列やスロットを入力しています。

これまでデータフレーム等の要素にアクセスする際、例えばhoge$fugaのように入力する場合で考えると、hogeは補完できる可能性があったわけですが、fugaについてはほとんど手入力であったと思います。今回のオムニ補完では、$や@がトリガーとなってauto-completeするので、$より後ろのfugaについても補完候補から選択できるようになります。さらに、オムニ補完時にはトリガー($とか@)より前の文字列に応じて補完候補が絞られるので、効率的に候補選択することができるようになりました。
実は、もっとオムニ補完を働かせる場面ってあるかなーと思ってたのですが、意外にありませんでした。いやあ、CLOSって素晴らしいですね!!
ほかにオムニ補完をきかせたいシチュエーションがあったら、教えてもらえると嬉しいです。
最後になりましたが、補完候補を表示している際にいきなり補完候補が落ちる現象を確認しています。あまり気になるようでしたらac-source-omni-pre-list-essacrとac-source-omni-pre-class-essacrを情報源からはずしてみてください。

それでは Enjoy R and Emacs!!

Rでいちいちstr()とかhelp()とかしたくない人は、tooltip使ったらいいかもしれない

恐らく、Rを使うなかで、かなりお世話になる関数にstr()とかsummary()があります。こいつらは、オブジェクトをいろいろと要約してくれる関数で、オブジェクトに期待するデータ型が入ってるかとか、変な計算してないか、とか確認するときにひっじょーに便利な関数。
この他にもRには便利な関数がたくさんあります。いろいろ便利な関数を手足のように使いこなしてこそ、Rの本領が発揮されるのですが、ここのところ加齢のせいか、引数が全く覚えられなくなってきました。んでいちいちhelpを引くのですが、結構面倒くさい。今回はこいつらを解決してくれるtipsをご紹介。

何が嫌なのか

str()なりsummary()は便利なので多用しているのですが、関数が吐き出すオブジェクトを使うことはほとんどありません。つまり、使い捨てな関数なのです。なのでESSではiESSとかにちょろっと書いて確認したりするわけですが、コーディングのテンポが乱されてなんか嫌です。できれば打たずにどうにかしたいなーとかナマグサなことを考えてしまいます。
もう一つの問題は、かなり深刻。昔は暗記ものが大の得意でした。日本史の年号の暗記なんてお手の物。覚えることが負担になることなんてほとんどなかったのです。んが、最近はほんと覚えられなくなってきました。もしかしたら覚えようという気力が起きてないのかもしれませんが、しょっちゅう忘れてしまいます。
特にRの関数は、膨大な引数持ってたりします。ま、引数忘れてもhelp引けば良いのですが、help()引くのも面倒くさい。そう、何もかもが面倒くさい、ダメな人間なのです。

ESSの眠れる機能を呼び起こす

SVNなどを追っていくとよくわかるのですが、ESSはかなり開発のペースが早いです。ので、どんどん新しい機能が追加されていくのですが、誰もドキュメントを整備しないので、実装されても誰か使ってるかわからない機能が結構あったりします。そんな機能の一つにtooltip-show-at-pointがあります。これはESS5.4からサポートされているようなのですが、デフォルトではnilになっているので、多分ほとんどの人が使っていないと思われる眠れる機能。どんなものかちょっと使ってみたいという方は.emacsに次のように書いてみましょう。んで、ふつーにRで書いていくと、変ったところがすぐわかると思います。
[cpp]
(add-hook ‘ess-mode-hook
(lambda ()
(define-key ess-mode-map "(" ‘ess-r-args-auto-show)))
(setq ess-r-args-show-as ‘tooltip)
[/cpp]

この機能を使って、先ほどの面倒くさいことをどうにかします。

まずはgithubからess-R-object-tooltip.elを落としてきて、ロードパスが通ったところに置いておきます。
ess-R-object-tooltip.el at master from myuhe’s auto-complete-acr.el – GitHub
そして、.emacsに(require ‘ess-R-object-tooltip)って書いておきましょう。
んで、ESSでてきとーに書いて関数の上でC-c C-gとしてみましょう。
例えば、データフレームの上ですると次のようなツールチップが上がってきます。
dataframe
これがrnorm()だとこんな感じで引数の説明が上がってきます。
func
つまり、これまで手入力してた面倒くさいこと諸々がキーバインド一発でできちゃうわけですね。ちなみに、すでに書いたようにこの関数はESS5.4以降でないと動きませんので、バージョンが古い人はこの機会にあげておきましょう。

参考リンク

今回のelispはミネアポリス在住の名無しさんからいただきました。便利なelispありがとうございました。
R Object Tooltips in ESS « Blogistic Reflections
ほかにもこんなのあるよーとかあったら教えてください。

ggplot2使ってみた

ここんとこ、日記も書かず、来週の支部会に向けて黙々とRと格闘しておりました。んで、解析作業はようやく目処がつきました。ヤター3連休丸潰しで格闘した甲斐があったというものです。

でわ、何をしていたかというと、まずMCMCの計算でWinBUGSが全く動き出さず、ウンウンうなって結局測定値をゴニョゴニョしてMCMCにかけることに。今度は推定したパラメータをまたしてもゴニョゴニョして正しい値に戻すとゆーもーほんとに面倒くさい方法でやってたわけです。

その成果がこちら。
test
今回は、ggplot2使って作図してみました。いちお作図コードの一部。
[cpp]
library(ggplot2)
sp.totalNO <- 31
ggdat <- data.frame(out.spmat)
gp1 <- ggplot(ggdat, aes(x=seq(82,118,len=100)))
gp1 <- gp1 + geom_area(aes(y=ggdat[,1],colour= topo.colors(sp.totalNO)[31],fill= topo.colors(sp.totalNO)[31]))

sp1 <- ggdat[,1]
assign(paste("sp",2,sep=""),ggdat[,2] + get(paste("sp",2-1,sep="")))
gp2 <- geom_ribbon(aes(ymin=get(paste("sp",2-1,sep=""))
,ymax=get(paste("sp",2,sep=""))
,colour= topo.colors(sp.totalNO)[2]
,fill= topo.colors(sp.totalNO)[2]
))
gp1+gp2
[/cpp]

ggplot2はすごくシンプルな作りになっていて、ベースをggplot()で吐き出して、後はレイヤーオブジェクト(今回はgeom_ribbonとgeom_area)を足していけばいいとゆー非常にシンプルな作り。なんですけど、いろいろといぢるとなんか変です。例えば、ループの中で大量にレイヤーオブジェクトを吐き出してもすべて同じオブジェクトに上書きされてしまいました。うーん書き方が悪かったのかしらん。今回はとりあえず大量のリボンを直書きで書いて足していくっつー泥縄な書き方で手打ちに。出来ちゃえばいいし、スマートとか考えてる場合じゃないし。

とりあえず、山場は越したですな。けど、来週は天草でお泊まり調査で、その次は東京出張。うむーなかなかひまにならない。

ループを使わずに1から10までの総和を計算させるのが流行っているらしい

ループを使わずに1から10までの総和を表示するプログラム – Bug Catharsis
vallog: 流行っていると聞いて – ループを使わずに1から10までの総和を表示するプログラム
ループを使わずに1から10までの総和を表示するプログラム – 医者を志す妻を応援する夫の日記

ので、Rで書いてみます。
[cpp]
sum(c(1:10))
[1] 55
[/cpp]
あ…

これだと面白くないので再帰っぽく。
再帰はループではないそうなので。

[cpp]
func <- function(n) {
if (n <= 1) return(1)
else n + func(n-1)
}
func(10)
[1] 55
[/cpp]

おーいい感じです。
お次はPythonで。
[python]
sum(range(1,11))
55
[/python]

…pythonお前ってやつは。。。

追記

@syou6162曰く
Twitter / Yasuhisa Yoshida: @myuhe Reduce辺りも追加をw

また、マニアックな。
でも、こっちの方がRっぽい気がする。
[cpp]
Reduce("+",c(1:10))
[1] 11
[/cpp]

始めて知ったんですが、Reduceってば単なるRの関数なんですね。
これ使ってる人見たことないですが、統計でご飯食べている人等にとっては便利な関数なんでしょう。

[cpp]
> Reduce
function (f, x, init, right = FALSE, accumulate = FALSE)
{
mis <- missing(init)
len <- length(x)
if (len == 0L)
return(if (mis) NULL else init)
f <- match.fun(f)
if (!is.vector(x) || is.object(x))
x <- as.list(x)
ind <- seq_len(len)
if (mis) {
if (right) {
init <- x[[len]]
ind <- ind[-len]
}
else {
init <- x[[1L]]
ind <- ind[-1L]
}
}
if (!accumulate) {
if (right) {
for (i in rev(ind)) init <- f(x[[i]], init)
}
else {
for (i in ind) init <- f(init, x[[i]])
}
init
}
else {
len <- length(ind) + 1L
out <- vector("list", len)
if (mis) {
if (right) {
out[[len]] <- init
for (i in rev(ind)) {
init <- f(x[[i]], init)
out[[i]] <- init
}
}
else {
out[[1L]] <- init
for (i in ind) {
init <- f(init, x[[i]])
out[[i]] <- init
}
}
}
else {
if (right) {
out[[len]] <- init
for (i in rev(ind)) {
init <- f(x[[i]], init)
out[[i]] <- init
}
}
else {
for (i in ind) {
out[[i]] <- init
init <- f(init, x[[i]])
}
out[[len]] <- init
}
}
if (all(sapply(out, length) == 1L))
out <- unlist(out, recursive = FALSE)
out
}
}
<environment: namespace:base>
[/cpp]

Rとscheme

こないだの勉強会では、valvallowさんがやっていたschemeのライブコーディングを見て、
「ああ、多分僕はlispを書いてはいけない人間なんだ。」、と落伍者の烙印を押された気分に浸っていたわけですが、今日改めてRについて勉強していたらショッキングな事実が判明。元ネタはR-FAQ日本語訳(2章) – RjpWiki。以下引用。

R のデザインは, 2つの既存の言語に強く影響を受けてきました. それは, Becker, Chambers および Wilks の S (「S とは何ですか?」 (What is S?) 参照) および Sussman の Scheme です. その結果として生まれた R の言語は, S と見た目が非常に似ているのに対して, R の実装およびセマンティクスは Scheme から派生しています. さらに詳しくは, 「R と S の違いは何ですか?」 (What are the differences between R and S?) を参照してください.

えーーーーまぢですか。
S言語は当たり前として、ご先祖さまがschemeだったなんて・・・似てるところが全くわからん。

Rのtwitterクライアント使ってフォローしている人を分類してみた

RjpWikiを見てたら、twitteRなる文字が。
予想どおり、twitterクライアントとして動作するパッケージでした。

とりあえず使ってみる

というわけで、早速インストール。ちなみに、ubuntuの場合は、Rcurlのインストールにlibcurl4-gnutls-devがいるみたいなので先にapt-getでいれときましょー

[cpp]
install.packages("RCurl")
install.packages("twitteR")
[/cpp]

使うときに少しハマったのがproxyの設定。実はtwitteRではproxyの設定ができないのです。ただ、twitteRで設定を保存するinitSession()は、RcurlのgetCurlHandle()のラッパーなので、getCurlHandle()でproxyを設定すれば良いようです。
後は、使うだけ。

[cpp]
library(twitteR)
init <- initSession("username", "password")
##proxy経由の場合は次のように書くと良い。
init <- getCurlHandle(userpwd=paste("username","password", sep=":"),.opts=curlOptions(httpHeader=c(Expect=""),proxy="http://proxy.com:8080"))
publicTimeline(init)#パブリックタイムライン表示
userTimeline("myuhe", init)#任意のユーザのタイムライン表示
mentions(init)#リプライを表示
friendsTimeline(init)
Rtweets(init)#R関連のつぶやきを表示
[/cpp]

画面には、こんな感じで出てきます。

Scrs.png

分析っぽいことしてみる

こんだけだと、すんごく使い勝手の悪いクライアントでしかないので、Rっぽいことやってみます。
とりあえず、お手軽なとこでクラスタ分析。僕がフォローしている人を勝手ながら標本にさせてもらいまして、その人たちのフォローしている人数、フォローされている人数、発言数使ってこれまた、勝手に分類させてもらいます。
こんな感じで必要なデータをぶっこ抜きまして、おもむろにhclust()

[cpp]
install.packages("Cairo")
library(Cairo)

followlist <- userFriends("myuhe",init)
followvector=NULL
followervector=NULL
status = NULL
name = NULL

for (i in 1: length(followlist)){
if ( followlist[[i]] @ screenName != "wordpress" && followlist[[i]] @ screenName != "google" && followlist[[i]] @ screenName != "rtm"){
followvector[i] <- followlist[[i]] @ friendsCount
followervector[i] <- followlist[[i]] @ followersCount
status[i] <- followlist[[i]] @ statusesCount
name[i] <- followlist[[i]] @ screenName
}
else{
followvector[i] <- NA
followervector[i] <- NA
status[i] <- NA
}
}

mat <- cbind(followvector,followervector)
mat <- cbind(mat,status)
rownames(mat) <- name
mat <- mat[complete.cases(mat),]
hc <- hclust(dist(mat))
Cairo(1000, 1000, file="plot.png", type="png", bg="white")
plot(hc)
dev.off()
[/cpp]

んで、出てきたデンドロビウム・・・じゃなかった、デンドログラムがこれ。

plot.png

ん~有名な人はちょこっと外れている節があるか??ま、サンプルが66ですからね。もっと多いとおもしろいカモですね。
Rユーザの人は自分のtwitterアカウントで遊んでみてはいかがでしょうか??
Rユーザでない人もこれを機にRを使ってみると新しい発見があるかもしれないですね。

auto-complete.elを拡張するR用のelisp書いた

※20090608追記:以下に書いているelispはさらに拡張されバージョンアップしています。詳しくは、auto-complete.elを拡張するR用のelispをさらに拡張してみた。をご覧ください。

Rでコードを書くときは、もっぱらESS on emacs(meadow)で、すこぶる便利に使わせてもらっているのですが、唯一不満な点が補完でした。abbrev、dabbrevといちおあることはあるのですが、補完候補が何かコマンドを押さないと出てこないし、しかも補完候補がすごく見にくい。

実は、昔perl書くときに使っていたxyzzyはsiteinit.lにごにょごにょすると、この問題はあっさり解決してしまうわけで。

ミニバッファの補完候補リストをチップヘルプにする

この機能になれていたというのもあって、emacsの補完のしょぼさはかなり苦痛でした。といいながらもESSの便利さに負けてしまい、xyzzyはサブのエディタに降格してしまったわけですが・・・

んで、学会終わって余裕もできたので、たまにはemacsの環境整理でもしてみるかと、Rjpwiki見てたら谷村さんから次のようなコメントが・・・

ESSのメーリングリストにも投げたのですが、便利すぎるauto-complete.elのR用拡張って誰か書いていませんか。rubyとかc++用の拡張のみならず、octave用の拡張も書かれているのに、R用のものがまだありません。私家版ですでに書いている人がいれば、是非公開してみてください。auto-complete-modeが初耳の方はDemoを見てみてください。ESSのR-modeで使ってみましたが、便利すぎます。後はあらかじめキーワードなどをリストアップした拡張を作れば、さらに便利になるのですが。。。 — 谷村 2009-03-04 (水) 12:02:57

そもそもauto-complete.el自体を知らなかったので、調べてみるとなんかxyzzyの補完候補リストみたいなものを出してくれるんだそうな。しかも容易に拡張可能ですでに各種言語の拡張やoctaveの拡張まで出てるんだそうな。詳しくは動画デモを見てもらえるとその便利さがわかってもらえると思うです。後、以下のリンクも参考になります。

テキスト入力中に補完候補を自動的に表示してくれる auto-complete.el をリリースしました — ありえるえりあ
情報源による拡張が可能な auto-complete 0.1.0 をリリースしました — ありえるえりあ
auto-complete.el にやてふ拡張を作ってみた – がべこれログ – Garbage Collection Log

お~知らなかった。便利すぎるやないですか。しかし、octaveの拡張elispが出ていて、なぜRが出てないんだろう・・・

というわけで、すでにどなたか書かれているかもしれませんが、R用のauto-complete.el拡張elisp書いてみました。
ソースコードはこちらから。今回は初めてgithubにホスティングしてみました。
ちなみにスクリプトを書くにあたってauto-complete-octave.elを参考(というかほとんどまんま)をさせてもらいました。thanks Yen-Chin,Lee!!

使い方
1.auto-complete.elとauto-complete-acr.elをロードパスが通ったフォルダに起きます。

2..emacsに次の3行を追加します。

[cpp]
(require ‘auto-complete)
(global-auto-complete-mode t)
(require ‘auto-complete-acr)
[/cpp]

4.emacsでESSを起動します。
5.M-x auto-complete-modeします。

その他
・ESSモードとiESSモードのhookに引っかけています。

ソースコード見ていただくとわかるのですが、ただRの関数をリストにしてぶち込んでいるだけです。芸がなくてごめんなさい。
関数もnamespace:statsのものを取り込んだだけで、まだかなり不完全です。最初はNAMESPACEファイルを動的に読み込んでリスト化させようかと試みましたが、そこまでの芸当はできませんでした(当方lisp歴3日)。
手力でNAMESPACEファイルからリストを生成しています。NAMESPACEからlispのリスト生成に使ったpythonスクリプトも恥ずかしながらgithubに置いときます。

[python]
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# NAMESPACEにある関数名をlispで使えるように変換させます。
# NAMESPACEのexpoort()の括弧内にある関数名をコピーしてファイルを指定してから使ってください。

import re
mer = []
f = open("D:RscriptessNAMESPACE")
for line in f:
line = line.strip()
mer.append(line)
f.close()
out = "".join(mer)
out = re.sub(" ","",out)
out = re.sub(‘"’,"",out)
out = re.sub(",","" "",out)
out = re.sub("^""",’"’,out)
out = re.sub("$""",’"’,out)
[/python]

ということで、このコードはgithubにてhostして凄腕lisperの方のさらなるハックを期待したい・・・なる他力本願なことを考えています。
改良したelispを作られた方いらっしゃったら、是非githubのアカウントを取ってidを教えてください。collborateしてbrunchできるようにします。
まだまだ知名度の低いgithubですが、使い出すと面白いですのでガンガンbrunch、pushしてもらえると凄くうれしいです。

後、答えられるか甚だ不安ですが、動かんぞ!!とかこんな馬鹿な書き方あるか!!みたいなお叱りはコメントでお願いします。

もう一回やり直し

たった今、今回の共同研究者であらせられるODさんからBUGSコードに対するダメだしメールが入りまして今回のモデルは個体群動態モデルにすることになりました。

 うん、今日は徹夜確定です・・・

MCMCの樹海で遭難中

こないだ、サクッと通ったMCMC。
実は結果がどうしようもないことになってしまっていて、ほとんどが0周辺をウロンコロンしているわけで。

となると、間伐しても発生する樹木の個体数なんかには影響ないよーってことに。
いやこれは・・・・と頭をリフレッシュさせて樹種数に丸め込むと、なんとなく言い訳ができそう。

つまり、種の定着、発生と言った現象には、skyfactorとかシカが食べたり、なんていうことで説明できるんだけど、個体の量的変化となると、レンジが小さすぎてそんなの入れても意味ない、ってなことをwinbugsは言いたいのだろう。

、とここまで行き着いたのが、昨日の話。

で、今日はスライドを鋭意作っていたわけですが、作っている途中にあることに気がつきまして。

種のあるなしレベルまでデータの質を落とすんであれば、それを樹種別にやってもいけるんでない??

つまり、種が出たか出ないかとゆー非常に単純な話に強引に持ち込んでいるわけなので、それをベルヌーイ分布で表せばいいのでは、と気づいてしまったのでした。
ま、この気付きには、単純にガウス分布をベルヌーイに変えてやって、元データをちょこっといぢればたいして時間はかかるまい、なる予断があったわけですが、MCMCの神様はイヂワルで、trapの嵐。

先日のやつは

「初期値、定義しとらんやんけ!!」

と怒られても、

「しょうがない、gen.initでてきとーにつくったるわ」

って感じでサクサクいってたわけですが、今回は「もーだめ~~~」と言ってるのか、どうかはしりませんがしょっぱなでtrap。今度は根こそぎ初期値入れて回すと最初は回りだすものの、途中であえなくtrap・・・・・

今度は、初期値操作で再チャレンジ。うちのメイン妖精さんがウンウンうなっている間、サブ妖精さんでこの文章を書いてたりするわけです。だが異様に時間かかってる・・・・

は~終わりが見えない。

1857追記
やっぱりだめ。計算は終わるものの、チェーンうまく絡みません。む~これがイバさんの言っていたマルチモーダルというやつ??
できないものは仕方ないので収束の悪い変数を減らして再々々チャレンジ。

人力モデル選択・・・・

1921追記
変数を1にまで絞るものの、やはりだめ。この時点で万策つきました。
む~ベルヌーイが固すぎるのかも。ま、一番の問題はデータが少なすぎるということなんですけどね。ここまでくるとあきらめるしかあるまい。

あ~半日ほどの努力が報われず。