# -*- mode: R -*-
##EmacsとESS使いの方に。モードラインの設定。
## 2011年 7月 8日(金)
##
##plot(1:10)
##setwd(choose.dir())

##「空間を2分割」を組み合わせ、の図。
maru <- matrix(
               c(13, 72,
                 11, 35,
                 37, 57,
                 23, 41,
                 26, 21,
                 40, 68,
                 43, 50,
                 68, 75,
                 89, 77),
               ncol=2 ,byrow=T)
batsu <- matrix(
                c(68, 50,
                  77, 43,
                  88, 46,
                  62, 27,
                  78, 29,
                  62, 31,
                  55, 13,
                  89, 16),
                ncol=2 ,byrow=T)

cbind(maru,1,0)
cbind(batsu,0,1)

mb <- rbind(cbind(maru,1,0),cbind(batsu,0,1))

plot(rbind(maru,batsu),
     pch=c(rep("o",9), rep("x",8)),
     col=c(rep("black",9), rep("red",8)))
##maruとbatsuを表示する。repは前にやった物。繰り返してくれる。
##pchは、マークとして使う"文字"

##全データの「X1の要素」
mb[,1]
##「X1の要素」を、小さい順に並べ替えると
sort(mb[,1])
##「X1の要素」をsortしたときに、各データは何番目?
##2番目の"11"が最小値なので、1
order(mb[,1])
##元のデータ(mb[,1])を「order()の順」に。つまりsortと同じ
(mb[,1])[order(mb[,1])]
##「maruかbatsuか」を「order()の順」に
(mb[,3])[order(mb[,1])]
##あたまから順に、n番目まで足してみる。
cumsum((mb[,3])[order(mb[,1])])
cumsum(1:10)
##cumulate【1】積み重ねる,積み上げる,集積する;累加する.

##maruの個数が閾値に従って、どの様に変化するか??
plot(sort(mb[,1]),cumsum((mb[,3])[order(mb[,1])]),type="b")
##batsuの個数が閾値に従って、どの様に変化するか??
plot(sort(mb[,1]),cumsum((mb[,4])[order(mb[,1])]),type="b")
##両方まとめて書いてみる
matplot(sort(mb[,1]),
        cbind(cumsum((mb[,3])[order(mb[,1])]),
              cumsum((mb[,4])[order(mb[,1])])),
        pch=c("o","x"), col=1:2
        ,type="b")
##X2でやってみる。
matplot(sort(mb[,2]),
        cbind(cumsum((mb[,3])[order(mb[,2])]),
              cumsum((mb[,4])[order(mb[,2])])),
        pch=c("o","x"), col=1:2
        ,type="b")
        )

##???はて、どこで区切れば良いのだろうか??X1軸、X2軸??
##天下り的に情報量
ent <- function(x){ if(x==0.0){ 0.0 } else {  -x*log(x) }}
plot(seq(0,1,0.1), 
     sapply(seq(0,1,0.1), ent)+
     sapply(1.0-seq(0,1,0.1), ent), type="b")

##例えば、例えば、「有る閾値」を決めた。
##それよりか「大きいか小さいか?」を教えて貰えたとする。。。

##1)閾値よりも小さい場合
m <-3 ##○の数
b <-2 ##×の数

pm <- m/(m+b) ##○割合。確率。
pb <- b/(m+b) ##×割合。確率。
ans <-(m+b)*(ent(pm)+ent(pb)) ##得られる情報(増える情報量の期待値)

##2)閾値よりも大きい場合
m <- 9-m
b <- 8-b

pm <- m/(m+b)
pb <- b/(m+b)
(m+b)*(ent(pm)+ent(pb))

##で、もし閾値で判断した時の増える情報量。
ans + (m+b)*(ent(pm)+ent(pb))

poi <- function(x){##(○の数,×の数)
  m <-x[1]
  b <-x[2]
  if(m+b==0){
    ans <- 0;
  } else {
    pm <- m/(m+b)
    pb <- b/(m+b)
    ans <- (m+b)*(ent(pm)+ent(pb))
  }
  
  m <- 9-x[1]
  b <- 8-x[2]
  if(m+b==0){
    ans;
  } else {
    pm <- m/(m+b)
    pb <- b/(m+b)
    ans + (m+b)*(ent(pm)+ent(pb))
  }
}
poi(c(3,2))

z <- cbind(cumsum((mb[,3])[order(mb[,1])]),
           cumsum((mb[,4])[order(mb[,1])]))

plot(sort(mb[,1]),apply(z, 1, poi),type="b")

z <- cbind(cumsum((mb[,3])[order(mb[,2])]),
           cumsum((mb[,4])[order(mb[,2])]))

plot(sort(mb[,2]),apply(z, 1, poi),type="b")
     
  
############################
######  Treeパッケージの準備
############################
## [パッケージ]>>[ローカルにあるZIPファイルから]
## ***.zipを選択
## [パッケージ]>>[パッケージの読み込み]
## treeを選択
##コマンドラインで実行するには、utils:::menuInstallLocal()
##これで、SelectFilesで、ローカルにあるZIPを指定し、インストールする。
##install.packages("tree_1.0-28.zip", repos="http://cran.md.tsukuba.ac.jp")
##最新のを入れる。
install.packages("tree", repos="http://cran.md.tsukuba.ac.jp")
##使う前に、以下を呼んでから
library(tree)

############################
######  決定木の説明　######
############################
splice.All <- read.table("splice.dat",header=TRUE)

splice <- splice.All[splice.All$label=="EI"|splice.All$label=="N",]
summary(splice)
##label    s6    s7    s8     s9     s10    s11   s12   s13    s14     
##EI: 764  A:675 A:846 A: 473 A: 377 A: 425 A:809 A:945 A: 433 A:528   
##IE:   0  C:702 C:533 C: 440 C: 403 C: 413 C:432 C:491 C: 472 C:569   
##N :1651  G:572 G:522 G:1054 G:1223 G: 400 G:772 G:542 G:1049 G:570   
##         T:466 T:514 T: 448 T: 412 T:1177 T:402 T:437 T: 461 T:748   

splice.mdl <- tree(label ~ ., data=splice)

plot(splice.mdl)
text(splice.mdl)

splice.mdl

##node), split, n, deviance, yval, (yprob) * denotes terminal node
##1) root 2415 3014.00 N ( 0.316356 0.000000 0.683644 )  
##  2) s10: A,C,G 1238   86.42 N ( 0.005654 0.000000 0.994346 ) *
##  3) s10: T 1177 1534.00 EI ( 0.643161 0.000000 0.356839 )  
##    6) s9: A,C,T 322    0.00 N ( 0.000000 0.000000 1.000000 ) *
##    7) s9: G 855  608.90 EI ( 0.885380 0.000000 0.114620 )  
##     14) s13: A,C,T 190  260.40 EI ( 0.563158 0.000000 0.436842 )  
##       28) s8: A,C,T 63   24.12 N ( 0.047619 0.000000 0.952381 ) *
##       29) s8: G 127  120.20 EI ( 0.818898 0.000000 0.181102 ) *
##     15) s13: G 665  143.40 EI ( 0.977444 0.000000 0.022556 )  
##       30) s12: A 488   14.38 EI ( 0.997951 0.000000 0.002049 ) *
##       31) s12: C,G,T 177   97.90 EI ( 0.920904 0.000000 0.079096 )  
##         62) s8: A,C,T 23   31.84 N ( 0.478261 0.000000 0.521739 ) *
##         63) s8: G 154   21.35 EI ( 0.987013 0.000000 0.012987 ) *

##予測と実際の違い。
splice.prd <- predict.tree(splice.mdl, newdata=splice, type="class")

##予測と本当のラベル
table(splice.prd,splice$label)

#############################
##全部のラベルでやってみよう。
#############################
splice <- splice.All
summary(splice)
splice.mdl <- tree(label ~ ., data=splice)
plot(splice.mdl)
text(splice.mdl)
splice.mdl

##予測
splice.prd <- predict.tree(splice.mdl, newdata=splice, type="class")

##予測と本当のラベル
table(splice.prd,splice$label)

##これでは、細かくならない。
splice.mdl <- tree(label ~ ., data=splice, mincut=1, minsize=2)
##これでも、細かくならない。
splice.mdl <- tree(label ~ ., data=splice, mindev=0.00000000001)
##これで、細かく出来る。
splice.mdl <- tree(label ~ ., data=splice, mincut=1, minsize=2, mindev=0)
splice.mdl <- tree(label ~ ., data=splice, mincut=1, mindev=0)

##これで、細かく出来る。
splice.mdl <- tree(label ~ ., data=splice, minsize=1, mindev=0)
splice.prd <- predict.tree(splice.mdl, newdata=splice, type="class")

plot(splice.mdl,type="uni")
text(splice.mdl)

##予測と本当のラベル
table(splice.prd,splice$label)


###些細なルールを作ってしまうと
##7割のデータを選んで、学習。3割のデータで、テスト。
lrn <- sort(sample(1:length(splice.All[,1]),length(splice.All[,1])*0.7))

splice <- splice.All[lrn,]
##summary(splice)
##splice.mdl <- tree(label ~ ., data=splice)
splice.mdl <- tree(label ~ ., data=splice, minsize=1, mindev=0)
##type="uni"で枝の長さを揃える
plot(splice.mdl, type="uni"); text(splice.mdl)

##3割のデータでテスト
splice.prd <- predict.tree(splice.mdl, newdata=splice.All[-lrn,], type="class")
##予測と本当のラベル
table(splice.prd,(splice.All[-lrn,])$label)

##ちょっとだけ、[-lrn]の使い方。
seq(10,100,by=10)
##5番目、3番目、8番目を選ぶ。
(seq(10,100,by=10))[c(5,3,8)]
##5番目、3番目、8番目を含まない物全部。
(seq(10,100,by=10))[-c(5,3,8)]
##つまり、splice.All[lrn,]番目を学習に用いて、splice.All[-lrn,]を検査に使う

a <- table(splice.prd,(splice.All[-lrn,])$label)
##対角線要素以外の総和
row(a)
col(a)
##対角線要素の抽出
a[row(a)==col(a)]
##対角線要素以外の抽出
a[row(a)!=col(a)]
##エラーの割合は、
sum(a[row(a)!=col(a)])/sum(a)


splice.mdl <- tree(label ~ ., data=splice, minsize=1, mindev=0)
plot(splice.mdl); text(splice.mdl)
##splice.mdl
##残りでテスト
splice.prd <- predict.tree(splice.mdl, newdata=splice.All[-lrn,], type="class")
##予測と本当のラベル
table(splice.prd,(splice.All[-lrn,])$label)


##植木屋さん。Trivialな枝から、切っていく。
lrn <- sort(sample(1:length(splice.All[,1]),length(splice.All[,1])*0.7))
splice <- splice.All[lrn,]
splice.mdl <- tree(label ~ ., data=splice, minsize=10, mindev=0.001)
summary(splice.mdl)

a <- table(splice.prd,(splice.All[-lrn,])$label)
##エラーの割合は、
sum(a[row(a)!=col(a)])/sum(a)


##んじゃ、どのくらいのサイズが良いのだろうか、虱潰し
splice <- splice.All
a.mdl <- tree(label ~ ., data=splice, minsize=1, mindev=0)
errOpen  <- c()
errClose <- c()
for(ts in 50:1){
  a.mdl <- prune.tree(a.mdl, best=ts) ##枝刈りして、全体大きさがtsの物を作る。
  ##どの位出来ているかテスト。
  splice.prd <- predict.tree(a.mdl, newdata=splice[-lrn,], type="class")
  a <- table(splice.prd,(splice[-lrn,])$label)
  ##対角線要素以外の総和
  errOpen <- c(errOpen,sum(a[row(a)!=col(a)])/sum(a))

  splice.prd <- predict.tree(a.mdl, newdata=splice[lrn,], type="class")
  a <- table(splice.prd,(splice[lrn,])$label)
  ##対角線要素以外の総和
  errClose <- c(errClose,sum(a[row(a)!=col(a)])/sum(a))
}
##1:Open:未学習データ。2:Close:学習データ。
matplot(50:1, cbind(errOpen,errClose), type="b")


##「どのくらいの大きさの木が良いのだろうか」をお任せで
#####cv.tree(splice.mdl,FUN=prune.tree)
lrn <- sort(sample(1:length(splice.All[,1]),length(splice.All[,1])*0.7))
splice <- splice.All[lrn,]
splice.mdl <- tree(label ~ ., data=splice, minsize=10, mindev=0.001)

splice.cv <- cv.tree(splice.mdl,FUN=prune.tree)

for(i in 2:5){
  splice.cv$dev <- splice.cv$dev + cv.tree(splice.mdl,, prune.tree)$dev
}
splice.cv$dev <- splice.cv$dev/5

plot(splice.cv)


identify(splice.cv$size, splice.cv$dev,splice.cv$size)
##右クリックで終了

##サイズが12が良さそう。#で木のサイズを12にしてみる
lrn <- sort(sample(1:length(splice.All[,1]),length(splice.All[,1])*0.7))
splice <- splice.All[lrn,]
splice.mdl <- tree(label ~ ., data=splice, minsize=10, mindev=0.001)
##これで、枝狩り。。
splice.mdl <- prune.tree(splice.mdl,best=12)
splice.prd <- predict.tree(splice.mdl, newdata=splice.All[-lrn,], type="class")
table(splice.prd,(splice.All[-lrn,])$label)
plot(splice.mdl); text(splice.mdl)

splice.prd <- predict.tree(splice.mdl, newdata=splice.All[-lrn,], type="class")
table(splice.prd,(splice.All[-lrn,])$label)


##判別分析との比較の為に
dat1 <- read.table("dat1.txt")
dat2 <- read.table("dat2.txt")

####dat2 <- dat2+c(-2,0)
dat2[,1]<- dat2[,1]-3

dat <- cbind(rbind(dat1,dat2),rep(c("red","blue"),c(1000,1000)))

##非常に細かい決定木
dat.mdl <- tree(dat[,3] ~ dat[,1]+dat[,2], minsize=1, mindev=0)

plot(dat.mdl)
text(dat.mdl)

plot(dat[,1], dat[,2], type="n")
text(dat[,1], dat[,2], c("b", "r")[dat[,3]], col=c("blue", "red")[dat[,3]])
partition.tree(dat.mdl, add = TRUE, cex = 1.5)

table(predict.tree(dat.mdl, type="class"), dat[,3])

##おまかせ決定木
dat.mdl <- tree(dat[,3] ~ dat[,1]+dat[,2])
table(predict.tree(dat.mdl, type="class"), dat[,3])
plot(dat[,1], dat[,2], type="n")
text(dat[,1], dat[,2], c("b", "r")[dat[,3]], col=c("blue", "red")[dat[,3]])
partition.tree(dat.mdl, add = TRUE, cex = 1.5)

##判別分析
dat.lda <- lda(dat[,3] ~ dat[,1]+dat[,2])
table(predict(dat.lda, dat)$class,dat[,3])