在高斯隐马尔可夫模型中解码序列

Question

我正在研究用隐马尔可夫模型解决股市预测问题。我的数据矩阵包含了特定证券的各种特征：

01-01-2001, .025, .012, .01
01-02-2001, -.005, -.023, .02

我拟合了一个简单的高斯HMM模型：

from hmmlearn import GaussianHMM
mdl = GaussianHMM(n_components=3,covariance_type='diag',n_iter=1000)
mdl.fit(train[:,1:])

使用模型(λ)，我可以解码观察向量，找到与观察向量对应的最可能的隐藏状态序列。

print mdl.decode(test[0:5,1:])
(72.75, array([2, 1, 2, 0, 0]))

在上面，我已经解码了一个观测向量O_t = (O₁, O₂, ..., O_d)的隐藏状态序列，其中包含测试集中前五个实例。我想要估计测试集中第六个实例的隐藏状态。思路是迭代地遍历第六个实例可能的离散特征值集合，并选择具有最高似然度argmax = P(O₁, O₂, ..., O_d+1 | λ )的观测序列O_t+1。一旦我们观察到O_d+1的真实特征值，我们可以将长度为5的序列向后移动一个位置，然后再次进行所有操作：

    l = 5
    for i in xrange(len(test)-l):
        values = []
        for a in arange(-0.05,0.05,.01):
            for b in arange(-0.05,0.05,.01):
                for c in arange(-0.05,0.05,.01):
                    values.append(mdl.decode(vstack((test[i:i+l,1:],array([a,b,c])))))
     print max(enumerate(values),key=lambda x: x[1])

问题在于，当我解码观察向量O_t+1时，最高可能性的预测几乎总是相同的（例如，最高可能性的估计总是具有O_d+1特征值等于[0.04 0.04 0.04]和隐藏状态[0]）：

(555, (74.71248518927949, array([2, 1, 2, 0, 0, 0]))) [ 0.04  0.04  0.04]
(555, (69.41963358191555, array([2, 2, 0, 0, 0, 0]))) [ 0.04  0.04  0.04]
(555, (77.11516871816922, array([2, 0, 0, 0, 0, 0]))) [ 0.04  0.04  0.04]

我完全有可能误解了mdl.decode的用途，因此使用不正确。如果是这种情况，我该如何最好地迭代可能的O_d+1值，然后最大化P(O₁, O₂, ..., O_d+1 | λ)？

- datasci

---

你能否提供完整的代码，而不仅仅是片段呢？这将使问题更加清晰明了。 - Jaffer Wilson

如果您在每次迭代中重置“values”，则最终的“print”语句将仅引用最后一个循环中的“values”。 - andrew_reece

JafferWilson：好的，我会做。@andrew_reece：我认为这是预期的行为。“i”代表了5个观察序列test [i:i+l]的开头。想法是对于第六个（未被观察到的）实例[a，b，c]的可能特征值进行离散迭代，使用包括第六个实例估计的新序列进行解码，并选择具有最高概率的最终序列。完成后，我们重置value对象，将序列位置向前移动一位，并估计下一个观测集。也许我只需要用前向-后向算法来做到这一点... - datasci

Answer 1

你的实际值是否被限制在 [-0.05, 0.05) 范围内？

当我最初构建一个样本数据集来研究你的问题时，我在 [0,1] 中随机生成浮点数。当我这样做时，我也得到了你观察到的相同结果 - 对于每个序列的第六个条目，总是选择 (a,b,c) 的最大值，并且总是相同的预测类别。但考虑到我的数据分布（均匀分布在 0 和 1 之间）比第六个条目的网格搜索值（在 -.05 和 .05 之间）具有更大的中心倾向性，因此 HMM 始终选择最高值的三元组 (.04,.04,.04)，因为它最接近其训练的分布的主密度。

当我使用与我们允许第六个条目相同范围内的可能值的均匀分布进行抽取时，输出变化更多：每个序列的 O_t+1 选择和类别预测都显示出合理的差异。从你的示例数据中，似乎你至少拥有正值和负值，但你可以尝试绘制每个特征的分布，并查看你可能的第六个条目值的范围是否都是合理的。

这里有一些样本数据和评估代码。每次出现新的最优 (a,b,c) 序列或第六项观测的预测发生变化时，它都会打印出一条消息（只是为了表明它们并不全都相同）。每个 6 元素序列的最高可能性，以及预测和最佳的第六个数据点存储在 best_per_span 中。

首先，构建一个样本数据集：

import numpy as np
import pandas as pd

dates = pd.date_range(start="01-01-2001", end="31-12-2001", freq='D')
n_obs = len(dates)
n_feat = 3
values = np.random.uniform(-.05, .05, size=n_obs*n_feat).reshape((n_obs,n_feat))
df = pd.DataFrame(values, index=dates)

df.head()
                   0         1         2
2001-01-01  0.020891 -0.048750 -0.027131
2001-01-02  0.013571 -0.011283  0.041322
2001-01-03 -0.008102  0.034088 -0.029202
2001-01-04 -0.019666 -0.005705 -0.003531
2001-01-05 -0.000238 -0.039251  0.029307

现在将数据集分成训练集和测试集：

train_pct = 0.7
train_size = round(train_pct*n_obs)
train_ix = np.random.choice(range(n_obs), size=train_size, replace=False)
train_dates = df.index[train_ix]

train = df.loc[train_dates]
test = df.loc[~df.index.isin(train_dates)]

train.shape # (255, 3)
test.shape # (110, 3)

在训练数据上拟合3状态HMM：

# hmm throws a lot of deprecation warnings, we'll suppress them.
import warnings
with warnings.catch_warnings():
    warnings.filterwarnings("ignore",category=DeprecationWarning)
    # in the most recent hmmlearn we can't import GaussianHMM directly anymore.
    from hmmlearn import hmm

mdl = hmm.GaussianHMM(n_components=3, covariance_type='diag', n_iter=1000)
mdl.fit(train)

现在进行网格搜索以找到最佳的第六个(t+1)观测值：

# length of O_t
span = 5

# final store of optimal configurations per O_t+1 sequence
best_per_span = []

# update these to demonstrate heterogenous outcomes
current_abc = None
current_pred = None

for start in range(len(test)-span):
    flag = False
    end = start + span
    first_five = test.iloc[start:end].values
    output = []
    for a in np.arange(-0.05,0.05,.01):
        for b in np.arange(-0.05,0.05,.01):
            for c in np.arange(-0.05,0.05,.01):
                sixth = np.array([a, b, c])[:, np.newaxis].T
                all_six = np.append(first_five, sixth, axis=0)
                output.append((mdl.decode(all_six), (a,b,c)))

    best = max(output, key=lambda x: x[0][0])

    best_dict = {"start":start,
                 "end":end,
                 "sixth":best[1],
                 "preds":best[0][1],
                 "lik":best[0][0]}  
    best_per_span.append(best_dict)

    # below here is all reporting
    if best_dict["sixth"] != current_abc:
        current_abc = best_dict["sixth"]
        flag = True
        print("New abc for range {}:{} = {}".format(start, end, current_abc))

    if best_dict["preds"][-1] != current_pred:
        current_pred = best_dict["preds"][-1]
        flag = True
        print("New pred for 6th position: {}".format(current_pred))

    if flag:
        print("Test sequence StartIx: {}, EndIx: {}".format(start, end))
        print("Best 6th value: {}".format(best_dict["sixth"]))
        print("Predicted hidden state sequence: {}".format(best_dict["preds"]))
        print("Likelihood: {}\n".format(best_dict["nLL"]))

循环运行时的报告输出示例：

New abc for range 3:8 = [-0.01, 0.01, 0.0]
New pred for 6th position: 1
Test sequence StartIx: 3, EndIx: 8
Best 6th value: [-0.01, 0.01, 0.0]
Predicted hidden state sequence: [0 2 2 1 0 1]
Likelihood: 35.30144407374163

New abc for range 18:23 = [-0.01, -0.01, -0.01]
New pred for 6th position: 2
Test sequence StartIx: 18, EndIx: 23
Best 6th value: [-0.01, -0.01, -0.01]
Predicted hidden state sequence: [0 0 0 1 2 2]
Likelihood: 34.31813078939214

best_per_span 的示例输出如下：

[{'end': 5,
  'lik': 33.791537281734904,
  'preds': array([0, 2, 0, 1, 2, 2]),
  'sixth': [-0.01, -0.01, -0.01],
  'start': 0},
 {'end': 6,
  'lik': 33.28967307589143,
  'preds': array([0, 0, 1, 2, 2, 2]),
  'sixth': [-0.01, -0.01, -0.01],
  'start': 1},
 {'end': 7,
  'lik': 34.446813870838156,
  'preds': array([0, 1, 2, 2, 2, 2]),
  'sixth': [-0.01, -0.01, -0.01],
  'start': 2}]

除了报告元素外，这并不是对您最初方法的重大更改，但它似乎按预期工作，而不会每次都达到最大值。

- andrew_reece