# -*- mode: R -*-
##EmacsとESS使いの方に。モードラインの設定。
##2014年8月1日(金曜日)

##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) {
  ##次の気分を決める。
  ##sample(x, size, replace = FALSE, prob = NULL)
  state <- sample(1:nState, 1, prob = tpm[state, ])###状態からの変遷
  ##probは、選択する確率。
  history <- c( state, history )
}
###各気分にどの位の確率で居るか
hist(history)

###気分遷移の履歴
plot(history, xlim=c(1,30), type="b")
title("history of transition")

###どの位続けて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 <- hist(lenhistory, breaks=range(lenhistory)[2]+1)
plot(a$counts, log="y")
title(sprintf("log(counts at state %d)",1))
##a$counts[-1]/a$counts[-length(a$counts)]
##par(new=T)
##plot(0.4^(0:20)*600,type="b")

###どの位続けて2の気分(ブルーな気分)にいるか
len <- 0
lenhistory <- c()
for( i in 1:length(history)){
  if(history[i]==2){
    len <- len+1;
  } else {
    if(len!=0){
      lenhistory <- c(len, lenhistory)
    }
    len <- 0;
  }
}
1;
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パッケージの準備
##install.packages("hmm.discnp", repo="http://cran.md.tsukuba.ac.jp", dependencies=TRUE)
library(hmm.discnp)

##連続9日間の表情
##sim.hmm(nsim,tpm,Rho,ispd=NULL,yval=NULL,verb=FALSE)
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

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