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

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

##sample(x, size, replace = FALSE, prob = NULL)
##配列xから、sizeを取り出す。配列の各要素が取り出される確率prob。

##barplot(c(0.1,0.7,0.2),col=c("yellow","blue","red"),names.arg=1:3)
##それぞれのボールの出現確率が、(0.1,0.7,0.2)
##「ボール取り出しの試行」を10000回繰り返した時の、出現頻度。

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:不安定状態。30%の確率で安定状態に変化する(green)
## 2:安定状態。20%の確率で不安定状態に変化する(red)

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

tpm <- matrix(c(0.7, 0.3,
                0.2, 0.8),2,2,byrow=TRUE)
tpm[1, ]
#[1] 0.7 0.3
##1は「安定状態」。
##次の状態は、0.7の確率で安定状態(1)に留まる、0.3の確率で不安定状態(2)に変化。

tpm[2, ]
#[1] 0.2 0.8
##2は「不安定状態」。
##次の状態は、0.2の確率で安定状態(1)に変化、0.8の確率で不安定状態(2)に留まる。

barplot(t(tpm),col=c("green","red"),names.arg=c("stable","unstable"))


history <- c()
nState <- ncol(tpm) ### 状態の数
state <- 1 ##最初は、安定からStart
for (i in 1:10000) {
  state <- sample(1:nState, 1, prob = tpm[state, ])###状態からの変遷
  ##probは、選択する確率。
  history <- c( state, history )
}
###各状態にどの位の確率で居るか
hist(history)


###どの位続けて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)
##a$counts[-1]/a$counts[-length(a$counts)]
##par(new=T)
##plot(0.4^(0:20)*600,type="b")
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パッケージの準備
install.packages("hmm.discnp", repo="http://cran.md.tsukuba.ac.jp", dependencies=TRUE)
library(hmm.discnp)

##連続9日間の天気
y.sim <- sim.hmm(rep(9,10),tpm,Rho,ispd=c(1,0))
y.sim
## [[1]]
## [1] 1 1 2 3 1 2 1 2 2
## [[2]]
## [1] 1 3 1 1 1 2 1 3 3
##.......
## [[10]]
## [1] 1 1 1 1 1 3 2 3 3


##「100年分の夏休みの宿題」データより推定
y.sim <- sim.hmm(rep(40,100),tpm,Rho,ispd=c(1,0))
##MarkovModelのステートの数Kが2だとすると。
try <- hmm(y.sim,K=2,verb=T,stationary=T)
try$Rho
try$tpm

barplot(t(rbind(tpm,try$tpm)),
        names.arg=t(outer(c("true", "estimated"),1:2,paste)),
        col=c("green","red"))

barplot(cbind(Rho,try$Rho),
        names.arg=t(outer(c("true", "estimated"),1:2,paste)))


1;

##MarkovModelのステートの数Kが4だとすると。
try <- hmm(y.sim,K=4,verb=T,stationary=T)
t(try$tpm)
try$Rho

##いかさまに、状態、2，3，4をまとめてみる。
barplot(t(try$tpm),
        names.arg=paste("estimated",1:4))
tmp <- t(try$tpm)
apply(t(try$tpm)[2:4,],2,sum)
tmp[2,] <- apply(t(try$tpm)[2:4,],2,sum)
tmp[3:4,] <- 0
tmp

tmp[,2] <- apply(tmp[,2:4],1,mean)
tmp[,3:4] <- 0
tmp
barplot(tmp[1:2,1:2],
        names.arg=paste("estimated",1:2))


barplot(cbind(t(tpm),tmp[1:2,1:2]),
        names.arg=t(outer(c("true", "estimated"),1:2,paste)),
        col=c("green","red"))
##なんと無く似ている??


##状態数4のラベルの出現頻度は、???。
barplot(try$Rho,names.arg=paste("estimated",1:4))

##いかさまに、状態、2，3，4をまとめてみる。
tmp <- try$Rho
tmp[,2] <- apply(tmp[,2:4],1,mean)
tmp[,3:4] <- 0
tmp

##ラベルの出現頻度は、割と似ている??
barplot(cbind(Rho,tmp[,1:2]),
        names.arg=t(outer(c("true", "estimated"),1:2,paste)))

barplot(c(0.1,0.4,0.5),col=c("yellow","blue","red"),names.arg=1:3)
barplot(c(0.8,0.1,0.1),col=c("yellow","blue","red"),names.arg=1:3)
barplot(c(0.3,0.3,0.4),col=c("yellow","blue","red"),names.arg=1:3)
        
##他のモデルに当てはめると、どうなるか??
##受理の確率を知りたい。。。
y.sim <- sim.hmm(rep(40,100),tpm,Rho,ispd=c(1,0))

s.vit <- viterbi(y.sim, tpm=tpm, Rho=Rho, ispd=c(1,0))
s.vit##状態のリスト。ispd=c(1,0)なので、すべての配列は状態1から始まる。
rownames(Rho) <- 1:3 ##no 'dimnames' attribute for arrayは、これだった。
mean(log(pr(s.vit, y.sim, tpm=tpm, Rho=Rho)))
## -43.22726

###
tpm.uso <- matrix(c(0.5, 0.5,
                0.5, 0.5),2,2,byrow=TRUE)
Rho.uso <- matrix(c(0.3, 0.3, 0.4,
                0.3, 0.4, 0.3), 3,2,byrow=FALSE)

s.vit <- viterbi(y.sim, tpm=tpm.uso, Rho=Rho.uso, ispd=c(1,0))
s.vit##状態のリスト。ispd=c(1,0)なので、すべての配列は状態1から始まる。
rownames(Rho.uso) <- 1:3 
mean(log(pr(s.vit, y.sim, tpm=tpm.uso, Rho=Rho.uso)))
## -24.34478
#####################################
######  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)