## -*- mode: R -*-
## 2015年10月16日(金)
##EmacsとESS使いの方に、モードラインの設定。

##plot(sin)
##setwd(choose.dir())

##13Pageの図。
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)

##svmのマージン最大化の説明
plot(rbind(maru,batsu),type="n")
##pch=21は色ありの●。
points(maru[,],pch=21,col=1,cex=2,lw=3,bg="orange")
points(batsu[,],pch=4,col=2,cex=2,lw=3)

col2rgb("orange")/255
plot(rbind(maru,batsu),type="n")
##21は色ありの●。
for( i in 1:1000 ){
  points(maru[,]+5*rnorm(length(as.numeric(maru))),
         col=NULL, bg=rgb(1, 0.64, 0, alpha=0.1), pch=21, cex=2)
}
points(maru[,],pch=21,col=1,cex=2,bg="orange",lw=3)

for( i in 1:1000 ){
  points(batsu[,]+5*rnorm(length(as.numeric(batsu))),
         pch=4,col=1,cex=2)
}
points(batsu[,],pch=4,col=2,cex=2,lw=2)


##さて、X1軸だけの値を頼りに分割
##plot=Fにしておくと、ヒストグラムを書かない。
brks <- hist(rbind(maru,batsu)[,1],plot=F)$breaks
a <- hist(rbind(maru,batsu)[,1], breaks=brks,plot=F)
##type="n"にしておくと、枠だけ。データの線を書かない
plot(a$mids, a$density, ylim=c(0, 0.04),type="n")
a <- hist(maru[,1], breaks=brks, plot=F)
lines(a$mids, a$density, type="b",col="orange",lw=3)
par(new=T)
a <- hist(batsu[,1], breaks=brks, plot=F)
lines(a$mids, a$density, type="b", col=1, lw=3)

     
plot(rbind(maru,batsu),type="n")
points(maru[,],pch=21,col=1,cex=2,bg="orange",lw=3)
points(batsu[,],pch=4,col=1,cex=2,lw=2)


mb <- rbind(cbind(maru,1,0),cbind(batsu,0,1))
##X1の要素
mb[,1]
##X1の要素を、小さい順に並べ替えると
sort(mb[,1])
##sortしたときに何番目に来るか??
order(mb[,1])
##その順に並べてみると、「小さい順」
(mb[,1])[order(mb[,1])]
##その順に並べてみると、maruか??
(mb[,3])[order(mb[,1])]
##端から順に足してみる。
cumsum((mb[,3])[order(mb[,1])])
##cf::cumsum(1:10)
##maruの個数が閾値に従って、どの様に変化するか??
plot(sort(mb[,1]),cumsum((mb[,3])[order(mb[,1])]))
##batsuの個数が閾値に従って、どの様に変化するか??
plot(sort(mb[,1]),cumsum((mb[,4])[order(mb[,1])]))
##両方まとめて書いてみる
matplot(sort(mb[,1]),
        cbind(cumsum((mb[,3])[order(mb[,1])]),
              cumsum((mb[,4])[order(mb[,1])])),
        pch=c("o","x"),
        type="b"
        )
##X2でやってみる。
matplot(sort(mb[,2]),
        cbind(cumsum((mb[,3])[order(mb[,2])]),
              cumsum((mb[,4])[order(mb[,2])])),
        pch=c("o","x"),
        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")

tm <-3
tb <-2

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

m <- 9-tm
b <- 8-tb
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()
## install.packages("tree_1.0-25.zip", repo=NULL)
## plot(sin)
##setwd(choose.dir())
##utils:::menuInstallLocal()
## library(tree)

############################
######  決定木の説明  ######
############################
##setwd(choose.dir())
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
##localにある、「tree_1.0-36.zip」を読み込んで使う。
##utils:::menuInstallLocal()。これじゃダメ。
library(tree)
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(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(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(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(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(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:2){
  a.mdl <- prune.tree(a.mdl, best=ts) ##枝刈りして、全体大きさがtsの物を作る。
  ##どの位出来ているかテスト。
  splice.prd <- predict(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(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:2, 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(splice.mdl, newdata=splice.All[-lrn,], type="class")
table(splice.prd,(splice.All[-lrn,])$label)
plot(splice.mdl); text(splice.mdl)

splice.prd <- predict(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(dat.mdl, type="class"), dat[,3])

##おまかせ決定木
dat.mdl <- tree(dat[,3] ~ dat[,1]+dat[,2])
table(predict(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)

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