嗯，終于迎來了終極大Boss。

1. 概論

想要理解清楚隱馬爾可夫模型確實(shí)比之前要難一些，但是我盡量講清楚。先來看一個(gè)問題：
假設(shè)天氣的狀況分為：晴天、多云、雨天。
我想預(yù)測明天的天氣狀況。
這里我們假設(shè)：每天的天氣狀況只跟前幾天的天氣狀況有關(guān)系，而不去考慮其他影響因素，例如風(fēng)力、氣壓等等。—— 這就是著名的馬爾可夫假設(shè)。即：假設(shè)模型的當(dāng)前狀態(tài)僅僅依賴于前面的幾個(gè)狀態(tài)?！?dāng)然，這種假設(shè)非常粗糙，并且因此可能將一些非常重要的信息丟失，但是它缺極大的簡化了問題。
一個(gè)馬爾科夫過程是狀態(tài)間的轉(zhuǎn)移僅依賴于前n個(gè)狀態(tài)的過程。這個(gè)過程被稱之為n階馬爾科夫模型，其中n是影響下一個(gè)狀態(tài)選擇的（前）n個(gè)狀態(tài)。最簡單的馬爾科夫過程是一階模型，它的狀態(tài)選擇僅與前一個(gè)狀態(tài)有關(guān)。為了把問題更加簡化，我們后面都將基于一階馬爾可夫過程來討論。

根據(jù)一階馬爾可夫過程，明天的天氣狀況，只依賴于今天的天氣狀況。由于從晴天到多云，從晴天到雨天，從晴天到晴天的變化不是確定的，也就是說，晴天的后面不一定還是晴天，也不一定是雨天，也不一定是多云，任何一個(gè)狀態(tài)都有可能是所有狀態(tài)的下一個(gè)轉(zhuǎn)移狀態(tài)，我們假設(shè)每一個(gè)狀態(tài)轉(zhuǎn)移都有一個(gè)概率值，稱為狀態(tài)轉(zhuǎn)移概率——這是從一個(gè)狀態(tài)轉(zhuǎn)移到另一個(gè)狀態(tài)的概率，那么我們有多少種狀態(tài)轉(zhuǎn)移的可能性呢？按上面的例子，我們總共有如下幾種：
晴天->晴天；
晴天->雨天；
晴天->多云；
雨天->晴天；
雨天->雨天；
雨天->多云；
多云->晴天；
多云->雨天；
多云->多云；
所以，總共是種狀態(tài)轉(zhuǎn)移過程。我們把上面這幾種過程的概率都列出來，就形成了一個(gè)狀態(tài)轉(zhuǎn)移概率矩陣，如下所示：

注意，每一行的總和為1。第一行表示，今天是晴天，那么明天也是晴天的概率為0.5，今天是晴天，明天是多云的概率是0.375，今天是晴天，明天是雨天的概率是0.125。
上面的概率轉(zhuǎn)移矩陣，我們用數(shù)學(xué)公式來表示為：

注意，此時(shí)是按列來看的。第一列代表今天是晴天，明天是晴天、多云、雨天的概率，依次類推。

假設(shè)今天是晴天，我們初始化一個(gè)向量，來表示今天是晴天：
，

那么

就代表明天是晴天、多云、雨天的概率。

有了明天的天氣概率，我們要預(yù)測后天的概率要怎么算呢？顯然，后天的概率就是用明天的天氣概率去乘以概率轉(zhuǎn)移矩陣，也即：

同理，大后天的天氣概率為：

可以看到，我們?nèi)绻A(yù)測未來n天的天氣，我們是通過今天的天氣一步一步向前預(yù)測推導(dǎo)的。后面每一天的天氣狀態(tài)都只跟前一天的天氣狀態(tài)相關(guān)。如此形成了一個(gè)鏈條，這就是馬爾可夫鏈。

好了，現(xiàn)在我們?cè)賮砜词裁词请[馬爾可夫鏈。

假設(shè)有一個(gè)被關(guān)在地牢里的家伙，他無法看到外面的天氣狀況。但是他可以看到地牢里的干濕狀態(tài)：干燥、稍干、潮濕、濕潤。地牢的干濕狀態(tài)與天氣的狀況有一定的概率關(guān)系。那么這個(gè)被關(guān)在地牢里的家伙他能否預(yù)測未來的天氣情況呢？
我們把地牢的干濕程度記為觀測狀態(tài)（用O表示），對(duì)于地牢的人來說，天氣狀態(tài)就是隱藏狀態(tài)（用S表示），用圖表示為：

這就是隱馬爾可夫模型。

從上面的介紹中，我們知道，一個(gè)隱馬爾可夫模型由如下五個(gè)元素：

1. 隱含狀態(tài) S。這些狀態(tài)之間滿足馬爾可夫性質(zhì)，是馬爾可夫模型中實(shí)際所隱含的狀態(tài)。這些狀態(tài)通常無法通過直接觀測而得到。
   （例如上面的天氣狀態(tài)）
2. 可觀測狀態(tài) O。在模型中與隱含狀態(tài)相關(guān)聯(lián)，可通過直接觀測而得到。（例如上面的地牢狀態(tài)）
3. 初始狀態(tài)概率矩陣 。表示隱含狀態(tài)在初始時(shí)刻t=1的概率矩陣
4. 隱含狀態(tài)轉(zhuǎn)移概率矩陣 A。
5. 觀測狀態(tài)轉(zhuǎn)移概率矩陣 B。有的地方也叫混淆矩陣，表示給定一個(gè)隱藏狀態(tài)后得到的觀察狀態(tài)的概率。

一般的，可以用三元組來簡潔的表示一個(gè)隱馬爾可夫模型。隱馬爾可夫模型實(shí)際上是標(biāo)準(zhǔn)馬爾可夫模型的擴(kuò)展，添加了可觀測狀態(tài)集合和這些狀態(tài)與隱含狀態(tài)之間的概率關(guān)系。

對(duì)于隱馬爾可夫模型，我們根據(jù)實(shí)際情況的不同，一般會(huì)遇到3類問題：
1) 評(píng)估模型
給定觀測序列 和模型參數(shù)，怎樣有效計(jì)算某一觀測序列的概率，進(jìn)而可對(duì)該HMM做出相關(guān)評(píng)估。例如，已有一些模型參數(shù)各異的HMM，給定觀測序列 ，我們想知道哪個(gè)HMM模型最可能生成該觀測序列。通常我們利用forward算法分別計(jì)算每個(gè)HMM產(chǎn)生給定觀測序列O的概率，然后從中選出最優(yōu)的HMM模型。
這類評(píng)估的問題的一個(gè)經(jīng)典例子是語音識(shí)別。在描述語言識(shí)別的隱馬爾科夫模型中，每個(gè)單詞生成一個(gè)對(duì)應(yīng)的HMM，每個(gè)觀測序列由一個(gè)單詞的語音構(gòu)成，單詞的識(shí)別是通過評(píng)估進(jìn)而選出最有可能產(chǎn)生觀測序列所代表的讀音的HMM而實(shí)現(xiàn)的。
2) 解碼問題
給定觀測序列 和模型參數(shù)，怎樣尋找某種意義上最優(yōu)的隱狀態(tài)序列。在這類問題中，我們感興趣的是馬爾科夫模型中隱含狀態(tài)，這些狀態(tài)不能直接觀測但卻更具有價(jià)值，通常利用Viterbi算法來尋找。
這類問題的一個(gè)實(shí)際例子是中文分詞，即把一個(gè)句子如何劃分其構(gòu)成才合適。例如，句子“發(fā)展中國家”是劃分成“發(fā)展-中-國家”，還是“發(fā)展-中國-家”。這個(gè)問題可以用隱馬爾科夫模型來解決。句子的分詞方法可以看成是隱含狀態(tài)，而句子則可以看成是給定的可觀測狀態(tài)，從而通過建HMM來尋找出最可能正確的分詞方法。
3) 學(xué)習(xí)問題
即HMM的模型參數(shù)未知，如何調(diào)整這些參數(shù)以使觀測序列 的概率盡可能的大。通常使用Baum-Welch算法以及Reversed Viterbi算法解決。

2. 評(píng)估模型

我們先來看第一個(gè)問題，即評(píng)估模型。評(píng)估模型主要解決的問題是：給定模型和觀測序列 ，計(jì)算模型下觀測到序列O出現(xiàn)的概率。
對(duì)于求解這個(gè)問題，我們可以使用窮舉法。還是以上面那個(gè)例子，天氣的狀態(tài)轉(zhuǎn)移矩陣我們已經(jīng)知道了，為：

假設(shè)，這三種天氣導(dǎo)致地牢干燥和潮濕的概率為（為了簡化問題，我們這里只假設(shè)地牢只有兩種狀態(tài)：干燥和潮濕）：

如果我們要求地牢的干濕程度為O={干燥、潮濕、潮濕}的概率，我們要這么計(jì)算：
如果三天都為晴天，得到上述觀測序列的概率為
如果三天的天氣狀況為：晴天，晴天，多云，則得到上面觀測序列的概率為：
如果三天的天氣狀況為：晴天，晴天，雨天，則得到上面觀測序列的概率為：
如果三天的天氣狀況為：晴天，多云，晴天，則得到上面觀測序列的概率為：
……
最后的得到概率應(yīng)該為上面窮舉的所有概率之和。
當(dāng)然，這里我們還沒有計(jì)算分別得到上面這所有天氣情況下的概率。這樣計(jì)算的話，將會(huì)非常復(fù)雜。在觀察值較多的情況下，顯然是不太現(xiàn)實(shí)的。因此我們要用到《數(shù)據(jù)驅(qū)動(dòng)的決策分析》中學(xué)到的動(dòng)態(tài)規(guī)劃思想。

動(dòng)態(tài)規(guī)劃
動(dòng)態(tài)規(guī)劃的思想，其實(shí)就是一個(gè)遞歸思想。

以最短路徑問題為例，來說明動(dòng)態(tài)規(guī)劃的基本思想。
從A到E的最短路徑S(A)，可以轉(zhuǎn)化為三個(gè)性質(zhì)完全相同，但規(guī)模較小的子問題，即分別從B1、B2、B3到E的最短路徑，記為S(B1), S(B2), S(B3)。

同樣，計(jì)算S(B1)、 S(B2) 、S(B3)又可以歸結(jié)為性質(zhì)完全相同，但規(guī)模更小的問題，即分別求C1，C2，C3到E的最短路徑問題S(Ci) (i=1, 2, 3)。

S(D1)和S(D2)是以知的，它們分別是：S(D1)=5，S(D2)=2。因而，可以從這兩個(gè)值開始，逆向遞歸計(jì)算S(A)的值。

向前算法
現(xiàn)在回到我們的問題。我們現(xiàn)在要評(píng)估觀測數(shù)據(jù)O出現(xiàn)的概率。我們可以這樣分解：
先看第一天地牢為干燥的概率。因?yàn)槲覀冇腥N天氣狀態(tài)，每種天氣狀態(tài)都有可能會(huì)導(dǎo)致地牢干燥。因此第一天地牢干燥的概率為：

假設(shè)：，也即第一天是晴天，那么。根據(jù)上面的狀態(tài)轉(zhuǎn)移矩陣，我們知道，，代表晴天到干燥的概率：P(干燥|晴天)， 分別代表P(干燥|多云)，P(干燥|雨天)。
我們記：
因此，

再來看第二天地牢為潮濕的概率。
當(dāng)時(shí)，輸出為。假設(shè)第二天是晴天（記為），那么地牢為潮濕的概率為：
，記為：
注意，這里的是狀態(tài)轉(zhuǎn)移矩陣A中的元素，代表從晴天、多云、雨天轉(zhuǎn)移到晴天的概率，不是我們上面定義的。
同理，假設(shè)第二天是多云，那么地牢為潮濕的概率為：
，記為：
假設(shè)第二天是雨天，那么地牢為潮濕的概率為：
，記為：
那么，第二天為潮濕的概率就為：

同理，我們來看第三天為潮濕的概率
當(dāng)時(shí)，輸出為。
當(dāng)?shù)谌鞛榍缣鞎r(shí)：

當(dāng)?shù)谌鞛槎嘣茣r(shí)：

當(dāng)?shù)谌鞛槎嘣茣r(shí)：

最后，我們得到最終的概率：

由此，我們得到向前算法的一般性步驟（公式）：

step1: 初始化。（i表示隱藏狀態(tài)的計(jì)數(shù)，例如上面的天氣狀態(tài)：晴天，多云，雨天）

step2: 計(jì)算

step3:

3. 解碼問題

給定觀測序列 和模型參數(shù)，找到最可能的狀態(tài)序列
同上面的思路，我們可以使用窮舉法，列出來所有可能導(dǎo)致地牢干燥、潮濕、潮濕的天氣狀態(tài)組合以及產(chǎn)生地牢干燥、潮濕、潮濕狀態(tài)的概率。（在上面一節(jié)中提到的），然后通過比較這些概率的大小，找到概率最大的那個(gè)天氣狀態(tài)組合即可。
顯然，當(dāng)數(shù)據(jù)量較大時(shí)，窮舉法在時(shí)間上也是不可接受的。同樣，我們也考慮使用遞歸算法來解決上述問題。這就是維比特算法。

維比特算法
維比特算法的思路跟上一節(jié)中，我們介紹的動(dòng)態(tài)規(guī)劃中的思路是一致的，只不過在動(dòng)態(tài)規(guī)劃中，我們求的是路徑最小值，現(xiàn)在我們求的是概率最大值。因此這里就不再贅述了。直接給出步驟和公式：

step1: 初始化。即計(jì)算在t=1的時(shí)刻，達(dá)到某狀態(tài)的最大概率（這個(gè)局部概率我們用來表示）。i表示隱藏狀態(tài)的計(jì)數(shù)，例如上面的天氣狀態(tài)：晴天，多云，雨天：。

step2: 遞推。計(jì)算t>1時(shí)刻的局部概率。在一階馬爾可夫假設(shè)下，狀態(tài)X在一個(gè)狀態(tài)序列后發(fā)生的概率只取決于之前的一個(gè)狀態(tài)。還是以上面的那個(gè)例子來說，在t時(shí)刻，天氣狀態(tài)為X的概率，只取決于t-1時(shí)刻的狀態(tài)。對(duì)于t-1時(shí)刻的狀態(tài)，我們有三種：晴天、多云、雨天。那么，從晴天到達(dá)t時(shí)刻X的概率為：
在t時(shí)刻X狀態(tài)下的觀察概率為：
同樣的，從雨天到達(dá)t時(shí)刻X的觀察概率為：
所以，

step3: 計(jì)算的最大值。獲得此最大值的觀測序列就是我們要找的觀測序列。

4. 學(xué)習(xí)問題

好，終于來到了我們要重點(diǎn)關(guān)注，也就是機(jī)器學(xué)習(xí)中經(jīng)常用到，但是也是最復(fù)雜的一個(gè)問題，就是已知觀測序列 ，估計(jì)模型參數(shù)，使得在該模型參數(shù)下出現(xiàn)該觀測序列的概率最大。
學(xué)習(xí)問題要根據(jù)觀測序列來推導(dǎo)模型參數(shù)，這一類的問題對(duì)應(yīng)到概率論中的極大似估計(jì)問題。但是這里是有隱含變量的極大似然估計(jì)，因此直接無法通過直接求導(dǎo)進(jìn)行求解，而要通過EM算法來求解這一類問題。
EM算法是一類算法，用于解決有隱含變量的概率模型參數(shù)的極大似然估計(jì)，具體到隱馬爾可夫模型中的具體算法是 Baum-Welch 算法（也叫向前-向后算法）。
直觀地理解EM算法，它也可被看作為一個(gè)逐次逼近算法：事先并不知道模型的參數(shù)，可以隨機(jī)的選擇一套參數(shù)或者事先粗略地給定某個(gè)初始參數(shù)λ0 ，確定出對(duì)應(yīng)于這組參數(shù)的最可能的狀態(tài)，計(jì)算每個(gè)訓(xùn)練樣本的可能結(jié)果的概率，在當(dāng)前的狀態(tài)下再由樣本對(duì)參數(shù)修正，重新估計(jì)參數(shù)λ ，并在新的參數(shù)下重新確定模型的狀態(tài)，這樣，通過多次的迭代，循環(huán)直至某個(gè)收斂條件滿足為止，就可以使得模型的參數(shù)逐漸逼近真實(shí)參數(shù)。
由于這個(gè)算法實(shí)在太復(fù)雜，這里就不細(xì)展開了。感興趣的同學(xué)可以在網(wǎng)上找一些相資料進(jìn)行學(xué)習(xí)。

如果是要做一個(gè)調(diào)包俠，個(gè)人覺得只看概論部分就應(yīng)該可以了。
好了，還是老規(guī)矩，最后一部分，讓我們用R代碼來結(jié)束本章內(nèi)容，也是機(jī)器學(xué)習(xí)系列的最后內(nèi)容。

5.R代碼的實(shí)現(xiàn)

例子1: 預(yù)測基因序列中的啟動(dòng)因子序列。
這個(gè)例子是：首先給定一堆已經(jīng)區(qū)分好了啟動(dòng)子和非啟動(dòng)子的基因序列，我們要找到某種規(guī)律或者模型，使得下次我們?cè)偬峁┮环N基因序列時(shí)，這個(gè)模型能夠很好的告訴我們，該序列是否屬于啟動(dòng)子序列。
這個(gè)例子中，不存在隱藏狀態(tài)。我們?yōu)榱耸褂玫紿MM模型，人為定義了一個(gè)隱藏狀態(tài)。該隱藏狀態(tài)的序列和觀測狀態(tài)的序列一樣，都是樣本中的基因序列。且，從隱藏狀態(tài)到觀測狀態(tài)，相同的堿基之間轉(zhuǎn)換的概率為1，不同的堿基之間轉(zhuǎn)換的概率為0。
另，基因序列是由四個(gè)不同的堿基通過不同的排列組合構(gòu)成。即，一個(gè)基因序列有四個(gè)狀態(tài)，分別為：a、c、g、t。為了區(qū)分基因序列的起止，我們認(rèn)為加入了兩個(gè)狀態(tài)。S代表基因序列的開始，X代表基因序列的結(jié)束。對(duì)于初始狀態(tài)，我們認(rèn)為只有S。因此初始狀態(tài)矩陣為c=(1,0,0,0,0,0)。
現(xiàn)在HMM模型中的五個(gè)要素，我們已經(jīng)有了四個(gè)。剩下的，我們只需要計(jì)算出來隱藏狀態(tài)轉(zhuǎn)移矩陣A即可。
有了模型之后，我們只需要使用向前算法，即可計(jì)算出觀測序列出現(xiàn)的概率。理論上，我們?nèi)绻褂玫姆菃?dòng)因子估計(jì)出的模型參數(shù)，那么在非啟動(dòng)因子的基因序列上，出現(xiàn)該基因序列的概率會(huì)比較大，而對(duì)于啟動(dòng)因子的基因序列，出現(xiàn)該基因序列的概率會(huì)比較小。由此我們可以使用訓(xùn)練數(shù)據(jù)，采用交叉檢驗(yàn)法來衡量模型的好壞。
好了，以上就是下面這個(gè)例子的總體思想。現(xiàn)在我們來看代碼：

#預(yù)測啟動(dòng)子基因序列。
promoters <- read.csv("study/graduate/機(jī)器學(xué)習(xí)/promoters.data", header=F, dec=",", strip.white=TRUE, stringsAsFactors = FALSE)
#V1:+表示啟動(dòng)子，-表示非啟動(dòng)子
#V2:特定序列的識(shí)別符
#V3: 核苷酸序列本身（即觀測序列，同時(shí)也是隱藏序列）
promoters[1,]

#數(shù)據(jù)拆分成啟動(dòng)子和非啟動(dòng)子。
#subset表示根據(jù)條件來篩選子集。這里的3表示選擇第3列，丟棄掉1，2列
positive_observations <- subset(promoters, V1=='+', 3)
negative_observations <- subset(promoters, V1=='-', 3)

#在每個(gè)序列之前加上S，表示數(shù)據(jù)的開始，在每個(gè)結(jié)尾加上X，表示數(shù)據(jù)的結(jié)束
positive_observations<-sapply(positive_observations, function(x) paste("S",x,"X",sep=""))
negative_observations<-sapply(negative_observations, function(x) paste("S",x,"X",sep=""))
positive_observations[1]

#將一串字符串，按每個(gè)字母分開,得到一個(gè)列表
positive_observations<-strsplit(positive_observations,"")
negative_observations<-strsplit(negative_observations,"")
head(positive_observations[1])

#定義隱藏狀態(tài)。我們的隱藏狀態(tài)有以下6個(gè)（S,X,a,c,g,t）
states <- c("S", "X", "a", "c", "g", "t")
#定義觀測狀態(tài)。我們的觀測狀態(tài)也有6個(gè)
symbols <- c("S", "X", "a", "c", "g", "t")
#定義初始概率，也即pi，表示剛開始的時(shí)候，只有S，即從S開始。
startingProbabilities<-c(1,0,0,0,0,0)
#定義發(fā)射概率，也即觀測狀態(tài)轉(zhuǎn)移概率矩陣 B。這里是一個(gè)對(duì)角矩陣，對(duì)角線都為1，其余為0
#代表只能從S到S，從a到a，從c到c，從g到g，從t到t
emissionProbabilities<-diag(6)
colnames(emissionProbabilities)<-states #列名為狀態(tài)值
rownames(emissionProbabilities)<-symbols #行名為觀測值
emissionProbabilities

#根據(jù)我們的樣本，來計(jì)算轉(zhuǎn)移概率矩陣。
calculateTransitionProbabilities <- function(data, states) {
  #初始化轉(zhuǎn)移狀態(tài)矩陣，剛開始都為0
    transitionProbabilities<-matrix(0,length(states),length(states))
    colnames(transitionProbabilities)<-states
    rownames(transitionProbabilities)<-states
  #按順序循環(huán)數(shù)據(jù)中的每一個(gè)字母，統(tǒng)計(jì)從狀態(tài)A狀態(tài)B的次數(shù)，并放到矩陣中。即匯總出每一個(gè)狀態(tài)轉(zhuǎn)移的總數(shù)
    for (index in 1:(length(data)-1)) {
        current_state <- data[index]
        next_state <- data[index+1]
        transitionProbabilities[current_state, next_state] <- transitionProbabilities[current_state, next_state] + 1
    }
    #根據(jù)次數(shù)矩陣，使用sweep函數(shù)來計(jì)算概率。
    #swepp中的第二個(gè)參數(shù)1，代表按行計(jì)算；2，代表按列計(jì)算
    #rowSums函數(shù)代表按行進(jìn)行求和
    #整個(gè)函數(shù)的意思就是，用每一行的元素除以這一行的數(shù)據(jù)總和，得到歸一化的概率；
    transitionProbabilities<-sweep(transitionProbabilities, 1, rowSums(transitionProbabilities), FUN="/")
    return(transitionProbabilities)
}

#這里相當(dāng)于把negative_observations從原來的列表轉(zhuǎn)換為了一個(gè)向量
#即，把原來每一行數(shù)據(jù)，都進(jìn)行了連接，形成了一個(gè)長向量，向量長度為原來的行數(shù)*每行的字母個(gè)數(shù)
negative_observation<-Reduce(function(x,y) c(x, y), negative_observations, c())
#使用非啟動(dòng)子基因序列的樣本來訓(xùn)練轉(zhuǎn)移概率矩陣。
(transitionProbabilitiesNeg  <- calculateTransitionProbabilities(negative_observation, states))

library("HMM")
#使用HMM的五個(gè)參數(shù)來初始化隱馬爾可夫模型
negative_hmm = initHMM(states, symbols, startProbs=startingProbabilities, transProbs=transitionProbabilitiesNeg, emissionProbs=emissionProbabilities)

#交叉檢驗(yàn)。每次我們從啟動(dòng)子基因序列的樣本中拿出一條數(shù)據(jù)作為檢驗(yàn)數(shù)據(jù)，剩下的作為訓(xùn)練集數(shù)據(jù)
#先用剩下的數(shù)據(jù)來構(gòu)建HMM模型，此時(shí)的模型是根據(jù)啟動(dòng)子基因序列來構(gòu)建的
#然后，我們用剩下的這條檢驗(yàn)數(shù)據(jù)，分別帶入到之前negative_hmm模型和現(xiàn)在的positive_hmm模型中，采用向前算法來計(jì)算出現(xiàn)該基因序列的概率
#理論上，因?yàn)檫@條檢驗(yàn)數(shù)據(jù)是來自于啟動(dòng)子基因序列，因此，通過negative_hmm模型計(jì)算出來的概率應(yīng)該小于使用positive_hmm模型計(jì)算出來的概率。
#統(tǒng)計(jì)錯(cuò)誤率，我們即可知道模型的好壞
incorrect<-0 #定義錯(cuò)誤次數(shù)統(tǒng)計(jì)變量
for (obs in 1:length(positive_observations)) {
  #使用啟動(dòng)子基因序列，構(gòu)造啟動(dòng)子基因序列模型，方便后面進(jìn)行對(duì)比  
  positive_observation<-Reduce(function(x,y) c(x, y), positive_observations[-obs], c())
  transitionProbabilitiesPos  <- calculateTransitionProbabilities(positive_observation, states)
  positive_hmm = initHMM(states, symbols, startProbs=startingProbabilities, transProbs=transitionProbabilitiesPos, emissionProbs=emissionProbabilities)

  #獲取檢驗(yàn)數(shù)據(jù)
  test_observation<-positive_observations[[obs]]
  final_index<-length(test_observation)
  
  #將檢驗(yàn)數(shù)據(jù)分別帶入到兩個(gè)模型中，使用向前算法，分別計(jì)算使用這兩個(gè)模型得到該序列的概率
  pos_probs<-exp(forward(positive_hmm,test_observation))
  neg_probs<-exp(forward(negative_hmm,test_observation))
  
  #計(jì)算出現(xiàn)該序列的概率之和。（上面得到的是序列中每個(gè)堿基出現(xiàn)的概率）
  pos_seq_prob<-sum(pos_probs[,final_index])
  neg_seq_prob<-sum(neg_probs[,final_index])
  
  #如果使用positive模型得到的概率小于使用negative模型得到概率，則錯(cuò)誤數(shù)+1
  if (pos_seq_prob<neg_seq_prob) incorrect<-incorrect+1
}

#反過來，使用positive數(shù)據(jù)構(gòu)建模型，使negative數(shù)據(jù)進(jìn)行檢驗(yàn)。因此以下代碼不再做解釋
positive_observation <- Reduce(function(x,y) c(x, y), positive_observations, c())
transitionProbabilitiesPos  <- calculateTransitionProbabilities(positive_observation, states)
positive_hmm = initHMM(states, symbols, startProbs=startingProbabilities, transProbs=transitionProbabilitiesPos, emissionProbs=emissionProbabilities)
  for (obs in 1:length(negative_observations)) {
  negative_observation<-Reduce(function(x,y) c(x, y), negative_observations[-obs], c())
  transitionProbabilitiesNeg <- calculateTransitionProbabilities(negative_observation, states)
  negative_hmm = initHMM(states, symbols, startProbs=startingProbabilities, transProbs=transitionProbabilitiesNeg, emissionProbs=emissionProbabilities)
  
  test_observation<-negative_observations[[obs]]
  final_index<-length(test_observation)
  
  pos_probs<-exp(forward(positive_hmm,test_observation))
  neg_probs<-exp(forward(negative_hmm,test_observation))
  
  pos_seq_prob<-sum(pos_probs[,final_index])
  neg_seq_prob<-sum(neg_probs[,final_index])
  
  if (pos_seq_prob>neg_seq_prob) incorrect<-incorrect+1
}

#最后計(jì)算交叉檢驗(yàn)后的準(zhǔn)確度
(cross_validation_accuracy <- 1 - (incorrect/nrow(promoters)))

例子2: 分析英語單詞的字母順序模式（尋找隱藏序列的patten）。
每個(gè)英語單詞，其字母的組成順序都是有一定規(guī)律的。我們很少看到在一個(gè)單詞中，有兩個(gè)元音字母連在一起。我們是否可以找到一個(gè)模型，來體現(xiàn)出來這種規(guī)律呢？
我們還是使用我們?cè)谏弦恢v中用到的文本數(shù)據(jù)。

#分析英語單詞的字母順序模式（尋找隱藏序列的patten）
library(ggplot2)
library("tm")
#讀取文件
nb_pos <- VCorpus(DirSource(path_to_pos_folder), readerControl = list(reader = reader(DirSource(path_to_pos_folder)), language = "en"))
nb_neg <- VCorpus(DirSource(path_to_neg_folder), readerControl = list(reader = reader(DirSource(path_to_neg_folder)), language = "en"))
nb_all <- c(nb_pos,nb_neg) #合并兩個(gè)數(shù)據(jù)
nb_all <- tm_map(nb_all, content_transformer(tolower)) #全部轉(zhuǎn)換為小寫

#因?yàn)閚b_all的格式不好處理，我們將nb_all中的文本信息提取出來，存入到texts變量中
texts <- sapply(1 : length(nb_all), function(x) nb_all[[x]])

#文本處理
texts<-sapply(texts, function(x) gsub("\\s","W", x)) #把空格用W代替
texts<-sapply(texts, function(x) gsub("[0-9]","N", x)) #把數(shù)字用N代替
texts<-sapply(texts, function(x) gsub("[[:punct:]]","P", x)) #把標(biāo)點(diǎn)符號(hào)用P代替
texts<-sapply(texts, function(x) gsub("[^a-zWNP]","O", x)) #其他符號(hào)用O代替
texts
#因?yàn)閠exts比較大，我們只取前40條數(shù)據(jù)，并將其每個(gè)字母打散，存入到列表中
big_text_splits<-lapply(texts[1:40], function(x) strsplit(x, ""))
big_text_splits<-unlist(big_text_splits, use.names = F) #將列表轉(zhuǎn)換為向量

#由于我們不知道單詞構(gòu)造到隱藏狀態(tài)，因此這里強(qiáng)行安了3個(gè)狀態(tài)
states <- c("s1", "s2", "s3")
numstates <- length(states)
symbols <- c(letters,"W","N","P","O") #觀測變量就是26個(gè)英文字母加上我們前面定義的幾個(gè)特殊字母
numsymbols <- length(symbols)

#初始化一個(gè)初始狀態(tài)概率矩陣pi。因?yàn)槲覀兊碾[藏狀態(tài)也是隨意給的，因此這里的初始狀態(tài)概率矩陣也可以隨機(jī)給一個(gè)
set.seed(124124)
startingProbabilities <- matrix(runif(numstates), 1, numstates)
startingProbabilities<-sweep(startingProbabilities, 1, rowSums(startingProbabilities), FUN="/")

#初始化一個(gè)轉(zhuǎn)移狀態(tài)矩陣
set.seed(454235)
transitionProbabilities<-matrix(runif(numstates*numstates),numstates,numstates)
transitionProbabilities<-sweep(transitionProbabilities, 1, rowSums(transitionProbabilities), FUN="/")

#初始化一個(gè)觀測狀態(tài)轉(zhuǎn)移概率矩陣（發(fā)射矩陣）
set.seed(923501)
emissionProbabilities<-matrix(runif(numstates*numsymbols),numstates,numsymbols)
emissionProbabilities<-sweep(emissionProbabilities, 1, rowSums(emissionProbabilities), FUN="/")

#初始化隱馬爾可夫模型
library("HMM")
hmm <- initHMM(states, symbols, startProbs = startingProbabilities, transProbs = transitionProbabilities, emissionProbs = emissionProbabilities)
#使用Baum-Welch 算法來學(xué)習(xí)參數(shù)。
hmm_trained <- baumWelch(hmm, big_text_splits)

hmm_trained$hmm

#根據(jù)學(xué)習(xí)后的參數(shù)，我們來畫出轉(zhuǎn)移概率矩陣。轉(zhuǎn)移概率矩陣體現(xiàn)了在一個(gè)單詞中，每個(gè)字母出現(xiàn)的概率
#畫出從狀態(tài)S1到各字母之間的概率
p1 <- ggplot(data=data.frame(hmm_trained$hmm$emissionProbs[1,]), aes(x = names(hmm_trained$hmm$emissionProbs[1,]), y = hmm_trained$hmm$emissionProbs[1,]))
p1 <- p1 + geom_bar(stat="identity")
p1 <- p1 + ggtitle("Symbol Emission Probabilities for State 1")
p1 <- p1 + theme(plot.title = element_text(lineheight=.8, face="bold", vjust=2))
p1 <- p1 + xlab("State")  
p1 <- p1 + ylab("Emission Probability") 
p1

#畫出從狀態(tài)S2到各字母之間的概率 
p2 <- ggplot(data=data.frame(hmm_trained$hmm$emissionProbs[2,]), aes(x = names(hmm_trained$hmm$emissionProbs[2,]), y = hmm_trained$hmm$emissionProbs[2,]))
p2 <- p2 + geom_bar(stat="identity")
p2 <- p2 + ggtitle("Symbol Emission Probabilities for State 2")
p2 <- p2 + theme(plot.title = element_text(lineheight=.8, face="bold", vjust=2))
p2 <- p2 + xlab("State")  
p2 <- p2 + ylab("Emission Probability") 
p2

#畫出從狀態(tài)S3到各字母之間的概率
p3 <- ggplot(data=data.frame(hmm_trained$hmm$emissionProbs[3,]), aes(x = names(hmm_trained$hmm$emissionProbs[3,]), y = hmm_trained$hmm$emissionProbs[3,]))
p3 <- p3 + geom_bar(stat="identity")
p3 <- p3 + ggtitle("Symbol Emission Probabilities for State 3")
p3 <- p3 + theme(plot.title = element_text(lineheight=.8, face="bold", vjust=2))
p3 <- p3 + xlab("State")  
p3 <- p3 + ylab("Emission Probability") 
p3

#獲得隱含狀態(tài)轉(zhuǎn)移概率矩陣
(trained_transition_probabilities<-hmm_trained$hmm$transProbs)

#模擬給定隱馬爾可夫模型的狀態(tài)和觀測路徑。即根據(jù)模型構(gòu)造一個(gè)句子（單詞）
set.seed(9898)
simHMM(hmm_trained$hmm, 30)

完結(jié)～

【參考文章】
機(jī)器學(xué)習(xí)——HMM(隱馬爾可夫模型的基本概念)
HMM學(xué)習(xí)最佳范例一：介紹
隱馬爾可夫模型(HMM) - 1 - 基本概念
通俗理解，隱馬爾可夫模型的前前后后
HMM-前向后向算法理解與實(shí)現(xiàn)
【NLP】隱馬爾可夫模型三個(gè)基本問題相關(guān)算法實(shí)現(xiàn)