## # -*- mode: R -*- ##EmacsとESS使いの方に。モードラインの設定。 ##2019年10月11日(金曜日) ##参照ファイルは、 ## forPrint6.pdf ## command.txt##これ本体。 ## ・read.table ## ・install ## ・source ############ ######################### ###### 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)##countを書く ##何故、この様に歯抜けなの??? 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)を作ってみると、 ## 0.5 1.5 2.5 3.5 hist(poi, breaks=1:4-0.5, label=T) hist(poi, breaks=1:4-0.5, plot=F) hist(poi, breaks=1:4-0.5, plot=F)$counts ##plot=Fは、ヒストグラムを描画しない。 ##「histogramを描画するに必要だったデータ」を教えてくれる。 ## 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, ] ## 0.7 0.3 ##1は「いい気分」。 ##次の状態は、0.7の確率でいい気分(1)に留まる、0.3の確率でブルーな気分(2)に変化。 tpm[2, ] ## 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 (transition probability matrix) 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では、カラムが1つの折れ線になる。 matplot(rbind(tpm,tpm), col=c("green","blue"), type="b") ##気分の移り変わりの歴史 history <- c() nState <- ncol(tpm) ### 状態の数 state <- 1 ##最初は、いい気分からStart for (i in 1:1000) { ##今の気分(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") text(c(10,10), c(1,2), c("happy", "blue")) title("histgram of kibun (happy or blue)") ##どの位続けて、「いい気分」状態に留まっているか ##ホントは、もっと、賢い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パッケージの準備 ##インターネットからインストールする方法。 install.packages("hmm.discnp", repo="http://cran.ism.ac.jp/", dependencies=TRUE) ##メニューバーの[パッケージ]-[ローカルにあるzipファイルからのパッケージ ##のインストール]を選択し、パッケージのzipファイルを指定 library(hmm.discnp) ##Simulate discrete data from a hidden Markov model. ##連続40回の表情の変遷。を100試行。 ##第一引数ylengthsは、「連続40回」を、100試行する。 ##スタートは、「いい気分」状態から始める。ispd faceS <- c("Laugh", "Cloudy", "Weep") y.sim <- rhmm(ylengths=rep(40,100),tpm=tpm,Rho=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 ##....... ** ./100 ** ## [1] 1 1 1 1 1 3 2 3 3 ##listを辞めて、全部つなげたvectorに変換する。 unlist(y.sim) hist(unlist(y.sim)) ##各状態にいる確率は、正しいのかね。 ##historyは、さっきの「気分状態」のシミュレーション結果だった。 barplot(hist(history,n=2)$counts, col=c("green","blue"), names.arg=c("happy","blue")) (kibunProb <- (hist(history,n=2)$counts)/length(history)) ##n=2は、BINの数を2ケにする。$countsで各BINの中の個数。 ##/length(history)で、確率になる。 ## Rhoは、気分によって、笑顔、曇り顔、泣き顔、が出てくる確率。 ## これで、なべての、笑顔、曇り顔、泣き顔、が出てくる確率。 Rho %*% kibunProb hist(unlist(y.sim)) length(unlist(y.sim)) ##全部で4000回だから、*4000なのよ。 ##一致する。ねと。 lines(1:3, (Rho %*% kibunProb)*4000, type="b",col=2) ##さて、本題、、、「毎日の表情から、心を読むぞ」 ##同じ様にデータを作る。連続40回の表情の変遷。を100試行。 y.sim <- rhmm(ylengths=rep(40,100),tpm=tpm,Rho=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, newstyle=FALSE) try$ispd try$Rho try$tpm ##「尤もらしい、気分の遷移」を推定。 ##アンドリュー・ビタビが考えた、アルゴリズム。 try.vit <- viterbi(y.sim, try, 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, newstyle=FALSE) t(try$tpm) try$Rho par(mfrow=c(1,1)) barplot(t(tpm), col=c("green","blue"), names.arg=paste("true",1:2)) ##気分数4の場合の遷移行列。tpm (transition probability matrix) barplot(t(try$tpm), names.arg=paste("estimated",1:4)) title("try$tpm for K=4") ##気分数4の場合の各状態でのラベルの出現頻度は、???。 barplot(try$Rho,names.arg=paste("estimated",1:4)) title("observation probability matrix") title("try$Rho for K=4") ##気分数4の場合の遷移行列で、 ##いかさまに、気分、2,3,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)) outer(c("true", "estimated"),1:2,paste) 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 <- rhmm(ylengths=rep(40,100),tpm=tpm,Rho=Rho,ispd=c(1,0), yval=1:3) ##本当のパラメータ 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.49, 0.51, 0.48, 0.52),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