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

##plot(1:10)
##setwd(choose.dir())
##参照ファイルは、
## forPrint6.pdf
## command-6.txt##これ本体。
## ・read.table
##mojiDat <- as.matrix(read.table("t10k-images.txt"))/255
##mojiLab <- unlist(read.table("t10k-labels.txt"))
## ・install
##ZIPファイルには、注意すべし。
##install.packages("hmm.discnp_0.2-4.zip", repos=NULL) 
## ・source
##source("snnMoji.txt")
############


#########################
######  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は、1,2,3のいずれか。probは、それぞれが出現する確率。
  ##前回、「何が出たか」には関係しない。毎回独立事象。
  state <- sample(1:3, 1, prob=c(0.1,0.7,0.2))
  poi <- c( state, poi )##state
}
hist(poi, label=T)

##何故、この様に歯抜けなの???
hist(poi)$counts
##  [1] 1008 0 0 0 0 0 0 0 0 7007 0 0 0 0 0 0 0 0 0 1985
##って、なっている。この0の大群は、、
hist(poi)$breaks
str(hist(poi)$breaks)##21個
str(hist(poi)$counts)##20個
##breaksは、BINの境目を示している。countsの1つめが、1008なのは、
##「1.0以上で、1.1未満」のデータの個数はいくつか?ってこと。
1:4-0.5##新しい境目(breaks)を作ってみると、
hist(poi, breaks=1:4-0.5, label=T)
hist(poi, breaks=1:4-0.5, plot=F)$counts
##plot=Fは、ヒストグラムを描画しない。



## 1:いい気分(happy)。30%の確率でブルーな気分(blue)に変化する
## 2:ブルーな気分(blue)。20%の確率でいい気分(happy)に変化する

###気分の遷移関数。
###遷移行列。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","blue"),names.arg=c("happy","blue"))
title("transition probability matrix")

##おまけ、、、
##RのGraphichDeviceの所の、
##History>Recordingをチェックしておくと、
##PageUpとPageDownで昔のグラフが見られる。
##こうやっても、出来る。
graphics.off() 
windows(record = TRUE)
.SavedPlots <- NULL ##ってすると、Historyがクリアされる。らしい。

##BarPlotでは、カラム方向に書かれる。
tpm
##      [,1] [,2]
## [1,]  0.7  0.3
## [2,]  0.2  0.8
barplot(tpm[1,], col=c("green","blue"), names.arg=c("happy","blue"))
title("barplot(tpm[1,])")
barplot(tpm[2,], col=c("green","blue"), names.arg=c("happy","blue"))
title("barplot(tpm[2,])")
##ちょっと、変だね。「棒」を足しても、1.0にならん。
barplot(tpm, col=c("green","blue"), names.arg=c("happy","blue"))
title("barplot(tpm)")
##こうすると、「棒」が足して1.0。
barplot(t(tpm), col=c("green","blue"), names.arg=c("happy","blue"))
title("barplot(t(tpm))")

##さあ、PageUpとPageDownで昔のを見てみよう。…

##おまけのおまけ、、、
rbind(tpm,tpm)
##      [,1] [,2]
## [1,]  0.7  0.3
## [2,]  0.2  0.8
## [3,]  0.7  0.3
## [4,]  0.2  0.8
matplot(rbind(tpm,tpm), col=c("green","blue"), type="b")
        


history <- c()
nState <- ncol(tpm) ### 状態の数
state <- 1 ##最初は、いい気分からStart
for (i in 1:10000) {
  ##今の気分(state)から、次の気分(newState)を決める。
  ##sampleの引数は、(x, size, replace = FALSE, prob = NULL)
  newState <- sample(1:nState, 1, prob = tpm[state, ])
  ##cf:いい気分からは、次の状態の確率は、0.7, 0.3
  ##probは、選択する確率。
  ##ちょっと表示。
  barplot(t(tpm), col=c("green","blue"), names.arg=c("happy","blue"))
  text(state, c(0.1,0.8)[newState], "*",col="white",cex=10)
  title(paste("traial", i))
  state <- newState
  history <- c( state, history )
}
history
###各気分にどの位の確率でいるか。
hist(history)

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


##どの位続けて、「いい気分」状態に留まっているか
##ホントは、もっと、賢いProgrammingはありそうだが、、、
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$mids, log(a$counts))
title(sprintf("log(counts at state %d)",1))
## $midsは、BINの中心値。$countsは、BIN中の要素数。

## 重みづけ線形フィット。
lsf <- lsfit(a$mids[a$counts!=0], log(a$counts)[a$counts!=0],
             wt=log(a$counts)[a$counts!=0])
##log(0)は、うまくないので省く。x軸もy軸も。
abline(coef=lsf$coef)
lsf$coef
##0.6896095。
##「いい気分がn回連続した」後に「もう一度いい気分にいるのは」
##0.69近いので、元々の0.7に近いでしょ。
exp(lsf$coef[2])
##Rでのlogは自然対数。
log(2.718281828)

##どの位続けて、「ブルーな気分」状態に留まっているか
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)
a <- hist(lenhistory, breaks=range(lenhistory)[2]+1)
plot(a$mids, log(a$counts))
title(sprintf("log(counts at state %d)",2))
## $midsは、BINの中心値。$countsは、BIN中の要素数。
## 重みづけ線形フィット。
lsf <- lsfit(a$mids[a$counts!=0], log(a$counts)[a$counts!=0],
             wt=log(a$counts)[a$counts!=0])
##log(0)は、うまくないので省く。x軸もy軸も。
abline(coef=lsf$coef)
lsf$coef

##でと、外から観測できるのは表情。笑顔1、曇り顔2、泣き顔3。ラベル。
##ラベルの出現頻度は、気分によって変わるはず
## observation probability
###「遷移行列」じゃ、無い。

Rho <- matrix(c(0.7, 0.2, 0.1,
                0.1, 0.4, 0.5), 3,2,byrow=FALSE)
faceS <- c("Laugh", "Cloudy", "Weep")
rownames(Rho) <- 1:3
colnames(Rho) <- c("happy", "blue")
(Rho)

## いい気分:笑顔(0.7)、曇り顔(0.2)、泣き顔(0.1)
## ブルーな気分:笑顔(0.1)、曇り顔(0.4)、泣き顔(0.5)

##コマンドでHMMパッケージの準備
##setwd(choose.dir())
##install.packages("hmm.discnp_0.2-4.zip", repos=NULL) ## repos=NULLがlocalからってこと。
##これは、インターネットからインストールする方法。
##install.packages("hmm.discnp", repo="http://cran.ism.ac.jp/", dependencies=TRUE)


##メニューバーの[パッケージ]-[ローカルにあるzipファイルからのパッケージ
##のインストール]を選択し、パッケージのzipファイルを指定
library(hmm.discnp)

##sim.hmm(nsim,tpm,Rho,ispd=NULL,yval=NULL,verb=FALSE)
##Simulate discrete data from a hidden Markov model.
##連続9回の表情の変遷。を10試行。
##第一引数は、「連続9回」を、10試行する。
##スタートは、「いい気分」状態から始める。ispd
faceS <- c("Laugh", "Cloudy", "Weep")
y.sim <- sim.hmm(rep(9,10),tpm,Rho,ispd=c(1,0), yval=1:3)
y.sim
##出力は、表情。笑顔(1)、曇り顔(2)、泣き顔(3)。
## [[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

##listを辞めて、全部つなげたvectorに変換する。
unlist(y.sim)
hist(unlist(y.sim))

##各状態にいる確率は、正しいのかね。

##historyは、さっきの1万回の「気分状態」のシミュレーション結果だった。
barplot(hist(history,n=2)$counts,
        col=c("green","blue"), names.arg=c("happy","blue"))
(hist(history,n=2)$counts)/length(history)
##n=2は、BINの数を2ケにする。$countsで各BINの中の個数。
##/length(history)で、確率になる。
##[1] 0.3955 0.6045

## Rhoは、気分によって、笑顔、曇り顔、泣き顔、が出てくる確率。
Rho %*% c(0.3955, 0.6045)

hist(unlist(y.sim))
length(unlist(y.sim))
##全部で90回だから、*90なのよ。
lines(1:3, Rho %*% c(0.3955, 0.6045)*90, type="b",col=2)


##さて、毎日の表情から、心を読むぞ
##同じ様にデータを作る。連続40回の表情の変遷。を100試行。
y.sim <- sim.hmm(rep(40,100),tpm,Rho,ispd=c(1,0), yval=1:3)

##MarkovModelのステートの数(気分の種類数)Kが2だとすると。
##try <- hmm(y.sim,faceS,K=2,verb=T,stationary=T)
try <- hmm(y.sim, yval=1:3, K=2, verb=T, stationary=T)
try$ispd
try$Rho
try$tpm

##「尤もらしい、気分の遷移」を推定。
##アンドリュー・ビタビが考えた、アルゴリズム。
try.vit <- viterbi(y.sim, object=try, ispd=c(1,0))
##try.vit <- viterbi(y.sim, tpm=try$tpm, Rho=try$Rho, ispd=c(1,0))
hist(unlist(try.vit))
##なんか変だな。。

barplot(t(rbind(tpm,try$tpm)),
        names.arg=t(outer(c("true", "estimated"),1:2,paste)),
        col=c("green","blue"))
title("transition probability matrix")
      
barplot(cbind(Rho,try$Rho),
        names.arg=t(outer(c("true", "estimated"),1:2,paste)))
title("observation probability matrix")


##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","blue"))

##なんと無く似ている??

##気分数4の場合の各状態でのラベルの出現頻度は、???。
barplot(try$Rho,names.arg=paste("estimated",1:4))
title("observation probability matrix")

##いかさまに、気分、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)))
title("observation probability matrix")



##他のモデルに当てはめると、どうなるか??
##受理の確率を知りたい。。。
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から始まる。
mean(log(pr(s.vit, y.sim, tpm=tpm, Rho=Rho, ispd=c(1,0))))
hist(log(pr(s.vit, y.sim, tpm=tpm, Rho=Rho, ispd=c(1,0))))
## -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)
rownames(Rho.uso) <- 1:3

s.vit.uso <- viterbi(y.sim, tpm=tpm.uso, Rho=Rho.uso, ispd=c(1,0))
s.vit.uso##気分のリスト。ispd=c(1,0)なので、すべての配列は気分1から始まる。
s.vit##気分のリスト。ispd=c(1,0)なので、すべての配列は気分1から始まる。

mean(log(pr(s.vit.uso, y.sim, tpm=tpm, Rho=Rho, ispd=c(1,0))))
hist(log(pr(s.vit.uso, y.sim, tpm=tpm, Rho=Rho, ispd=c(1,0))))
## -24.34478


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

####################################
######  DeepLearningの前ふり  ######
####################################
##setwd(choose.dir())
##文字イメージデータの読み込み
mojiDat <- as.matrix(read.table("t10k-images.txt"))/255
range(mojiDat)##0-255
dim(mojiDat)
## [1] 10000   784 ;;10000個の文字データ、28x28のます。784=28*28

mojiLab <- unlist(read.table("t10k-labels.txt"))
length(mojiLab)
## [1] 10000 ;;10000個の文字データのラベル。0~9
table(mojiLab)
##   0    1    2    3    4    5    6    7    8    9 
## 980 1135 1032 1010  982  892  958 1028  974 1009 

for( i in 1:100 ){
  ##i番目の文字データ(728の実数)を読み込んで。
  tmp <- matrix(mojiDat[i,],ncol=28)
  ##縦方向が下から積みあがっているので、ちょっと換える。
  image(x=1:28, y=1:28, z=tmp[,28:1],col=gray(seq(0,1,len=256)))
  par(new=T)
  contour(x=1:28, y=1:28, z=tmp[,28:1],col=2)
  title(sprintf("%d  (%d bannme)", mojiLab[i], i))
}

source("snnMoji.txt")
##実は前回まで使っていたのと、全く同じのニューラルネット
##違いはニューロンの数の定義のところのみ。
## # network definition;
## # layers
## Input <- layer(28*28);##入力層は、28*28
## Hidden <- layer(30);##中間層は、30個。
## Output <- layer(10);##出力層は、0~9を当てるので、10個。

dim(mojiDat)
## [1] 10000   784 ;;10000個の文字データ、28x28のます。784=28*28

##入力バッファの形式は、
##BufInput[<入力層の個数分>,<サンプルの個数分>]。なので、転置。
BufInput <- t(mojiDat)
dim(BufInput)


##出力バッファの形式は、
##BufTarget[<出力層の個数分>,<サンプルの個数分>]。
BufTarget <- matrix(0, nrow=10,ncol=dim(BufInput)[2])
dim(BufTarget)##全部、0にしたっつ。
BufTarget[,1:4]
mojiLab[1:4]

for( i in seq(dim(BufInput)[2]) ){
  if(mojiLab[i]==0){##おバカ。「0」の時に10にする
    BufTarget[10,i] <- 1; 
  } else {
    BufTarget[mojiLab[i],i] <- 1
  }
}
BufTarget[,1:4]
mojiLab[1:4]

##以上で準備完了っ。

e<-c()
gcc(0)
##いざ学習。
date()
for(i in 1:10){
  for(iSet in order(seq(dim(BufTarget)[2]))){
    learn(iSet);
    if(iSet%%100==0){print(iSet)};
    e<-c(e,Error);
  }
}
date()

plot(e)

aa <- seq(dim(BufInput)[2])
aa[sapply(aa, activate)>0.1]
plot(sapply(aa, activate))

i <- 592
activate(i)
matplot(cbind(Output$out,target), type="b")

##i番目の文字データ(728の実数)を読み込んで。
tmp <- matrix(mojiDat[i,],ncol=28)
##縦方向が下から積みあがっているので、ちょっと換える。
image(x=1:28, y=1:28, z=tmp[,28:1],col=gray(seq(0,1,len=256)))
par(new=T)
contour(x=1:28, y=1:28, z=tmp[,28:1],col=2)
title(sprintf("%d  (%d bannme)", mojiLab[i], i))



##ちょっとズレ、「1」を例にして。
par(mfrow=c(1,1))
for( i in seq(dim(BufInput)[2])[1:1000] ){
  ##if(mojiLab[i]==1){
  if(mojiLab[i]==7){
    ##i番目の文字データ(728の実数)を読み込んで。
    tmp <- matrix(mojiDat[i,],ncol=28)
    ##縦方向が下から積みあがっているので、ちょっと換える。
    image(x=1:28, y=1:28, z=tmp[,28:1],col=gray(seq(0,1,len=256)))
    par(new=T)
    contour(x=1:28, y=1:28, z=tmp[,28:1],col=2)
    title(sprintf("%d  (%d bannme)", mojiLab[i], i))
    ##画像のモーメントを計算する。
    m00 <- sum(tmp)
    bb1 <- matrix(1:28,ncol=28,nrow=28)
    bb1t <- t(bb1)

    ##生の1次モーメントを計算して
    m10 <- sum(tmp*bb1)
    m01 <- sum(tmp*bb1t)

    ##生の2次モーメントを計算して
    bb2 <- matrix((1:28)^2,ncol=28,nrow=28)
    m20 <- sum(tmp*bb2)
    bb2t <- t(bb2)
    m02 <- sum(tmp*bb2t)
    m11 <- sum(tmp*bb1*bb1t)

    ##重心周りの2次モーメントに変換して
    u20 <- m20/m00 - (m10/m00)^2
    u02 <- m02/m00 - (m01/m00)^2
    u11 <- m11/m00 - (m10/m00)*(m01/m00)

    b=1/(1/2*atan(2*u11/(u20-u02)))##傾きを計算して、
    y0 <- m01/m00
    x0 <- m10/m00

    abline(v=m10/m00, col="white")
    abline(h=m01/m00, col="white")
    ##(x0,y0)を通り、傾きbの直線。
    ##ablineは、「a+b*x」である。
    abline(a=y0-b*x0, b=1/(1/2*atan(2*u11/(u20-u02))),col="yellow",lwd=3)
  }
}