HMM

最新推荐文章于 2023-02-17 19:14:54 发布

天青如水

最新推荐文章于 2023-02-17 19:14:54 发布

阅读量752

点赞数 2

分类专栏： # 机器学习

本文链接：https://blog.csdn.net/qq_16829085/article/details/105970330

版权

机器学习专栏收录该内容

29 篇文章 6 订阅

订阅专栏

HMM

公式推导

在 HMM 中，有两个基本假设：

齐次 Markov 假设（未来只依赖于当前）：
$p(i_{t+1}|i_t,i_{t-1},\cdots,i_1,o_t,o_{t-1},\cdots,o_1)=p(i_{t+1}|i_t)$
观测独立假设：
$p(o_t|i_t,i_{t-1},\cdots,i_1,o_{t-1},\cdots,o_1)=p(o_t|i_t)$

HMM 要解决三个问题：

Evaluation： $p(O|\lambda)$ ，Forward-Backward 算法
Learning： $\lambda=\mathop{argmax}\limits_{\lambda}p(O|\lambda)$ ，EM 算法（Baum-Welch）
Decoding： $I=\mathop{argmax}\limits_{I}p(I|O,\lambda)$ ，Vierbi 算法
1. 预测问题： $p(i_{t+1}|o_1,o_2,\cdots,o_t)$
2. 滤波问题： $p(i_t|o_1,o_2,\cdots,o_t)$

Evaluation

$p(O|\lambda)=\sum\limits_{I}p(I,O|\lambda)=\sum\limits_{I}p(O|I,\lambda)p(I|\lambda)$

$p(I|\lambda)=p(i_1,i_2,\cdots,i_t|\lambda)=p(i_t|i_1,i_2,\cdots,i_{t-1},\lambda)p(i_1,i_2,\cdots,i_{t-1}|\lambda)$

根据齐次 Markov 假设：
$p(i_t|i_1,i_2,\cdots,i_{t-1},\lambda)=p(i_t|i_{t-1})=a_{i_{t-1}i_t}$
所以：
$p(I|\lambda)=\pi_1\prod\limits_{t=2}^Ta_{i_{t-1},i_t}$
又由于：
$p(O|I,\lambda)=\prod\limits_{t=1}^Tb_{i_t}(o_t)$
于是：
$p(O|\lambda)=\sum\limits_{I}\pi_{i_1}\prod\limits_{t=2}^Ta_{i_{t-1},i_t}\prod\limits_{t=1}^Tb_{i_t}(o_t)$
我们看到，上面的式子中的求和符号是对所有的观测变量求和，于是复杂度为 $O(N^T)$ 。

下面，记 $\alpha_t(i)=p(o_1,o_2,\cdots,o_t,i_t=q_i|\lambda)$ ，所以， $\alpha_T(i)=p(O,i_T=q_i|\lambda)$ 。我们看到：
$p(O|\lambda)=\sum\limits_{i=1}^Np(O,i_T=q_i|\lambda)=\sum\limits_{i=1}^N\alpha_T(i)$
对 $\alpha_{t+1}(j)$ ：
$\begin{aligned}\alpha_{t+1}(j)&=p(o_1,o_2,\cdots,o_{t+1},i_{t+1}=q_j|\lambda)\\ &=\sum\limits_{i=1}^Np(o_1,o_2,\cdots,o_{t+1},i_{t+1}=q_j,i_t=q_i|\lambda)\\ &=\sum\limits_{i=1}^Np(o_{t+1}|o_1,o_2,\cdots,i_{t+1}=q_j,i_t=q_i|\lambda)p(o_1,\cdots,o_t,i_t=q_i,i_{t+1}=q_j|\lambda) \end{aligned}$
利用观测独立假设：
$\begin{aligned}\alpha_{t+1}(j)&=\sum\limits_{i=1}^Np(o_{t+1}|i_{t+1}=q_j)p(o_1,\cdots,o_t,i_t=q_i,i_{t+1}=q_j|\lambda)\\ &=\sum\limits_{i=1}^Np(o_{t+1}|i_{t+1}=q_j)p(i_{t+1}=q_j|o_1,\cdots,o_t,i_t=q_i,\lambda)p(o_1,\cdots,o_t,i_t=q_i|\lambda)\\ &=\sum\limits_{i=1}^Nb_{j}(o_t)a_{ij}\alpha_t(i) \end{aligned}$
上面利用了齐次 Markov 假设得到了一个递推公式，这个算法叫做前向算法。

还有一种算法叫做后向算法，定义 $\beta_t(i)=p(o_{t+1},o_{t+1},\cdots，o_T|i_t=i,\lambda)$ ：
$\begin{aligned}p(O|\lambda)&=p(o_1,\cdots,o_T|\lambda)\\ &=\sum\limits_{i=1}^Np(o_1,o_2,\cdots,o_T,i_1=q_i|\lambda)\\ &=\sum\limits_{i=1}^Np(o_1,o_2,\cdots,o_T|i_1=q_i,\lambda)\pi_i\\ &=\sum\limits_{i=1}^Np(o_1|o_2,\cdots,o_T,i_1=q_i,\lambda)p(o_2,\cdots,o_T|i_1=q_i,\lambda)\pi_i\\ &=\sum\limits_{i=1}^Nb_i(o_1)\pi_i\beta_1(i) \end{aligned}$
对于这个 $\beta_1(i)$ ：
$\begin{aligned}\beta_t(i)&=p(o_{t+1},\cdots,o_T|i_t=q_i)\\ &=\sum\limits_{j=1}^Np(o_{t+1},o_{t+2},\cdots,o_T,i_{t+1}=q_j|i_t=q_i)\\ &=\sum\limits_{j=1}^Np(o_{t+1},\cdots,o_T|i_{t+1}=q_j,i_t=q_i)p(i_{t+1}=q_j|i_t=q_i)\\ &=\sum\limits_{j=1}^Np(o_{t+1},\cdots,o_T|i_{t+1}=q_j)a_{ij}\\ &=\sum\limits_{j=1}^Np(o_{t+1}|o_{t+2},\cdots,o_T,i_{t+1}=q_j)p(o_{t+2},\cdots,o_T|i_{t+1}=q_j)a_{ij}\\ &=\sum\limits_{j=1}^Nb_j(o_{t+1})a_{ij}\beta_{t+1}(j) \end{aligned}$
于是后向地得到了第一项。

Learning

为了学习得到参数的最优值，在 MLE 中：
$\lambda_{MLE}=\mathop{argmax}_\lambda p(O|\lambda)$
我们采用 EM 算法（在这里也叫 Baum Welch 算法），用上标表示迭代：
$\theta^{t+1}=\mathop{argmax}_{\theta}\int_z\log p(X,Z|\theta)p(Z|X,\theta^t)dz$
其中， $X$ 是观测变量， $Z$ 是隐变量序列。于是：
$\lambda^{t+1}=\mathop{argmax}_\lambda\sum\limits_I\log p(O,I|\lambda)p(I|O,\lambda^t)\\ =\mathop{argmax}_\lambda\sum\limits_I\log p(O,I|\lambda)p(O,I|\lambda^t)$
这里利用了 $p(O|\lambda^t)$ 和 $\lambda$ 无关。将 Evaluation 中的式子代入：
$\sum\limits_I\log p(O,I|\lambda)p(O,I|\lambda^t)=\sum\limits_I[\log \pi_{i_1}+\sum\limits_{t=2}^T\log a_{i_{t-1},i_t}+\sum\limits_{t=1}^T\log b_{i_t}(o_t)]p(O,I|\lambda^t)$
对 $\pi^{t+1}$ ：
$\begin{aligned}\pi^{t+1}&=\mathop{argmax}_\pi\sum\limits_I[\log \pi_{i_1}p(O,I|\lambda^t)]\\ &=\mathop{argmax}_\pi\sum\limits_I[\log \pi_{i_1}\cdot p(O,i_1,i_2,\cdots,i_T|\lambda^t)] \end{aligned}$
上面的式子中，对 $i_2,i_2,\cdots,i_T$ 求和可以将这些参数消掉：
$\pi^{t+1}=\mathop{argmax}_\pi\sum\limits_{i_1}[\log \pi_{i_1}\cdot p(O,i_1|\lambda^t)]$
上面的式子还有对 $\pi$ 的约束 $\sum\limits_i\pi_i=1$ 。定义 Lagrange 函数：
$L(\pi,\eta)=\sum\limits_{i=1}^N\log \pi_i\cdot p(O,i_1=q_i|\lambda^t)+\eta(\sum\limits_{i=1}^N\pi_i-1)$
于是：
$\frac{\partial L}{\partial\pi_i}=\frac{1}{\pi_i}p(O,i_1=q_i|\lambda^t)+\eta=0$
对上式求和：
$\sum\limits_{i=1}^Np(O,i_1=q_i|\lambda^t)+\pi_i\eta=0\Rightarrow\eta=-p(O|\lambda^t)$
所以：
$\pi_i^{t+1}=\frac{p(O,i_1=q_i|\lambda^t)}{p(O|\lambda^t)}$

Decoding

Decoding 问题表述为：
$I=\mathop{argmax}\limits_{I}p(I|O,\lambda)$
我们需要找到一个序列，其概率最大，这个序列就是在参数空间中的一个路径，可以采用动态规划的思想。

定义：
$\delta_{t}(j)=\max\limits_{i_1,\cdots,i_{t-1}}p(o_1,\cdots,o_t,i_1,\cdots,i_{t-1},i_t=q_i)$
于是：
$\delta_{t+1}(j)=\max\limits_{1\le i\le N}\delta_t(i)a_{ij}b_j(o_{t+1})$
这个式子就是从上一步到下一步的概率再求最大值。记这个路径为：
$\psi_{t+1}(j)=\mathop{argmax}\limits_{1\le i\le N}\delta_t(i)a_{ij}$

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寻找最可能的状态序列。
- 参数：
  - X：array-like，形状为 (n_samples, n_features)。指定了观测的样本。
  - lengths：array-like，形状为 (n_sequences, )。指定了观测样本中，每个观测序列的长度，其累加值必须等于n_samples 。
  - algorithm：字符串，指定解码算法。必须是’viterbi’（维特比）或者’map’。
- 返回值：
  - logprob：浮点数，代表产生的状态序列的对数似然函数。
  - state_sequence：一个数组，形状为(n_samples, )，代表状态序列。
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)：从当前模型中生成随机样本。
- 参数：
  - n_samples：生成样本的数量。
  - random_state：指定随机数。如果为None，则使用构造函数中的random_state。
- 返回值：
  - X：观测序列，长度为n_samples 。
  - state_sequence：状态序列，长度为n_samples 。
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 inline


data = pd.read_csv("./sample_stock.csv")
print(data.shape)
data.head()

(1258, 7)

	Date	Open	High	Low	Close	Adj Close	Volume
0	21/09/2012	9.85	10.09	9.80	9.82	9.528357	66200
1	24/09/2012	9.83	10.07	9.61	9.90	9.605980	44500
2	25/09/2012	9.93	9.99	9.57	9.58	9.295483	54400
3	26/09/2012	9.63	9.72	9.55	9.59	9.305186	29100
4	27/09/2012	9.60	9.71	9.44	9.57	9.285780	55300

原始数据是以日为单位的某股票基础交易数据，各属性列表示:

日期
开盘价
当日最高价
当日最低价
收盘价
调整收盘
当日交易量

# 日期格式转换
data['Date'] = pd.to_datetime(data['Date'], format = "%d/%m/%Y", errors = 'coerce')
data.info()
data.head()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1258 entries, 0 to 1257
Data columns (total 7 columns):
Date         1258 non-null datetime64[ns]
Open         1258 non-null float64
High         1258 non-null float64
Low          1258 non-null float64
Close        1258 non-null float64
Adj Close    1258 non-null float64
Volume       1258 non-null int64
dtypes: datetime64[ns](1), float64(5), int64(1)
memory usage: 68.9 KB

	Date	Open	High	Low	Close	Adj Close	Volume
0	2012-09-21	9.85	10.09	9.80	9.82	9.528357	66200
1	2012-09-24	9.83	10.07	9.61	9.90	9.605980	44500
2	2012-09-25	9.93	9.99	9.57	9.58	9.295483	54400
3	2012-09-26	9.63	9.72	9.55	9.59	9.305186	29100
4	2012-09-27	9.60	9.71	9.44	9.57	9.285780	55300

特征准备

日期和交易量去除第一天的数据,因为第一天会被用来计算第二天的涨跌值特征，所以马尔可夫链实际是从第二天开始训练的。

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的最大迭代次数

model =  GaussianHMM(n_components=3, covariance_type="diag",n_iter=1000)
model.fit(X)

GaussianHMM(algorithm='viterbi', covariance_type='diag', covars_prior=0.01,
      covars_weight=1, init_params='stmc', means_prior=0, means_weight=0,
      min_covar=0.001, n_components=3, n_iter=1000, params='stmc',
      random_state=None, startprob_prior=1.0, tol=0.01, transmat_prior=1.0,
      verbose=False)

print("Means and vars of each hidden state")
means_ = []
vars_ = []
for i in range(model.n_components):
    means_.append(model.means_[i][0])
    vars_.append(np.diag(model.covars_[i])[0])
    print("第{0}号隐藏状态".format(i))
    print("mean = ", model.means_[i])
    print("var = ", np.diag(model.covars_[i]))
    print()

Means and vars of each hidden state
第0号隐藏状态
mean =  [-6.63217241e-03  5.88893426e+04]
var =  [4.84610308e-02 4.85498280e+08]

第1号隐藏状态
mean =  [-7.25238348e-02  4.22133974e+05]
var =  [1.10870852e+00 1.09613199e+11]

第2号隐藏状态
mean =  [4.41582696e-02 1.42342273e+05]
var =  [1.42504830e-01 3.27804399e+09]

因为可见层输入采用了两个输入特征（涨跌幅、成交量），所以每个结点有两个元素，分别代表该状态的涨跌幅均值、成交量均值。

由于系统的目标预测涨跌幅，所以这里只关心第一个特征的均值，有如下结论：

状态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，是三个状态中方差最大的，即该状态的变化范围较大，并非很可信。

print("Transition matrix")
print(model.transmat_)

Transition matrix
[[0.82001546 0.01739021 0.16259432]
 [0.03793618 0.22931597 0.73274785]
 [0.23556112 0.04332874 0.72111014]]

第一行最大的数值是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.0066321724101791124, -0.07252383481713848, 0.04415826960924146]

vars_

[0.048461030777300955, 1.1087085156555492, 0.14250482994747]

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.0066321724101791124, 0.04415826960924146, 0.04415826960924146]
predict_vars [0.048461030777300955, 0.14250482994747, 0.14250482994747]

# 画图
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 == i
    yesterday_mask = np.concatenate(([False], mask[:-1]))
    if len(dates[mask]) <= 0:
        continue

    if 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 day
    ax.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）。

参考资料

天青如水

关注

2
点赞
踩
7

收藏

觉得还不错? 一键收藏
2
评论
HMM

HMM公式推导在 HMM 中，有两个基本假设：齐次 Markov 假设（未来只依赖于当前）：p(it+1∣it,it−1,⋯ ,i1,ot,ot−1,⋯ ,o1)=p(it+1∣it)p(i_{t+1}|i_t,i_{t-1},\cdots,i_1,o_t,o_{t-1},\cdots,o_1)=p(i_{t+1}|i_t)p(it+1∣it,it−1,⋯,i1,ot,o...
复制链接

扫一扫