# -*- mode: R -*- ##EmacsとESS使いの方に。モードラインの設定。 ##2016年10月14日(金曜日) ##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は、1,2,3のいずれか。 state <- sample(1:3, 1, prob=c(0.1,0.7,0.2)) poi <- c( state, poi ) } hist(poi) ## 1:いい気分。20%の確率でブルーな気分(red)に変化する ## 2:ブルーな気分。30%の確率で安定状態(green)に変化する ###気分の遷移関数。 ###遷移行列。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("happy","blue")) title("transition probability matrix") ##BarPlotでは、カラム方向に書かれる。 tpm ## [,1] [,2] ## [1,] 0.7 0.3 ## [2,] 0.2 0.8 barplot(tpm[1,], col=c("green","red"), names.arg=c("happy","blue")) title("barplot(tpm[1,])") barplot(tpm[2,], col=c("green","red"), names.arg=c("happy","blue")) title("barplot(tpm[2,])") barplot(tpm, col=c("green","red"), names.arg=c("happy","blue")) title("barplot(tpm)") ##おまけ、、、 matplot(tpm, col=c("green","red"), type="b") 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, ]) ## 状態からの変遷 ## いい気分では、次の状態の確率は、[1,] 0.7 0.3 ##probは、選択する確率。 history <- c( state, history ) } 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$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 exp(lsf$coef[2]) log(2.718281828) ###どの位続けて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) 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 exp(lsf$coef[2]) ## 0.7987281 ##でと、外から観測できるのは表情。笑顔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パッケージの準備 ##CUIでやる時。 ##setwd(choose.dir()) ##install.packages("hmm.discnp_0.2-4.zip", repos=NULL) ## repos=NULLがlocalからってこと。 ##GUIからやる時。 ##メニューバーの[パッケージ]-[ローカルにある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日間の表情。 ##第一引数は、「連続9日間」を、10試行する。 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 ##listを辞めて、全部つなげたvectorに変換する。 unlist(y.sim) hist(unlist(y.sim)) ##各状態にいる確率は、正しいのかね。 (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)) 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)) ##MarkovModelのステートの数(気分の種類数)Kが2だとすると。 ##こっから、が、 ##try <- hmm(y.sim,K=2,verb=T,stationary=T) try <- hmm(y.sim,K=2,verb=T,stationary=F, cis=F) try$Rho try$tpm ##dim(try$ispd) ##matplot(t(try$ispd),type="b") barplot(t(rbind(tpm,try$tpm)), names.arg=t(outer(c("true", "estimated"),1:2,paste)), col=c("green","red")) title("tpm,try$tpm") barplot(cbind(Rho,try$Rho), names.arg=t(outer(c("true", "estimated"),1:2,paste))) title("Rho,try$Rho") ## tpmとRhoがなんとなく合っている!! ##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))) ##他のモデルに当てはめると、どうなるか?? ##受理の確率を知りたい。。。 ##ちょっと多目にデータを作って、、、 y.sim <- sim.hmm(rep(40,100),tpm,Rho,ispd=c(1,0)) ##ビタビ・アルゴリズム ##「観察された表情の時系列(y.sim)」から、気分の時系列を推定 s.vit <- viterbi(y.sim, tpm=tpm, Rho=Rho, ispd=c(1,0)) table(unlist(y.sim))## ## 1 2 3 ##1506 1236 1258 table(unlist(s.vit)) ## 1 2 ##1681 2319 s.vit##気分のリスト。ispd=c(1,0)なので、気分の時系列は1から始まる。 y.sim[[1]]##「観察された表情の時系列(y.sim)」の1試行目 ## [1] 1 3 1 1 2 2 3 3 1 2 1 1 1 2 1 3 .... ## やる気になれば、ずっと「いい気分」状態で、 ##笑顔、曇り顔、泣き顔を出せる。でも、泣き顔は「ブルーな気分」で多い。 rownames(Rho) <- 1:3 ##no 'dimnames' attribute for arrayは、これだった。 mean(log(pr(s.vit, y.sim, tpm=tpm, Rho=Rho))) ##prは、y.simが受け入れられる確率。 ## -6.81935 ##別の表情機械。tpm.uso, Rho.uso 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)) mean(log(pr(s.vit.uso, y.sim, tpm=tpm.uso, Rho=Rho.uso))) ## -24.34478 ##2種類の塩基配列を生成する ## state1->state2は、0.5, と、state1->state5も、0.5。 ## 次の位置>> 1 2 3 4 5 6 7 8 tpm <- matrix(c(0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, ##1<<元の位置 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, ##2 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, ##3 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, ##4 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, ##5 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, ##6 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, ##7 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0), ##8 8,8,byrow=TRUE) barplot(t(tpm),col=1:8, names.arg=1:8) ## 1A 2T 3G 4C 5s 6e a <- 0.2 Rho <- matrix(c(0.0, 0.0, 0.0, 0.0, 1.0, 0.0, ##1 1/4-a, 1/4-a, 1/4+a, 1/4+a, 0.0, 0.0, ##2 ストリーム@GCリッチ 1/4-a, 1/4-a, 1/4+a, 1/4+a, 0.0, 0.0, ##3 ストリーム@GCリッチ 1/4-a, 1/4-a, 1/4+a, 1/4+a, 0.0, 0.0, ##4 ストリーム@GCリッチ 1/4+a, 1/4+a, 1/4-a, 1/4-a, 0.0, 0.0, ##5 ストリームAATリッチ 1/4+a, 1/4+a, 1/4-a, 1/4-a, 0.0, 0.0, ##6 ストリームAATリッチ 1/4+a, 1/4+a, 1/4-a, 1/4-a, 0.0, 0.0, ##7 ストリームAATリッチ 0.0, 0.0, 0.0, 0.0, 0.0, 1.0),##8 6,8,byrow=FALSE) ## ##スタートとエンドが有る。 ##s→*→*→*→e、てな感じの配列。 y <- sim.hmm(rep(5,1000), tpm, ispd=c(1.0, rep(0,7)), Rho, yval=c("A","T","G","C","s","e")) ## リストの行列化 sim <- do.call("rbind", y) table(sim) sim ## [,1] [,2] [,3] [,4] [,5] ## [1,] "s" "G" "C" "G" "e" ## [2,] "s" "C" "T" "C" "e" dim(sim) ## sim ## A C e G s T ## 768 739 1000 740 1000 753 sum(table(sim)) ##2種類の塩基配列が出来ているか?? ##試しに、ATなら-1点,GCなら1点。合計が正ならGCリッチ。 zz<-function(x){switch(EXPR = x, "A"=-1, "T"=-1, "G"=1, "C"=1, 0)}##ATGC以外の場合は0 zzz<-function(xx){mean(sapply(xx,zz))} hist(apply(sim,1,zzz),main="Histgram of GC score") 1; ##################################### ###### HMMをつかった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) ##すみません。時間切れ、、、 ・3種類学習させて、 ・未知の配列をいれたら、どのHMMが一番受理の確率が高いか?を計算して、判 定する。 ・ハズ。