### 2010年 7月23日(金)

##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)

plot(rbind(maru,batsu),pch=c(rep("o",9), rep("x",8)))  

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

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")
        )
##X2でやってみる。
matplot(sort(mb[,2]),
        cbind(cumsum((mb[,3])[order(mb[,2])]),
              cumsum((mb[,4])[order(mb[,2])])),
        pch=c("o","x")
        )

##???はて、どこで区切れば良いのだろうか??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)
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])

#########################
######  HMMの説明　######
#########################

## 1:不安定状態。30%の確率で安定状態に変化する
## 2:安定状態。20%の確率で不安定状態に変化する

###大気状態の遷移関数。
###遷移行列。tpm (transition probability matrix)

tpm <- matrix(c(0.7, 0.3,
                0.2, 0.8),2,2,byrow=TRUE)
history <- c()
nState <- ncol(tpm) ### 状態の数
state <- 1
for (i in 1:10000) {
  state <- sample(1:nState, 1, prob = tpm[state, ])###状態からの変遷
  history <- c( state, history )
}
###各状態にどの位の確率で居るか
hist(history)
###revise.ispd(tpm)

###
sample(x, size, replace = FALSE, prob = NULL)
###配列xから、sizeを取り出す。配列の各要素が取り出される確率prob。
poi <- c()
for (i in 1:10000) {
  state <- sample(1:3, 1, prob=c(0.1,0.7,0.2))
  poi <- c( state, poi )
}
hist(poi)

###どの位続けて1の状態にいるか
len <- 0
lenhistory <- c()
for( i in 1:length(history)){
  if(history[i]==1){
    len <- len+1;
  } else {
    if(len!=0){
      lenhistory <- c(len, lenhistory)
    }
    len <- 0;
  }
}
  
hist(lenhistory, breaks=range(lenhistory)[2]+1)

##外から観測できるのは、晴れ1、曇り2、雨3。ラベル。
##ラベルの出現頻度は、大気の状態によって変わるはず。

Rho <- matrix(c(0.7, 0.2, 0.1,
                0.1, 0.4, 0.5), 3,2,byrow=FALSE)

安定状態::晴れ(0.7)、曇り(0.2)、雨(0.1)
不安定状態::晴れ(0.1)、曇り(0.4)、雨(0.5)

HMMパッケージの準備
[パッケージ]>>[ローカルにあるZIPファイルから]
hmm.discnp_0.0-5.zipを選択
[パッケージ]>>[パッケージの読み込み]
hmm.discnpを選択

##コマンドで実行するには、
utils:::menuInstallLocal()##が実行された模様
library(hmm.discnp)

y.sim <- sim.hmm(1000,tpm,Rho,1)
hist(y.sim)
y.sim
   [1] 3 2 2 3 3 1 1 3 1 1 2 2 3 3 3 2 2 2 1 2 2 3 2 2 3 2 2 2.....
  [47] 3 1 2 1 1 1 3 3 3 3 1 2 1 1 2 3 3 3 3 2 3 1 1 1 2 1 2 2.....
  [93] 3 1 1 1 1 2 2 1 1 3 2 3 2 3 3 1 1 2 1 2 3 3 1 2 2 2 1 1.....
 [139] 2 2 1 2 2 1 1 1 1 2 1 1 1 3 3 2 3 1 1 2 1 3 2 2 3 3 3 2.....
 [185] 2 2 2 1 2 2 2 1 1 1 1 3 2 3 3 3 2 3 2 3 3 2 1 3 2 1 3 2.....
 [231] 1 1 2 3 2 3 3 2 2 1 2 1 1 2 1 1 1 2 2 1 2 2 3 3 2 2 3 2.....
 [277] 3 2 3 2 3 3 1 2 3 3 3 1 1 2 1 3 2 3 3 3 3 3 1 1 3 2 2 2.....
 [323] 2 3 3 2 3 2 3 2 3 3 3 2 2 3 1 3 1 2 1 1 1 3 1 1 3 2 1 3.....

###100年分の夏休みの宿題より
y.sim <- sim.hmm(1000,tpm,Rho,100)
###MarkovModelのステートの数Kが2だとすると。
try <- hmm(y.sim,K=2,verb=T)
try
$Rho
           [,1]      [,2]
[1,] 0.09772003 0.6932688
[2,] 0.40110030 0.2039947
[3,] 0.50117967 0.1027365

$tpm
          [,1]      [,2]
[1,] 0.7970685 0.2029315
[2,] 0.2997133 0.7002867

$ispd
[1] 0.5962725 0.4037275

$log.like
[1] -107697.4

$converged
[1] TRUE

$nstep
[1] 83

$data.name
[1] "y.sim"

###MarkovModelのステートの数Kが3だとすると。
try <- hmm(y.sim,K=3,verb=T)
try
$Rho
           [,1]      [,2]       [,3]
[1,] 0.06276731 0.2084462 0.71718686
[2,] 0.39263076 0.3974432 0.19061852
[3,] 0.54460193 0.3941106 0.09219462

$tpm
          [,1]      [,2]      [,3]
[1,] 0.6981252 0.1178648 0.1840099
[2,] 0.1497446 0.6883309 0.1619245
[3,] 0.1881359 0.1237247 0.6881394

$ispd
[1] 0.3620603 0.2792944 0.3586452

$log.like
[1] -107696.3

$converged
[1] TRUE

$nstep
[1] 72

$data.name
[1] "y.sim"

try <- hmm(y.sim,K=5,verb=T) 
try
$Rho
           [,1]       [,2]      [,3]      [,4]       [,5]
[1,] 0.03984065 0.09612057 0.1974515 0.4586146 0.75464584
[2,] 0.39262263 0.40397032 0.3961624 0.3081072 0.17215400
[3,] 0.56753673 0.49990911 0.4063861 0.2332782 0.07320016

$tpm
           [,1]       [,2]       [,3]       [,4]      [,5]
[1,] 0.63486340 0.07153536 0.07028649 0.07749362 0.1458211
[2,] 0.07993295 0.63867604 0.07152319 0.07534211 0.1345257
[3,] 0.09086917 0.08197891 0.63023560 0.07355949 0.1233568
[4,] 0.10978490 0.09209232 0.07719730 0.61147925 0.1094462
[5,] 0.11604565 0.09128368 0.07180533 0.06322356 0.6576418

$ispd
[1] 0.2163845 0.1891563 0.1636488 0.1556082 0.2752022

$log.like
[1] -107699.0

$converged
[1] TRUE

$nstep
[1] 76

$data.name
[1] "y.sim"

###受理の確率を知りたい。。。
##try <- hmm(y.sim,par0=list(tpm=tpm, Rho=Rho),verb=T)
##hmmから取り出した。学習なし。
###rp$llcは、yの分すべて出てくる。上記の例では、1000x2個
###logとってのsumは、かけ算。

calloglike <- function (y, tpm, Rho)
{
  fy <- ffun(y, Rho)
  nc <- ncol(y)
  rp <- recurse(fy, tpm, nc, 0.000000001)
  ll <- matrix(rp$llc, ncol = nc)
  apply(log(abs(ll)), 2, sum)
}

calloglike(sim.hmm(1000,tpm,Rho,2), tpm, Rho)

tpm1 <- matrix(c(0.5, 0.5,
                0.5, 0.5),2,2,byrow=TRUE)
Rho1 <- matrix(c(0.3, 0.3, 0.4,
                0.3, 0.4, 0.3), 3,2,byrow=FALSE)
calloglike(sim.hmm(1000,tpm1,Rho1,2), tpm, Rho)

#####################################
######  HMMをつかったEIの予測　######
#####################################


##splice.All <- read.table("splice.dat",header=TRUE)
##spliceEI <- splice.All[splice.All$label=="EI",]

poi <- read.table("tmp.txt",header=F)
poi <- poi[,c(1,7:15,21)]
tpoi <- t(poi)

try <- hmm(tpoi,K=11,verb=T)
contour(try$tpm)
barplot(try$Rho)
barplot(as.vector(try$Rho[1:4,]))

ss <- mps(tpoi, try)
ss <- viterbi(tpoi, try)