RPyで遊んでみる

一学期の授業にてデータ解析という授業を取っています。統計ソフトRを用いて*1、車市場の線形価格予想モデルのようなものを構築します…というのはとりあえずおいておいてw。

Rを使っていろいろやるわけですが、RをPythonから使えるRPyというライブラリがあるのでつかってみました。

まず、RPyを使うのに必要なNumericというライブラリをインストールします。しなくてもいいのかもしれないけど、Strongly Recommendされているので入れておきます。ここからどうぞ。次に、RPyのインストール。ここからどうぞ。

これで使えるようになってると思うので、使ってみましょう。デモがよくできているので、この付近でいろいろ試してみましょう。

普通はこれでうまく行くと思うんですが、R側にあるデータや、functionを使う時にちょっとつまった点があったので、自己メモ。

R側にあるデータは「r.load("~/.RData")」のような感じでロードしてから、「r("hib")」のようにして取得できます。「""」がないとRではなく、Pythonのオブジェクトだと思われてしまうので注意。

あと、Rのfunctionで「function.plot()」のようなものがあるとして、RPyでは「r.funcion_plot()」のように使います。「.」を「_」に置き換えて考えましょう。講義で作ったような関数とデータを使うには以下のような感じでできました。

from rpy import *
r.load("~/.RData")

r.png("a.png")
r.scatterplot_bordered_by_boxplots(r("hib"),r("age"))
r.dev_off()

r.png("b.png")
r.plot_strip_n_medians(r("hib"),r("age"),5)
r.dev_off()

r.png("c.png")
r.plot_strip_n_boxplots(r("hib"),r("age"),5)
r.dev_off()

グラフもしっかりとできた。



簡単すぎるのでもう少し複雑（？）なことをやってみる。線形モデルでリグレッションしたときのベータハットを計算する関数と散布図の組合せを右上の部分だけ書くプログラム*2。pngの名前も使用した変数を使って分かりやすくとかしてみた。

from rpy import *
r.load("~/.RData")

def plot(x,y):
    r.png("./plot/"+x+"&"+y+".png")
    r.plot(r(x),r(y),xlab=x,ylab=y)
    r.dev_off()

def beta_hat(a,b):
    a=r(a)
    b=r(b)
    r.assign("y",r.t(r.rbind(b)))
    r.assign("X",r.cbind(a))
    return r("solve(t(X)%*%X,t(X)%*%y)")[0][0]
array=["weight","length1","price","mpg"]
print beta_hat("weight","price")
for i in range(len(array)):
    for j in range(len(array)):
        if i!=j and i>j : plot(array[i],array[j])

とにかくはまったのが、r.assignってところ。http://rpy.sourceforge.net/rpy/doc/rpy_html/Caveat-and-bugs.html#Caveat-and-bugsを見てみると分かるのだが、「a bit trickier」とか書いてあるんだもん。はまりますorz

*1:http://www.okada.jp.org/RWiki/

*2:pairsみたいな関数使えば問題ないのだが、ここでは知らないことにしとこっとw。