Rとデータマイニング入門
Rは5年以上前に、ちょこっと触ったのだけど
今や完全に忘却の彼方だったので、復習も兼ねて纏めてみた。
schemeの影響を受けてる為か関数型言語使いには優しい(^^;
(無限があったり個人的にはHaskellっぽいかなと。。)
R言語の特徴
- データマイニングに利用されるOSS言語
- 関連するアルゴリズム実装が豊富
- 動的型付け
- 高階関数
- 無限OK
- ベクトル（行列）処理に特化
- 行列計算は意外と速い
- 説明変数が多くても大丈夫（複雑なモデルでもOK）
- そこそこのデータ規模で威力を発揮
- TBクラスのデータは扱い切れない
などなど、、
TIPS
コーディング
普通はRを起動するのでしょうが、僕はEmacs & Rscriptでコーディングします。
a.r
#!/usr/bin/evn Rscript message("Hellow world !")
起動
bashからでも良いですが、emacsから'compile で起動すると楽かと。。
$ ./a.r
結果
./a.r Hellow world !
Emacs（R-mode）
ここからESSを拾ってきてR-modeを導入すると幸せになれるかも！
(add-to-list 'load-path "~/.emacs.d/ess")
(setq auto-mode-alist
(cons (cons "\\.r$" 'R-mode) auto-mode-alist))
(autoload 'R-mode "ess-site" "Emacs Speaks Statistics mode" t)解凍してlisp以下を~/.emacs.d/essへコピーして~/.emacsに書くだけです。
Emacs使いなら何も問題ないかと。。
elcが入ってたので環境によっては消さないと駄目かも？
find ~/.emacs.d/ess -name '*.elc' | xargs -n 1 rm
ベクトル１
単純なベクトルを作る
c
#!/usr/bin/evn Rscript A<-c(1,2,3,4,5) print(A) print(length(A))
こっちでも可
A<-1:5 print(A) print(length(A))
結果
[1] 1 2 3 4 5 [1] 5
ちなみにprintは省略できる。
ラベリング
名前を付けておいた方が後から何かと楽になる。
A<-1:5 names(A)<-c("aa","bb","cc","dd","ee") print(A)
結果
aa bb cc dd ee 1 2 3 4 5
ベクトル２
2x5のベクトルを作って色々
matrix
B<-matrix(0:9,2,5) B[1,2]<-99 message("== dim(B) ==") dim(B) message("== B ==") B message("== B[1,] ==") B[1,] message("== B[,3] ==") B[,3] colnames(B)<-c('aa','bb','cc','dd','ee') rownames(B)<-c('XX','YY') message("== B(labeled) ==") B
結果
== dim(B) ==
[1] 2 5
== B ==
[,1] [,2] [,3] [,4] [,5]
[1,] 0 99 4 6 8
[2,] 1 3 5 7 9
== B[1,] ==
[1] 0 99 4 6 8
== B[,3] ==
[1] 4 5
== B(labeled) ==
aa bb cc dd ee
XX 0 99 4 6 8
YY 1 3 5 7 9特徴的なのはこの辺りでしょう
- B[1,]
- 1番目のrowを示す
- B[,3]
- 3番目のcolumnを示す
- B[1,2]
- 上記２つが合わさる
ベクトルを纏めて扱う
Array , List , data.frame などがあるが省略！
ファイル操作
詳細は略！
- write
- ベクトルの状態で書き出す(read.tableで読む)
- write.table
- データフレームの状態で書き出す(read.table or scan で読む)
- save
- マーシャル書き出し(loadで読む)
詳しくはこちら
Rはこの辺り癖がある。
しかしファイル（＆文字列）操作を頑張るなら他に幾らでも良い言語があるのでそちらでやるべき（と思う）
行列計算
詳しい話は面倒なので省略ｗ
大学の講義でも思い出しながら・・・
元行列
message("== B ==") B<-matrix(0:35,5,5) colnames(B)<-c('C1','C2','C3','C4','C5') rownames(B)<-c('R1','R2','R3','R4','R5') B[1,2]<-99 B
== B == C1 C2 C3 C4 C5 R1 0 99 10 15 20 R2 1 6 11 16 21 R3 2 7 12 17 22 R4 3 8 13 18 23 R5 4 9 14 19 24
足し算
message("== B+matrix(100:129,5,5) ==") B+matrix(100:129,5,5)
== B+matrix(100:129,5,5) ==
C1 C2 C3 C4 C5
R1 100 204 120 130 140
R2 102 112 122 132 142
R3 104 114 124 134 144
R4 106 116 126 136 146
R5 108 118 128 138 148スカラー値を掛ける
message("== B*2 ==") B*2
== B*2 == C1 C2 C3 C4 C5 R1 0 198 20 30 40 R2 2 12 22 32 42 R3 4 14 24 34 44 R4 6 16 26 36 46 R5 8 18 28 38 48
行列を掛ける１
message("== matrix(c(10,20,30,40,50),5,1)*c(1,2,3,4,5) ==") matrix(c(10,20,30,40,50),5,1)*c(1,2,3,4,5)
== matrix(c(10,20,30,40,50),5,1)*c(1,2,3,4,5) ==
[,1]
[1,] 10
[2,] 40
[3,] 90
[4,] 160
[5,] 250行列を掛ける２
message("== B*c(1,2,3,4,5) ==") B*c(1,2,3,4,5)
== B*c(1,2,3,4,5) == C1 C2 C3 C4 C5 R1 0 99 10 15 20 R2 2 12 22 32 42 R3 6 21 36 51 66 R4 12 32 52 72 92 R5 20 45 70 95 120
集約
Row毎に集計
message("== SUM col ==") apply(B,1,sum)
== SUM col == R1 R2 R3 R4 R5 144 55 60 65 70
Row毎に集計（平均）
message("== MEAN col ==") apply(B,1,mean)
== MEAN col == R1 R2 R3 R4 R5 28.8 11.0 12.0 13.0 14.0
column毎に集計
message("== SUM row ==") apply(B,2,sum)
== SUM row == C1 C2 C3 C4 C5 10 129 60 85 110
平均値との差
message("== MY row ==") my<-function(x){ x-mean(x) } apply(B,2,my)
== MY row == C1 C2 C3 C4 C5 R1 -2 73.2 -2 -2 -2 R2 -1 -19.8 -1 -1 -1 R3 0 -18.8 0 0 0 R4 1 -17.8 1 1 1 R5 2 -16.8 2 2 2
グラフ表示
Rには組み込み済みのデータセットが用意されている
- iris
- アヤメの３品種(setosa,versicolor,virginica)と茎や花の計測値１５０件
data(iris) iris
Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa ： 50 5.0 3.3 1.4 0.2 setosa 51 7.0 3.2 4.7 1.4 versicolor 52 6.4 3.2 4.5 1.5 versicolor ： 100 5.7 2.8 4.1 1.3 versicolor 101 6.3 3.3 6.0 2.5 virginica 102 5.8 2.7 5.1 1.9 virginica ： 150 5.9 3.0 5.1 1.8 virginica
これらのデータをプロットしてみる。
単純にプロット
plot(iris[,1],iris[,3])
値をプロット
plot(iris[,1],iris[,3],type='n') text(iris[,1],iris[,3])
種で色分け
plot(iris[,1],iris[,3],pch = 21, cex=0.5,bg = c(2, 3, 4)[unclass(iris$Species)],col=c(2,3,4)[unclass(iris$Species)])
立体的にプロット
library(lattice) cloud(Sepal.Length ~ Petal.Length * Petal.Width, data = iris,groups=Species,pch=c(16,17,18),cex=0.5,bg = c(2, 3, 4),col=c(2,3,4)) cloud(Sepal.Length ~ Petal.Length * Petal.Width, data = iris,groups=Species,pch=c(16,17,18),cex=0.5,bg = c(2, 3, 4),col=c(2,3,4),screen=list(z=0,x=10,y=-15))
解析
次回に続く
（続）Rとデータマイニング入門
前回（グラフプロット迄）の続き
データ解析は複雑なデータを扱うので視覚情報に頼る事が多い。
(模様的にいい感じ！といった具合)
Rはグラフ出力が優秀で助かる！！
データマイニング
rpart(決定木)
150件全部のデータで学習させ、決定木を作ってみる。
グラフを見ると色々解る
- "setosa"は単純にPetal.Lengthで分離できる。
- "versicolor","virginica"はPetal.Widthで大体分離できるが少しは混ざる(エラー)
rpart.r
#!/usr/bin/evn Rscript library(rpart.plot) data(iris) # Learn and create tree tree <- rpart(Species~.,data=iris,method="class") # Dump tree tree # Plots rpart.plot(tree,yesno=F,extra=101,type=4,faclen=15,varlen=15) rpart.plot(tree,yesno=F,extra=101,type=4,faclen=15,varlen=15,fallen.leaves=T) rpart.plot(tree,branch.type=4,yesno=F,faclen=15,varlen=15,extra=101)
単純な木だからDumpでもソコソコ理解できる
n= 150
node), split, n, loss, yval, (yprob)
* denotes terminal node
1) root 150 100 setosa (0.33333333 0.33333333 0.33333333)
2) Petal.Length< 2.45 50 0 setosa (1.00000000 0.00000000 0.00000000) *
3) Petal.Length>=2.45 100 50 versicolor (0.00000000 0.50000000 0.50000000)
6) Petal.Width< 1.75 54 5 versicolor (0.00000000 0.90740741 0.09259259) *
7) Petal.Width>=1.75 46 1 virginica (0.00000000 0.02173913 0.97826087) *
予測
木を作れば、未知のデータに対しても予測ができる。
predict.r
150件のデータ中半分(75件)をサンプリングして学習し残りのデータを予測する。
#!/usr/bin/evn Rscript library(rpart) library(rpart.plot) data(iris) d<-iris # Sampling s <- sample(nrow(d),nrow(d)*0.5) d.sample <- d[s,] d.test <- d[-s,] # Create tree tree <- rpart(Species~.,data=d.sample,method="class") tree rpart.plot(tree,yesno=F,extra=101,type=4,faclen=15,varlen=15,fallen.leaves=T) # Predict predtree <- predict(tree,d.test,type="class") # Results table(d.test$Species,predtree)
横軸が『正解』縦軸が『予測』を表している。
つまり本当は"virginica"なのに"versicolor"と予測してしまったモノが３件ある
上のTreeの結果の通り"versicolor"と"virginica"の分離は難しい模様・・・
n= 75
node), split, n, loss, yval, (yprob)
* denotes terminal node
1) root 75 46 versicolor (0.28000000 0.38666667 0.33333333)
2) Petal.Length< 2.45 21 0 setosa (1.00000000 0.00000000 0.00000000) *
3) Petal.Length>=2.45 54 25 versicolor (0.00000000 0.53703704 0.46296296)
6) Petal.Width< 1.65 28 1 versicolor (0.00000000 0.96428571 0.03571429) *
7) Petal.Width>=1.65 26 2 virginica (0.00000000 0.07692308 0.92307692) *
predtree
setosa versicolor virginica
setosa 29 0 0
versicolor 0 21 0
virginica 0 3 22
random forest(決定木)
やっと今回の目的まで来た！！
そう。大体のデータはまずrandom forestにかけるベシ！
random forestでは上の様な決定木を幾つも作り、多角的に評価する。
上の単純な決定木では拾いきれない連中を助ける事ができる。（ハズ）
randomForest.r
#!/usr/bin/evn Rscript library(randomForest) data(iris) d<-iris s <- sample(nrow(d),nrow(d)*0.5) d.sample <- d[s,] d.test <- d[-s,] forest <- randomForest(Species~.,data=d.sample,ntree=500,importance=T) forest plot(forest) varImpPlot(forest) importance(forest) predforest <- predict(forest,d.test,type="class") table(d.test$Species,predforest)
あれれ？あまり変わんなぁ
3/75 = 4%
学習のサンプリング数が少ないし、random forest君もサンプリングのエラー率が6.67%と言ってる。
まあ予測精度ってのはこんなモノである。
randomForest 4.6-6
Type rfNews() to see new features/changes/bug fixes.
Call:
randomForest(formula = Species ~ ., data = d.sample, ntree = 500, importance = T)
Type of random forest: classification
Number of trees: 500
No. of variables tried at each split: 2
OOB estimate of error rate: 6.67%
Confusion matrix:
setosa versicolor virginica class.error
setosa 26 0 0 0.00000000
versicolor 0 22 3 0.12000000
virginica 0 2 22 0.08333333
setosa versicolor virginica MeanDecreaseAccuracy
Sepal.Length 2.929630 0.06475861 3.383496 1.8000081
Sepal.Width 1.140845 0.13359321 1.432476 0.8392769
Petal.Length 5.436318 5.38457538 6.160392 3.4376186
Petal.Width 5.171155 5.61849111 6.673574 3.5327468
MeanDecreaseGini
Sepal.Length 5.264504
Sepal.Width 1.676426
Petal.Length 20.868405
Petal.Width 21.441419
predforest
setosa versicolor virginica
setosa 24 0 0
versicolor 0 25 0
virginica 0 3 23少し解説するとMeanDecreaseGiniは変数の重要度。
つまりコノ部分はこう言っている
『幾ら多角的に評価してもPetal.LengthとPetal.Widthが圧倒的に重要である』
rpartと結果が変わらないのも納得
MeanDecreaseGini Sepal.Length 5.264504 Sepal.Width 1.676426 Petal.Length 20.868405 Petal.Width 21.441419
クラスタリング
サンプリングで学習するのでは無く、データの集合から同じ特徴を持つグループを見つけ出す手法。
アルゴリズムは色々あるのだが、難しい話は抜きにして色々試してみた。
それでも詳しい話が気になる人向けに・・・
- http://www.kamishima.net/jp/clustering/
- http://www.slideshare.net/hamadakoichi/webr-r
- http://www1.doshisha.ac.jp/~mjin/R/28/28.html
- http://www1.doshisha.ac.jp/~mjin/R/29/29.html
Original data
機械学習：Random Forest
全データでサンプリングして、同じデータをテストしているので、完全に特徴を掴めている。
=== RandomForest ===
predforest
setosa versicolor virginica
setosa 50 0 0
versicolor 0 50 0
virginica 0 0 50
機械学習：Rpart
やはり"versicolor" "virginica" の部分がとり切れない。
データプロットで見ると致し方ない事がよく解る。
=== Rpart ===
predtree
setosa versicolor virginica
setosa 50 0 0
versicolor 0 49 1
virginica 0 5 45
ソフト：Cmeans (fuzzy k-means)
ソフトクラスタリングは複数のクラスタに所属する事ができる。（比重付）
比重が60%を下回っているモノは(11,12,13)としてtable出力した。(setosaは完璧に分離出来たので13は無い)
プロットでは比重は色分けによって表現した。(RGB)
境目の状態が良く解る
=== Cmeans ===
c
1 2 3 11 12
setosa 0 0 50 0 0
versicolor 42 2 0 5 1
virginica 9 35 0 4 2
ソフト：Cmeans (extends UFC)
UFCの拡張らしい。詳しい事は良く解らない！！
cmeansとの違いはなんだろう？
=== Ufcl ===
c
1 2 3 11 13
setosa 0 50 0 0 0
versicolor 1 0 42 2 5
virginica 35 0 11 1 3
k-means
hartigan-wongが一番良い結果が出るらしいのだが・・・
実際の答えとは外れているが、プロット的には綺麗だ。
きっちり特徴を捉えているとも言える。
=== Kmeans hartigan-wong ===
1 2 3
setosa 50 0 0
versicolor 0 48 2
virginica 0 14 36
階層型：average
階層型クラスタリング
次元が大きく（変数が多く）なったり、乱れが大きいデータでは役にたたない。
ユークリッド距離を平方距離とした場合(下)：
なるほど！確かにそこで分けても良いよな！！
=== average ===
1 2 3
setosa 50 0 0
versicolor 0 50 0
virginica 0 14 36
1 2 3
setosa 50 0 0
versicolor 0 50 0
virginica 0 38 12
階層型：complete
確かにソコで分けても良さそう！
=== complete ===
1 2 3
setosa 50 0 0
versicolor 0 23 27
virginica 0 49 1
1 2 3
setosa 50 0 0
versicolor 0 23 27
virginica 0 49 1
階層型：median
あまり綺麗に分かれてない・・・
=== median ===
1 2 3
setosa 50 0 0
versicolor 0 14 36
virginica 0 49 1
階層型：ward
これも微妙だ・・・
=== ward ===
1 2 3
setosa 50 0 0
versicolor 0 50 0
virginica 0 14 36
1 2 3
setosa 50 0 0
versicolor 0 49 1
virginica 0 15 35
