使用隱馬爾科夫模型生成資料

    為了理解HMMs,讓我們是靠一下例子，一個銷售人員因為其工作，需要轉輾在三個城市之間（倫敦，巴塞羅納，紐約）。他的目的是最小化旅行的時間，以便更高效的工作。考慮一下他的工作委託和計劃安排，我們能得到他從x城市到y城市的概率。在一下給出的資訊中，P(X->Y)表示從X城市到Y城市的概率：

P(倫敦->倫敦) =0.10

P(倫敦->巴塞羅那) = 0.70 

P(倫敦->紐約) = 0.20

P(巴塞羅那->巴塞羅那) =0.15

P(巴塞羅那->倫敦) = 0.75

P(巴塞羅那->紐約) = 0.10

P（紐約->紐約）=0.05

P（紐約->倫敦）= 0.60

P（紐約->巴塞羅那） = 0.35

將上面資料轉化為矩陣：

                倫敦    巴塞羅那    紐約

     倫敦     0.10     0.70        0.20

巴塞羅那    0.75    0.15        0.10

      紐約     0.60    0.35        0.05

我們現在得到了所有的資訊，讓我們繼續問題陳述.這個銷售人員在週二這天從倫敦出發開始了他的旅行，，並且它需要在週五計劃一些事情，但是這依靠他勝出何地。他將出現在巴塞羅那的概率是多少呢？這個表格將幫我們弄清楚：

如果我們沒有一個 馬爾可夫鏈來對這個問題建模，我們將不會知道他的日程安排是怎樣的。我們的目標是給出他在被給定日期他在特定城市的可能性。假設我們將轉換矩陣設為T,當前日期為X(i),則

X(i+1)=X(i)*T

在這個案例中，週五離週二是3天。這意味著我們必須計算X(i+3)

X(i+1)=X(i)*T

X(i+2)=X(i+1)*T

X(i+3)=X(i+2)*T

可推出：X(i+3)=X(i)*T^3

我們需要設定X(i) 

X（i）=[0.10 0.70 0.20]

下一步就是計算矩陣的立方。網上有很多線上對矩陣進行操作的工具例如 http://matrix.reshish.com/multiplication.php。計算結果如下：

T^3 =

則在週二的概率是

P(倫敦)= 0.31

P(巴塞羅那)=0.53

P(NY)=0.16

我們能夠考到巴塞羅那的概率最高，這也具有地理意義，因為巴塞羅那相比於紐約，離倫敦是更近的。那如何構建一個HMMS？

建立一個python文件，程式碼如下：

import datetime
import numpy as np
import matplotlib.pyplot as plt
from hmmlearn.hmm import GaussianHMM
from Slicing_time_series_data import read_data
#　　　hmmlearn實現了三種HMM模型類，按照觀測狀態是連續狀態還是離散狀態，可以分為兩類。GaussianHMM和GMMHMM是連續觀測狀態的HMM模型，而MultinomialHMM是離散觀測狀態的模型，也是我們在HMM原理系列篇裡面使用的模型。
#載入資料
data = np.loadtxt('data_1D.txt',delimiter=',')


#提取第三列進行訓練
x=np.column_stack([data[:,2]])


#建立GaussianHMM，引數：狀態的數量n_components=5,covariance_type='diag',“diag” - 每個狀態使用對角協方差矩陣。n_iter要執行的最大迭代次數。
num_components =5
hmm =GaussianHMM(n_components=num_components,covariance_type='diag',n_iter=1000)


#訓練HMM
print("\n正在訓練隱馬爾科夫模型....")
hmm.fit(x)
#輸出每個HMM狀態的平均值和方差
print("\n均值和方差：")
for i in range(hmm.n_components):
    print('\n隱狀態',i+1)
    print('均值 = ',round(hmm.means_[i][0],2))
    print("方差 = ",round(np.diag(hmm.covars_[i])[0],2))


#生成1200條資料訓練HMM模型並繪出
num_samples = 1200
generated_data,_=hmm.sample(num_samples)#_約定不關心數字的變數，後期不使用
plt.plot(np.arange(num_samples),generated_data[:,0],c='red')
plt.title('Gnenerate data')
plt.show()

1200資料如下

輸出結果如下:結果與書上不同，本人環境使用的是hmmlearn0.20版本，較之前訓練資料不同了