干货！隐马尔科夫模型
公式推導
Hmmlearn
GaussianHMM
GMMHMM
MultinomialHMM
股票走勢預測
特征準備
建立模型
可視化短線預測
參考資料
HMM
公式推導
在 HMM 中，有兩個基本假設：
齊次 Markov 假設（未來只依賴于當前）：?
?
觀測獨立假設：?
?
HMM 要解決三個問題：
Evaluation：，Forward-Backward 算法
Learning：，EM 算法（Baum-Welch）
Decoding：，Vierbi 算法
預測問題：
濾波問題：
Evaluation
??
??
根據齊次 Markov 假設：
??
所以：
又由于：
于是：
??
我們看到，上面的式子中的求和符號是對所有的觀測變量求和，于是復雜度為 ??。
下面，記 ??，所以，??。我們看到：
??
對 ??：
??
利用觀測獨立假設：
??
上面利用了齊次 Markov 假設得到了一個遞推公式，這個算法叫做前向算法。
還有一種算法叫做后向算法，定義?，?：
??
對于這個 ??：
??
于是后向地得到了第一項。
Learning
為了學習得到參數的最優值，在 MLE 中：
??
我們采用 EM 算法（在這里也叫 Baum Welch 算法），用上標表示迭代：
??
其中，???是觀測變量，??是隱變量序列。于是：
??
?這里利用了 ?和?無關。將 Evaluation 中的式子代入：
??
對 ??：
??
上面的式子中，對???求和可以將這些參數消掉：
??
上面的式子還有對???的約束??。定義 Lagrange 函數：
??
于是：
??
對上式求和：
??
所以：
??
Decoding
Decoding 問題表述為：
??
我們需要找到一個序列，其概率最大，這個序列就是在參數空間中的一個路徑，可以采用動態規劃的思想。
定義：
??
于是：
??
這個式子就是從上一步到下一步的概率再求最大值。記這個路徑為：
??
Hmmlearn
hmmlearn中有三種隱馬爾可夫模型:GaussianHMM、GMMHMM、MultinomialHMM。它們分別代表了觀測序列的不同分布類型。
GaussianHMM
適合用于可見層狀態是連續類型且假設輸出概率符合Gaussian分布的情況
class?hmmlearn.hmm.GaussianHMM(n_components=1,?covariance_type='diag',min_covar=0.001,startprob_prior=1.0,?transmat_prior=1.0,?means_prior=0,??????means_weight=0,covars_prior=0.01,covars_weight=1,algorithm='viterbi',random_state=None,?n_iter=10,?tol=0.01,verbose=False,?params='stmc',?init_params='stmc')參數:
n_components：整數，指定了隱藏層結點的狀態數量。
covariance_type：字符串，輸出概率的協方差矩陣類型。可選值：
'spherical'：球面協方差矩陣，即矩陣對角線上的元素都相等，且其他元素為零
'diag'：對角協方差矩陣，即對角線元素可以是任意值，其他元素為零。
'full'：完全矩陣，即任意元素可以是任意值。
'tied'：所有狀態都是用同一個普通的方差矩陣。
min_covar：浮點數，協方差矩陣中對角線上的最小數值，該值設置得越小模型對數據的擬合就越好，但更容易出現過度擬合。
startprob_prior：數組，形狀為(n_components, )。初始狀態的先驗概率分布。
transmat_prior：數組，形狀為(n_components, n_components )。先驗的狀態轉移矩陣。
means_prior,means_weight：先驗隱藏層均值矩陣。
covars_prior、covars_weight：先驗隱藏層協方差矩陣
algorithm：字符串，指定了Decoder算法。可以為 'viterbi'（維特比算法）或者'map'，map比viterbi更快速但非全局最優解。
random_state：隨機種子，用于在Baum-Welch算法中初始化模型參數。
n_iter:Baum-Welch算法最大迭代次數。該值越大，訓練模型對數據的擬合度越高，但訓練耗時越長。
tol：指定迭代收斂閾值。該值越小(必須＞=0),訓練模型對數據的擬合度越高，但訓練耗時越長。
verbose：是否打印Baum-Welch每次迭代的調試信息
params：字符串,在訓練過程中更新哪些HMM參數。可以是四個字母中的任意幾個組成的字符串。
's'：初始概率。
't'：轉移概率。
'm'：均值。
'c'：偏差。
init_params：字符串,在訓練開始之前使用哪些已有概率矩陣初始化模型。
's'：初始概率。
't'：轉移概率。
'm'：均值。
'c'：偏差。
屬性：
n_features：整數，特征維度。
monitor_：ConvergenceMonitor對象，可用它檢查EM算法的收斂性。
transmat_：矩陣，形狀為 (n_components, n_components)，是狀態之間的轉移概率矩陣。
startprob_：數組，形狀為(n_components, )，是初始狀態的概率分布。
means_：一個數組，形狀為(n_components,n_features )，每個狀態的均值參數。
covars_：數組，每個狀態的方差參數，其形狀取決于方差類型：
'spherical'：形狀為(n_components, ) 。
'diag'：形狀為(n_components,n_features ) 。
'full'：形狀為(n_components, n_features, n_features) 。
'tied'：形狀為(n_features,n_features ) 。
方法：
decode(X, lengths=None, algorithm=None)：已知觀測序列X尋找最可能的狀態序列。
logprob：浮點數，代表產生的狀態序列的對數似然函數。
state_sequence：一個數組，形狀為(n_samples, )，代表狀態序列。
X：array-like，形狀為 (n_samples, n_features)。指定了觀測的樣本。
lengths：array-like，形狀為 (n_sequences, )。指定了觀測樣本中，每個觀測序列的長度，其累加值必須等于n_samples 。
algorithm：字符串，指定解碼算法。必須是'viterbi'（維特比）或者'map'。
參數：
返回值：
fit(X, lengths=None)：根據觀測序列 X，來訓練模型參數。
在訓練之前會執行初始化的步驟。如果你想避開這一步，那么可以在構造函數中通過提供init_params關鍵字參數來避免。
參數：X，lengths 參考 decode() 方法。
返回值：self對象。
predict(X, lengths=None)：已知觀測序列X，尋找最可能的狀態序列。
參數：X，lengths 參考 decode() 方法。
返回值：數組，形狀為(n_samples, )，代表狀態序列。
predict_proba(X, lengths=None)：計算每個狀態的后驗概率。
參數：X，lengths 參考 decode() 方法。
返回值：數組，代表每個狀態的后驗概率。
sample(n_samples=1, random_state=None)：從當前模型中生成隨機樣本。
X：觀測序列，長度為n_samples 。
state_sequence：狀態序列，長度為n_samples 。
n_samples：生成樣本的數量。
random_state：指定隨機數。如果為None，則使用構造函數中的random_state。
參數：
返回值：
score(X, lengths=None)：計算預測結果的對數似然函數。
參數：X，lengths 參考 decode() 方法。
返回值：預測結果的對數似然函數。
GMMHMM
混合高斯分布的隱馬爾可夫模型，適合用于隱藏層狀態是連續類型且假設輸出概率符合GMM 分布（Gaussian Mixture Model，混合高斯分布）的情況
原型
hmmlearn.hmm.GMMHMM(n_components=1,n_mix=1,startprob_prior=1.0,transmat_prior=1.0,covariance_type='diag',?covars_prior=0.01,?algorithm='viterbi',random_state=None,n_iter=10,tol=0.01,?verbose=False,params='stmcw',init_params='stmcw')大部分參數與GaussianHMM中的含義一樣，不再重復說明，不同點有如下兩個:
n_mix：整數，指定了混合高斯分布中高斯分布的數量，如果n_min等于1，則GMMHMM退化為GaussianHMM。
means_prior,means_weight,covars_prior,covars_weight：這四個參數雖然名稱與GaussianHMM一致，但其維數隨著n_mix的設置而不同。
屬性：
n_features：整數，特征維度。
monitor_：ConvergenceMonitor對象，可用它檢查EM算法的收斂性。
transmat_：矩陣，形狀為 (n_components, n_components)，是狀態之間的轉移概率矩陣。
startprob_：數組，形狀為(n_components, )，是初始狀態的概率分布。
gmms_：列表，指定了每個狀態的混合高斯分布的分模型。
方法：參考hmmlearn.hmm.GaussianHMM 。
MultinomialHMM
多項式分布的隱馬爾可夫模型,適合用于可見層狀態是離散類型的情況。
原型為：
class?hmmlearn.hmm.MultinomialHMM(n_components=1,?startprob_prior=1.0,?transmat_prior=1.0,algorithm='viterbi',?random_state=None,?n_iter=10,tol=0.01,verbose=False,?params='ste',?init_params='ste')參數：整數，參考hmmlearn.hmm.GaussianHMM。
屬性：
n_features：一個整數，特征維度。
monitor_：一個ConvergenceMonitor對象，可用它檢查EM算法的收斂性。
transmat_：一個矩陣，形狀為 (n_components, n_components)，是狀態之間的轉移概率矩陣。
startprob_：一個數組，形狀為(n_components, )，是初始狀態的概率分布。
emissionprob_：一個數組，形狀為(n_components, n_features)，每個狀態的發射概率。
方法：參考hmmlearn.hmm.GaussianHMM 。
#?MultinomialHMM案例 import?numpy?as?np from?hmmlearn?import?hmm#?發射概率 emission_probability?=?np.array([[0.4,0.3,0.3],?#?晴:打球,讀書,訪友[0.2,0.3,0.5],?#?陰:打球,讀書,訪友[0.1,0.8,0.1]??#?雨:打球,讀書,訪友 ])#?轉移概率 transition_probability?=?np.array([[0.7,0.2,0.1],#?晴?陰?雨[0.3,0.5,0.2],#?陰[0.3,0.4,0.3]?#?雨 ])#?定義隱藏層各狀態的初始概率?晴、陰、雨 start_probability?=?np.array([0.5,0.3,0.2])#?建模 model?=?hmm.MultinomialHMM(n_components=3)?#?3種可見層狀態 model.startprob_?=?start_probability model.transmat_?=?transition_probability model.emissionprob_?=?emission_probability#?0為打球，1為讀書，2為訪友 observe_chain?=?np.array([0,2,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,0]).reshape(-1,1)#?0為晴，1為陰，2為雨 model.predict(observe_chain) array([0, 0, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0])股票走勢預測
本案例從雅虎金融網站獲取真實的歷史交易數據，訓練并分析歷史走勢，嘗試尋找出歷史走勢規律以預測未來。下面通過用HMM模型來預測走勢規律
HMM 時間軸：由于數據模型是日交易信息，所以本模型的時間軸以日為單位，即每一天是一個HMM狀態結點。
可見層特征：選取數據文件中的兩個重要數據作為可見層特征，即收盤漲跌值、交易量。因為這些屬性為連續值，所以選用Gaussian模型。收盤漲跌值特征并未在數據模型中直接體現，它要通過計算前后兩天收盤價的差值獲得
隱藏層狀態：定義三種狀態，對應于漲、平、跌。
預測方式：通過查看隱藏層轉換矩陣的轉換概率，根據最后一個結點的隱藏狀態預測未來一天的漲跌可能。
隱藏層的漲跌狀態只受當天及之前表示層特征的影響，所以其本身并不能決定下一天的走勢。而要預測下一天的隱藏層狀態，需要結合轉換概率矩陣進行分析。
from?datetime?import?datetime import?pandas?as?pd import?numpy?as?np import?matplotlib.pyplot?as?plt from?matplotlib.dates?import?DayLocator from?hmmlearn.hmm?import?GaussianHMM import?warnings warnings.filterwarnings("ignore")%matplotlib?inlinedata?=?pd.read_csv("sample_stock.csv") print(data.shape) data.head() (1258, 7)原始數據是以日為單位的某股票基礎交易數據，各屬性列表示:
日期
開盤價
當日最高價
當日最低價
收盤價
調整收盤
當日交易量
日期和交易量去除第一天的數據,因為第一天會被用來計算第二天的漲跌值特征，所以馬爾可夫鏈實際是從第二天開始訓練的。
dates?=?data['Date'].values[1:] close_v?=?data['Close'].values volume?=?data['Volume'].values[1:]#?計算數組中前后兩個值差額 diff??=?np.diff(close_v) close_v?=?close_v[1:]X?=?np.column_stack([diff,volume]) X.shape (1257, 2)建立模型
n_components 參數指定了使用3個隱藏層狀態
covariance_type定義了協方差矩陣類型為對角線類型，即每個特征的高斯分布有自己的方差參數，相互之間沒有影響
n_iter參數定義了Baum-Welch的最大迭代次數
因為可見層輸入采用了兩個輸入特征（漲跌幅、成交量），所以每個結點有兩個元素，分別代表該狀態的漲跌幅均值、成交量均值。
由于系統的目標預測漲跌幅，所以這里只關心第一個特征的均值，有如下結論：
狀態0的均值是-6.63217241e-03，約為?0.00663，認為該狀態是“平”。
狀態1的均值是4.41582696e-02，約為0.044（漲4分錢），得知該狀態是“漲”。
狀態2的均值是-7.25238348e-02，約為?0.072（跌7分錢），得知該狀態是“跌”。
由漲跌幅特征的方差，可以得到的信息為：
狀態“平”的方差為4.84610308e-02，是三個狀態中方差最小的，也就是說該狀態的預測非常可信。
狀態“漲”的方差為1.42504830e-01，可信度居中。
狀態“跌”的方差為1.10870852e+00，是三個狀態中方差最大的，即該狀態的變化范圍較大，并非很可信。
第一行最大的數值是0.82001546，因此“平”傾向于保持自己的狀態，即第二天仍舊為“平”。
第二行最大的數值是0.72111014，得知“漲”后第二天傾向于變為“漲”。
根據第三行狀態，“跌”后第二天仍舊傾向于“漲”。
根據以上轉換概率，可以看出該股票的長線態勢是非常看好的。
可視化短線預測
X?=?X[-26:] dates?=?dates[-26:] close_v?=?close_v[-26:]hidden_states?=?model.predict(X) hidden_states array([0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2,2, 2, 2, 2]) means_ [-0.0066321724101716766, -0.07252383481726057, 0.044158269609221555] vars_ [0.048461030777305854, 1.1087085156561962, 0.1425048299475001] model.transmat_ array([[0.82001546, 0.01739021, 0.16259432],[0.03793618, 0.22931597, 0.73274785],[0.23556112, 0.04332874, 0.72111014]]) #?計算當前狀態后一天的預測均值 predict_means?=?[] predict_vars?=?[]for?idx,?hid?in?enumerate(range(model.n_components)):comp?=?np.argmax(model.transmat_[idx])?#?最可能的第二天狀態predict_means.append(means_[comp])predict_vars.append(vars_[comp])print('predict_means',predict_means) print('predict_vars',predict_vars) predict_means [-0.0066321724101716766, 0.044158269609221555, 0.044158269609221555] predict_vars [0.048461030777305854, 0.1425048299475001, 0.1425048299475001] #?畫圖 fig,?axs?=?plt.subplots(model.n_components?+?1,?sharex=True,?sharey=True,figsize=(10,10))for?i,?ax?in?enumerate(axs[:-1]):ax.set_title("{0}th?hidden?state".format(i))#?Use?fancy?indexing?to?plot?data?in?each?state.mask?=?hidden_states?==?iyesterday_mask?=?np.concatenate(([False],?mask[:-1]))if?len(dates[mask])?<=?0:continueif?predict_means[i]?>?0.01:#?上漲預測ax.plot_date(dates[mask],?close_v[mask],?"^",?c="#FF0000")elif?predict_means[i]?<?-0.01:#?下跌預測ax.plot_date(dates[mask],?close_v[mask],?"v",?c="#00FF00")else:#?平ax.plot_date(dates[mask],?close_v[mask],?"+",?c="#000000")#?locate?specified?days?of?the?dayax.xaxis.set_minor_locator(DayLocator())ax.grid(True)ax.legend(["Mean:?%0.3f\nvar:?%0.3f"?%?(predict_means[i],?predict_vars[i])],loc='center?left',bbox_to_anchor=(1,?0.5))#?打印真實走勢,用作對比 axs[-1].plot_date(dates,?close_v,?"-",?c='#000000') axs[-1].grid(True) box?=?axs[-1].get_position() axs[-1].set_position([box.x0,?box.y0,?box.width?*?0.8,?box.height]) ax.xaxis.set_minor_locator(DayLocator())#?調整格式 fig.autofmt_xdate() plt.subplots_adjust(left=None,?bottom=None,?right=0.75,?top=None,?wspace=None,?hspace=0.43)plt.show()由圖可知，這個月中大部分隱藏狀態是“平”或者“漲”，其中“漲”的數量更多，由最下方真實走勢圖可見當月結果確實為上漲。最后一天的結點落在了“2th hidden state”中，其預測狀態為上漲（均值為0.044，方差為0.143）。
參考資料
從機器學習到深度學習
hmmlearn文檔
統計學習方法
總結
以上是生活随笔為你收集整理的干货！隐马尔科夫模型的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Python20行代码实现视频字符化
- 下一篇： 真香！一行Python代码，帮你制作小姐