R for Statistical Computing (on Fedora Core 1)

(まだ書き中。まとまってません)

Rさわってみました。備忘録程度。FC1に依存してることがあるかも。セットアップからとりあえず使えるようになるとこまで。

00 スタート

基本はmathematica使いだけれども、高価すぎということで統計環境Rを使ってみることに。パタン識別や画像処理もできますよ(The R Bookにちょこっと載っています)ということで、関連分野で使うにはー？という調査。ばりばり統計なひとはS言語やS-Plusとかでもいろんな本が出ているのでそちらを参考にしたほうがよろし。

まったく知らない状態から始めたんですが、学習ソースとしては、(1)R入門を読む(邦訳)。(2)このtipsを読む(pdf版もあるし)。(3)RjpWikiにある例とかを睨む。でだいたい事足りるような気も。わりとシンプルかも。(ただしmatlabと同時に使うと見た目似ているので少々こんがらがる…)

◎

01 セットアップ

RとESS(Emacs Speaks Statistics)セットアップ(on Fedora Core 1)は以下の感じ(だいたい)。◎ 

% cd tmp
% wget http://cran.md.tsukuba.ac.jp/src/base/R-1.9.1.tgz
% tar xvzf R-1.9.1.tgz
% cd R-1.9.1
% ./configure
% make
% su
% make install
% cd ..
% wget http://stat.ethz.ch/ESS/downloads/ess-5.2.2.tar.gz  
% tar xvzf ess-5.2.2.tar.gz -C /usr/local/share/emacs/site-lisp/
% exit

基本システムはこれで使えます。ktermに「R」とだけタイプしてください。

つぎ、パッケージのセットアップ。拡張パッケージをいれとかなくては便利性に欠けるということで。(ただし、Rをシステムディレクトリに入れた場合、管理者権限が必要であるので以下は適宜suしてください)。パッケージはCRANに登録されているものなら、Rセッション中でinstall.packages('パッケージ名')ってやれば勝手にとってきてインストールしてくれる。たとえばklaRをインストールするには以下。

% su
% R
 :
 略
 :
> install.packages('klaR')
>

ネットワークが使えるPCなら、このようにすれば勝手にCRANにとりにいって勝手にインストールされます。この際、インストールされたパッケージについて以下のような情報が表示されます。 

> install.packages('klaR')
trying URL `http://cran.r-project.org/src/contrib/PACKAGES'
Content type `text/plain; charset=iso-8859-1' length 191426 bytes
opened URL
.......... .......... .......... .......... ..........
.......... .......... .......... .......... ..........
.......... .......... .......... .......... ..........
.......... .......... .......... ......
downloaded 186Kb

Error in unique(pkgs) : Object "klaR" not found
> install.packages('klaR')
trying URL `http://cran.r-project.org/src/contrib/PACKAGES'
Content type `text/plain; charset=iso-8859-1' length 191426 bytes
opened URL
.......... .......... .......... .......... ..........
.......... .......... .......... .......... ..........
.......... .......... .......... .......... ..........
.......... .......... .......... ......
downloaded 186Kb

trying URL `http://cran.r-project.org/src/contrib/klaR_0.3-3.tar.gz'
Content type `application/x-tar' length 68798 bytes
opened URL
.......... .......... .......... .......... ..........
.......... .......
downloaded 67Kb

* Installing *source* package 'klaR' ...
** R
** data
** help
 >>> Building/Updating help pages for package 'klaR'
     Formats: text html latex example 
  B3                                text    html    latex   example
  EDAM                              text    html    latex   example
  NaiveBayes                        text    html    latex   example
  TopoS                             text    html    latex   example
  b.scal                            text    html    latex   example
  betascale                         text    html    latex   example
  calc.trans                        text    html    latex   example
  centerlines                       text    html    latex   example
  countries                         text    html    latex   example
  distmirr                          text    html    latex
  dkernel                           text    html    latex   example
  dmvnorm                           text    html    latex
  drawparti                         text    html    latex
  e.scal                            text    html    latex   example
  errormatrix                       text    html    latex   example
  hmm.sop                           text    html    latex   example
  partimat                          text    html    latex   example
  plot.NaiveBayes                   text    html    latex   example
  predict.NaiveBayes                text    html    latex   example
  predict.rda                       text    html    latex   example
  predict.sknn                      text    html    latex   example
  predict.svmlight                  text    html    latex   example
  quadplot                          text    html    latex   example
  quadtrafo                         text    html    latex   example
  rda                               text    html    latex   example
  shardsplot                        text    html    latex   example
  sknn                              text    html    latex   example
  stepclass                         text    html    latex   example
  svmlight                          text    html    latex   example
  triframe                          text    html    latex   example
  trigrid                           text    html    latex   example
  triperplines                      text    html    latex   example
  triplot                           text    html    latex   example
  tripoints                         text    html    latex   example
  tritrafo                          text    html    latex   example
  ucpm                              text    html    latex   example
* DONE (klaR)

Delete downloaded files (y/N)? y

>

またパッケージ(ライブラリ)の説明が見たいときはRセッション中で以下のようにすればよいです。library()でインストールされているライブラリのリスト、引数に「help=名前」とすればそのパッケージの説明が出ます。library(名前)とすれば実際にロードされます。ライブラリの各関数にはたいていexampleがついていて、example(関数名)とすれば実行結果を見れます(グラフィクスが出るものはばかっと表示されるので注意)。詳しくは最初にあげたサイトなど見てください。以下は一例。 

> library()
> library(help=kernlab)
> library(kernlab)
> example(ksvm)
> q()

あえて、tar ballから入れたいひとは(もしくはCRANにないけどwww上で配布されているRパッケージのインストール)、foo.tar.gzを持ってきて、例えば

% wget http://stat.ethz.ch/~dettling/LogitBoost_1.1.tar.gz
% R CMD INSTALL LogitBoost_1.1.tar.gz

とやればインストールされます。以下、わたしの興味あった(インストールした)パッケージ群。以下の説明ではこれらのインストールを仮定してるかもしれません。例えば、library(ほにゃ)として「ほにゃ」を呼ぶ場合は「ほにゃ」がインストールされている必要があります。

パッケージはCRANのパッケージリストにたくさんあるのでディスクに余裕があるならいろいろ試してみると良いかもしれません。僕が使う範囲で入れてみたパッケージは以下の感じです。他にもいろんなパッケージがあります。 日本語での説明がRjpwikiにありました。

パッケージ名説明
quadprogFunctions to solve Quadratic Programming Problems
kernlabKernel Methods Lab
klaRClassification and visualization
mlbenchMachine Learning Benchmark Problem
e1071Misc Functions of the Department of Statistics (e1071), TU Wien
bootBootstrap R (S-Plus) Functions (Canty)
ade4Analysis of Environmental Data
fieldsTools for spatial data
gregmiscGreg's Miscellaneous Functions
randomForestBreiman's random forest for classification and regression
treeClassification and regression trees
ipredImproved Predictors
gbmGeneralized Boosted Regression Models
fastICAFastICA algorithms to perform ICA and Projection Pursuit
classPPProjection Pursuit for supervised classification
knncatNearest-neighbor classification with categorical variables
knnTreek-nn classification with variable selection inside leaves of a tree
mdaMixture and flexible discriminant analysis
snaTools for Social Network Analysis

◎

02 ESSでR

ktermとかからつかっても良いですが、とりあえずEmacs内で動かしてます。基本的には初心者なんで以下は参考程度。

~/.emacs.elに以下を足して、M-x Rで起動します。 

(load "/usr/local/share/emacs/site-lisp/ess-5.2.2/lisp/ess-site.el")

ただし、パスは上のインストールみたく展開した場合で、展開したところにパス通せばよいだけ。 もしessにloadパスが通ってたら((setq load-path (cons "/path/to/ess/lisp-directory" load-path)する手もあり)、

(require 'ess-site)

だけで動くと思われます。以下によくつかうキーバインドを表にしてみました。でも、ESSのページにリファレンスカードがPDFでおいてありました。

ESS(.Rファイルの編集)
C-c M-r選択範囲を評価しiESSに飛ぶ
C-c M-bバッファを評価しiESSに飛ぶ
C-c M-jカーソル行を評価しiESSに飛ぶ
C-c M-f関数を評価しiESSに飛ぶ
C-c C-r選択範囲を評価
C-c C-bバッファを評価
C-c C-jカーソル行を評価
C-M-x関数を評価
C-M \範囲をインデント
C-M a関数の頭へ移動
C-M e関数の末尾へ移動
C-c C-lファイルをRにロード
iESS(Rとのインタラクティブセッション)
M-r/M-pまえのコマンド
C-cエラー箇所に飛ぶ
C-c C-lファイルをロード
C-c C-qRセッションを終了
C-c C-xオブジェクトのリスト表示
C-c C-vオブジェクトのヘルプ表示(「?名前」と同じ)

FC1だと.emacs.elに例えばてきとうにessのフックをかけて画面をわっておくとちょうどemacsを最大化してつかえば統合環境のように見えます。ESSでのRの起動は少々時間はかかります。helpが小窓に出ますがバッファの移動具合によっては以下のままだと変になることもあるかもしれません(ESSのせい)

アドホックなので、自分でわっても良いんですが(C-x 2とC-x 3で)、以下のコードを~/.emacs.elに加えておくと、「M-x R」でR起動したときに勝手に画面わって配置します(たぶん)。環境によってはスプリットする順番かえなきゃいかんかもしれんですけど。

;; ess
(load "/usr/local/share/emacs/site-lisp/ess-5.2.2/lisp/ess-site.el")
(defun ess:format-window-1 ()
  (split-window-horizontally)
  (other-window 1)
  (split-window)
  (other-window 2))
(add-hook 'ess-pre-run-hook 'ess:format-window-1)

あとはバッファに関数オブジェクトつくって、それをRセッションに送って評価(関数単位の評価、選択範囲の評価、バッファ全体の評価、行単位の評価、などいろいろでけます)しつつ、書き直します。またセッションのほうに移動して関数の利用手続きなどインタプリタ窓として使います。

いろいろ変えたい場合、/usr/local/share/emacs/site-lisp/ess-5.2.2/lisp/の.elを読むとだいたいわかります。ユーザが変えられる変数(メニューを出すかどうかとかいろいろ)はess-cust.elのUser changeable variablesて書いてあるところ以下に書いてあります。.emacs.elで読んでsetqすればよさげです。

変数はess-cust.elみてもらうとして、たとえば、ESSはワーキングディレクトリを聞いて来ますが、うぜー、毎回ーおれは同じディレクトリしかまず使わねーって場合とか、あとでsetwdしますがなーという場合とかは~/.emacs.elなどに

(setq ess-ask-for-ess-directory nil)
(setq ess-directory "/home/1gac/R/workspace/")

などと対応する変数をnilにしとけば良いです。ess-directoryをsetしない場合はemacsを起動したディレクトリになります。あとでsetwdする手もありです。

◎

03 gnuplotの代用

linuxだとデータをグラフのプロットにするのにgnuplot使うひとが多いと思うのですが、まあ、Rで代用したほうが便利なこともあるかもしれんし、ということでそこから。こんなデータファイル"hoge.dat"があるとします。

-1.343287 1.555834
1.873737 1.053686
-0.8624064 -0.5302879
-0.5959827 1.566078
-0.5541462 -0.4270813
1.253616 -0.5671337
-0.6297169 0.02841587
0.7199952 -0.8629212
-0.8054724 0.299316
0.5773122 1.687105

ちなみにこれはRだと以下のようにして作った標準正規乱数(rnormのnormは正規分布の乱数だということを表しています)の2x10行列です。xとやったときと同じにするためには一応t(x)と転置してwriteする必要があります。

> x=cbind(rnorm(10),rnorm(10))
> write(t(x),"hoge.dat",ncolumns=2)
> x

これをいろいろなグラフとして表示し、psやpdfとして出力してみるなんてことがgnuplotのもっとも多い使われ方(少なくともわたしの周辺では)ではないかと思うので、それをRでやってみますと。Rではデータは次元数×サンプル数の行列として扱われるみたいなのですが、各列(属性)と各行(サンプル番号)に名前つけて、その名前で扱ったりできるようにデータフレームつう形式も用意されてます。上のファイルをデータフレームで読みます。 

> y = read.table("hoge.dat")
> y
           V1          V2
1  -1.3432870  1.55583400
2   1.8737370  1.05368600
3  -0.8624064 -0.53028790
4  -0.5959827  1.56607800
5  -0.5541462 -0.42708130
6   1.2536160 -0.56713370
7  -0.6297169  0.02841587
8   0.7199952 -0.86292120
9  -0.8054724  0.29931600
10  0.5773122  1.68710500
> y[1,1]
[1] -1.343287
> y[1,]
         V1       V2
1 -1.343287 1.555834
> y$V1
 [1] -1.3432870  1.8737370 -0.8624064 -0.5959827 -0.5541462  1.2536160
 [7] -0.6297169  0.7199952 -0.8054724  0.5773122
> colnames(y) = c("x1","x2")
> y
           x1          x2
1  -1.3432870  1.55583400
2   1.8737370  1.05368600
3  -0.8624064 -0.53028790
4  -0.5959827  1.56607800
5  -0.5541462 -0.42708130
6   1.2536160 -0.56713370
7  -0.6297169  0.02841587
8   0.7199952 -0.86292120
9  -0.8054724  0.29931600
> y[3:7,]
          x1          x2
3 -0.8624064 -0.53028790
4 -0.5959827  1.56607800
5 -0.5541462 -0.42708130
6  1.2536160 -0.56713370
7 -0.6297169  0.02841587

上のように添字などでデータフレーム(列名・行名つきデータ行列)を調べたりできます。これをいろいろプロットしてみます。y$x1は一列目の数値列(1次元データ)です。read.tableはsepとしてセパレータを指定できるのでデリミタがスペース以外でも利用できます。「?read.table」でヘルプ見てください。

hoge
二次元データの散布図

plot(y)
色オプションと折れ線グラフ

plot(y$x1,col=rgb(1,0,0),type='l')
色オプションと折れ線グラフ+点プロット(type='b'だと丸が刺さらない)

plot(y$x1,col=rgb(1,0,0),type='o')
色オプションと折れ線グラフ+点プロット

plot(y$x1,col=rgb(1,0,0),type='l') points(y$x1,col=rgb(0,0,1))
plot(y$x1,col=rgb(1,0,0),type='h')
points(y$x1,col=rgb(0,0,1))
二種類のグラフの同時プロット

plot(y$x1, xlab="x", ylab="y", type='o',lty=1)
lines(y$x2,type='o',lty=2)
二種類のグラフの同時プロット(凡例の8,2は座標)

matplot(y,type='l')
matpoints(y,pch=1:2)
legend(8,2,legend=c("hoge1","hoge2"),
col=c,lwd=1,pch=1:2)
棒グラフ

barplot(abs(y$x1))
円グラフ

pie(abs(y$x1))

まあRのグラフィクスはいろいろできるのでもしもうちょい凝ったことをしたければ最初にあげたサイトなどをみてみてください。プロットした絵を出力したい場合は以下のようにするのが一番てっとりばやいです。

何かplot処理をして絵が出ている状態で
> dev.print(device=png,file="foo.png",width=400,height=400)
> dev.print(device=postscript,file="foo.ps")
> dev.print(device=pdf,file="foo.pdf")
など

◎

04 プログラミング

はじめにあげたところなどを見るとコツなどがつかめますが、いちおう、行列とベクトルやnnet(ニューラルネット)あるいはJavaとのインタフェースなどについてはこちらにちょこっと書きました。

ESS使わない場合、プログラムはfoo.Rに書いて、Rセッション中で以下のように実行する。

> source('foo.R')

時間かかりそうな数値実験とかバッチモードで動かす。R プログラムファイル myfile.R のバッチ実行(結果は myfile.Rout に記録される)は以下の感じ。

% R CMD BATCH myfile.R &

(行列・データフレーム。data.frame。factor。latex組版とか、いろいろ。未完)

> options(digits=23)
> .Machine$integer.max   # 最大可能正浮動小数 (=2^(31))
> .Machine$double.xmax   # 最大可能正整数
> system("命令名")
> getwd()
> setwd("directory")
> system("ls -F")

◎

05 自前パッケージ

とりあえず、他言語で書いた関数を外から呼ぶ方法はRjpwikiに分かり易いC言語の場合の簡単な例があります。重い処理などはCで書いたほうが良いんかもー。

あと、パッケージにするにはいろいろRの基準に従う必要があるらしく、これは"Writing R Extensions"という詳しいドキュメントがあります(和訳もあります。必要ならRjpwikiとかからたどってください)。

比較的新しいRでは最適化関数optimなど様々な関数がRのAPIを通してC言語のソースから呼べるようです。詳しくは上記ドキュメント参照。

◎

06 パタン識別(1)

パタン識別を実装する。mlbenchはUCIリポジトリなどにある有名なデータなど含んでいるため、これを使ってみる。まず必要なライブラリをロードして、使えるデータのリストを見る(ライブラリはinstall.packageなどでインストールされている必要がある)。

library(kernlab)
library(mlbench)
data()                 # 使えるデータ表示

「data()」で使えるデータの一覧と簡単な説明が表示される。また「data(package='mlbench')」などとすればパッケージ内のデータの一覧と簡単な説明を表示する。実際にロードするには次のように「data(データ名)」とする。

ここではSonarデータをロードしてpairsプロットで始めのほうの特徴だけを見てみることにする。

######### データロード　＆　データの詳細調査  ###############
data(Sonar)            # Sonarデータを(データフレームとして)ロード
dim(as.matrix(Sonar))  # データ行列のサイズ表示
summary(Sonar)         # サマリ表示
colnames(Sonar)        # 属性表示
unique(Sonar$Class)    # Classラベルを見る
# はじめの5属性(5属性)だけのサブセットを選択
pairs(Sonar[1:10], pch = 21,bg = c("red", "green3")[unclass(Sonar$Class)])

以下のようなはじめの5特徴についてのペアプロットが出る。

次にこのSonarデータからSupport Vector Machineを用いて識別器を構成する。ここではとりあえず半分を訓練データとし、残りの半分をテストデータとする。

######### ロードしたデータをSVMで識別してみる  ###############
# テストデータと訓練データを作成
num    = dim(Sonar)[1]     # データのサンプル数取得
num.tr = as.integer(num/2) # 半分を訓練に
train  = Sonar[1:num.tr,]
test   = Sonar[(num.tr+1):1,]
# library(kernlab)のksvmを訓練(識別器を作成)
# "Class"というのはSonarデータのクラスラベル属性についている名前
svm.classifier = ksvm(Class~.,data=train,kernel="rbfdot",kpar=list(sigma=1))
# 情報を見る
svm.classifier
# つくった識別器でテストデータをテスト
result = predict(svm.classifier,test)
# 結果を表示
paste("RESULT: correct ",length(result[result==test$Class]),"/",dim(test)[1])

各々のステップで以下のように情報が出力される。


> svm.classifier


Support Vector Machine object of class "ksvm" 

SV type: C-classification 
 parameter : cost C = 1 

Gaussian Radial Basis kernel function. 
 Hyperparameter : sigma =  1 

Number of Support Vectors : 104 
Training error : 0 


> cat(paste("RESULT: correct ",length(result[result==test$Class]),"/",dim(test)[1],"\n"))


RESULT: correct  104 / 105

サンプル数も多くないので、半分訓練、半分テスト、では識別率の推定が雑であるため、もっと凝った推定を使う。ここではipredパッケージを用いて10-fold cross validationと632+ bootstrapを利用する。

######### もっと詳細な方法を用いて誤識別率を推定 ######### 
# (ipredを使う)
# 今回使うksvmやlibrary(e1071)のsvmなどはerrorestの規格に沿っているけど
# 一応、自前の関数を渡すさいの参考のためにモデルと予測方式をあえて指定
mymodel = function(formula,data){
  ksvm(formula,data=data,kernel="rbfdot",kpar=list(sigma=0.1))
}
mypredict = function(object, newdata){
  predict(object,newdata)
}
# 10-fold cross validation
library(ipred)
est1 = errorest(Class~., data=Sonar, model=mymodel, predict=mypredict,estimator="cv", est.para=control.errorest(k = 10))
est1
# 632+ bootstrap estimation (100 bootstrap) 
est2 = errorest(Class~., data=Sonar, model=mymodel, predict=mypredict,estimator="632plus", est.para=control.errorest(nboot=100))
est2
# 結果を表示
str1 = paste("ERROR RATE: ",est1$error,"(by 10-fold cv)")
str2 = paste("ERROR RATE: ",est1$error,"(by 632+ bootstrap)")
cat(paste(str1,"\n",str2,"\n",sep=""))

以下のように簡単に誤式別率が得られる。この例では100サンプルのブートストラップと10-fold cvはまあ同じ精度で推定されている。


> cat(paste(str1,"\n",str2,"\n",sep=""))

ERROR RATE:  0.230769230769231 (by 10-fold cv) 
ERROR RATE:  0.230769230769231 (by 632+ bootstrap)

◎

07 パタン識別(2)

今度は問題を低次元(2次元)にすることで識別境界そのものを描画してみる。ここではふたつの正規分布からの正規乱数を用いてサンプルを生成し、それに対するsvmとknnの識別境界をプロットしてみる。以下のような図が得られる。

###################################
## 設定
library(e1071)    # e1071をロード
par(mfrow=c(1,3)) # 絵を三枚出すことにする(なくてもいい)
resol = 100       # グラフィクス解像度
n1    = 100       # サンプル数(クラス1)
n2    = 100       # サンプル数(クラス2)
# 各クラスの平均と分散行列
mean.c1 = c(1,1) 
var.c1  = matrix(c(1,.5,2,.5),nrow=2)
mean.c2 = c(-1,-1) 
var.c2  = .5*diag(2)
# まず、二次元の正規乱数を行列で作成する関数を定義
rnorm2d = function(num,mean=c(0,0),variance=diag(2)){
  var = svd(variance);
  U = var$u;  V = var$v;
  D = diag(sqrt(var$d));
  B = U %*% D %*% t(V);
  w = c();
  for(i in 1:num) w = append(w,mean + B %*% rnorm(2));
  matrix(w,ncol=2);
}
# データ作成
cat("Generating data ...\n")
cls1 = rnorm2d(n1,m=mean.c1,v=var.c1)
cls2 = rnorm2d(n2,m=mean.c2,v=var.c2)
x = rbind(cls1,cls2)
y = c(rep(1,nrow(cls1)),rep(-1,nrow(cls2)))
# plot範囲の計算
rng = apply(x,2,range)
matplot(rng[,1], rng[,2], type= "n", xlab = "feat.1", ylab = "feat.2")
points(cls1,pch=21,bg="black")
points(cls2,pch=21)
###################################
## svmモデル作成(訓練)
cat("Training ...\n")
model = svm(x,y)
# テスト点をグリッド上に作成
tx = seq(rng[1,1],rng[2,1],length=resol)
ty = seq(rng[1,2],rng[2,1],length=resol)
pnts = matrix(nrow=length(tx)*length(ty),ncol=2)
k = 1
for(j in 1:length(ty)){
 for(i in 1:length(tx)){
    pnts[k,] = c(tx[i],ty[j]);
    k = k+1;
  }
}
# svmモデル予測(識別)
cat("Classifying ...\n")
re = predict(model,pnts)
# データのプロット
z = matrix(re,nrow=length(tx),ncol=length(ty))
image(tx,ty,z,xlab="feat.1",ylab="feat.2")
contour(tx,ty,z,add=TRUE, drawlabels=FALSE, level=0)
title(main = "SVM classification", font.main = 4)
points(cls1,pch=21,bg="black")
points(cls2,pch=21)
box()
###################################
## 比較のためk-最近隣法のプロットを作成
k = 1
nn = function (a,b, k=1) {
  mean( y[ order( (x[,1]-a)^2 + (x[,2]-b)^2 )[1:k] ] )
}
res = matrix(re,nrow=length(tx),ncol=length(ty))
for(i in 1:length(tx))
  for(j in 1:length(ty)){
    res[i,j] <- nn(tx[i],ty[j])
  }
contour(tx,ty,res,drawlabels=FALSE, level=0)
title(main = paste(k,"-Nearest Neighbor Plot"), font.main = 4)
points(cls1,pch=21,bg="black")
points(cls2,pch=21)
points(matrix(tx, nr=length(tx), nc=length(tx)),
       matrix(ty,nr=length(ty),nc=length(ty),byrow=T),
       pch=".",col=c("red", "blue")[ as.numeric( res >0 )+1])
abline(h=rng[1,1]:rng[2,1],v=rng[1,2]:rng[2,2],lty=2)
## 出力
dev.print(device=png,file="test2.png",width=900,height=300)

◎