★ FEATURED ARTICLE 做时序数据分析的时候很多时候表面看到的数据背后其实藏着一层看不见的状态在驱动。比如语音信号里不同的音素交替出现比如股票价格背后牛熊切换再比如用户行为序列里那些没被直接观察到的意图变化。这种“观测可见、状态隐藏”的序列建模问题隐马尔可夫模型HMM是经典解法而 hmmlearn 就是 Python 里训练 HMM 最省心的库。这篇文章我直接讲实操从模型选型到数据格式从训练细节到调参避坑把 hmmlearn 跑通 HMM 的完整链路拆开揉碎。1. 为什么是 HMM为什么选 hmmlearn1.1 HMM 在解决什么问题先对齐一下概念。HMM 假设你看到的数据序列是由一条看不见的“状态链”生成的。每个时刻系统处于某个离散状态状态按照转移概率矩阵跳转每个状态下又会按照发射概率产生一个可观测的输出。我们的目标就是根据观测序列反推出状态序列或者估计出模型参数让模型尽可能好地解释数据。拿一个我实际做过的场景举例。某套设备的状态监测系统里传感器返回的是连续振动幅值但设备本身其实只有“正常”“磨损”“即将故障”三种离散工况。振动幅值是连续的观测工况就是隐藏状态。HMM 恰好能同时完成两件事一是估计出三种工况各自对应的振动分布二是根据一段振动序列判断当前最可能处于哪种工况。这类问题用 hmmlearn 来处理代码量非常小而且不需要你自己推导 EM 算法。1.2 hmmlearn 的定位与优势hmmlearn 最早是从 scikit-learn 里拆分出来的独立库专门做隐马尔可夫模型在 PyPI 上可以直接安装。它不是最通用的概率编程框架比如 PyMC 那种而是聚焦 HMM 这一件事开箱即用。我推荐它的核心原因有三个第一API 设计非常接近 scikit-learn 的习惯fit()、predict()、score()一脉相承熟悉 sklearn 的人几乎零学习成本。第二内置了三种常用的发射分布模型GaussianHMM适合观测值是连续向量的场景MultinomialHMM适合观测值是离散符号的场景GMMHMM在连续观测但单高斯拟合不够时可以上。第三模型训练用的是经典的 Baum-Welch 算法本质是 EM 算法前向-后向计算和参数更新全都封装好了你只需要把数据喂进去。安装也很简单pip install hmmlearn依赖的 numpy、scikit-learn 等库会自动带上。如果是在 conda 环境里直接conda install -c conda-forge hmmlearn也行。2. 动手前先选型三种模型怎么挑2.1 GaussianHMM连续观测值的主力大多数实际问题里观测值都不是简单的离散符号而是连续数值甚至是多维向量。比如语音的 MFCC 特征、股票收益率、传感器读数。GaussianHMM假设每个隐藏状态下的观测向量服从一个多元高斯分布每个状态都有自己的均值向量和协方差矩阵。它在 hmmlearn 里用的参数是covariance_type有四种选择spherical、diag、full、tied。我用一个表格把这几种协方差类型的差别和适用场景列清楚类型参数数量适用场景备注spherical每个状态一个方差标量特征各维度独立且方差相等参数最少训练最快但假设太强diag每个状态一个对角协方差特征各维度独立但方差不同最常用平衡了自由度和稳定性full每个状态一个完整协方差矩阵特征维度之间存在相关性参数多要求数据量大否则容易过拟合tied所有状态共享同一个完整协方差矩阵状态差异主要体现在均值上用协方差刻画全局噪声水平常见于信号处理我自己在项目里最常用的是diag因为多维特征之间通常没有那么强的耦合关系强行用full在小数据集上反而会把协方差矩阵估计得奇奇怪怪。2.2 MultinomialHMM离散观测值的场景如果观测值是离散符号比如 NLP 里的词性序列、DNA 序列里的碱基、用户点击行为里的页面编号那就用MultinomialHMM。它每个状态下对应一个类别分布也就是每个符号在该状态下出现的概率。需要注意的是MultinomialHMM 的输入数据必须是非负整数编码。也就是说你要先把原始符号映射成0, 1, 2, ...这样的索引否则模型会直接报错。2.3 GMMHMM连续观测值更复杂的分布有些场景下单个高斯分布描述不清楚状态内的观测分布。比如语音识别中同一个音素在不同人的发音下可能会呈现双峰或多峰分布单高斯拟合效果不好。GMMHMM在每个状态下用混合高斯模型GMM来建模发射概率每个状态内部可以有多个高斯成分。它的额外参数是n_mix表示每个状态用几个高斯成分。代价是参数数量暴涨训练时间明显变长而且更容易陷入局部最优。我的建议是先用GaussianHMM跑通如果score明显偏低或者业务上有理由相信分布是多峰的再换GMMHMM。3. 训练 HMM 的核心细节从数据准备到 fit3.1 数据格式最重要也最容易翻车的地方hmmlearn 对训练数据的格式要求非常严格这一点新手几乎必踩。fit()方法接受的X是一个二维数组形状是[n_samples, n_features]其中n_samples是所有观测时间步的总数n_features是每一步观测向量的维度。关键在于如果你有多个独立的观测序列不能简单地把它们拼接成一个大数组那样会让模型把上一个序列的结尾和下一个序列的开头连在一起学到根本不存在的转移关系。正确的做法是用numpy.nan作为序列之间的分隔符让模型知道每个序列在哪里结束、下一条从哪里开始。直接看代码import numpy as np from hmmlearn import hmm # 模拟两个独立序列每个序列 100 个时间步每步 2 维特征 seq1 np.random.randn(100, 2) seq2 np.random.randn(100, 2) # 中间用 nan 行分隔 X np.vstack([seq1, np.full((1, 2), np.nan), seq2]) # 对应的序列长度数组hmmlearn 0.3.0 之后也支持这种方式 lengths [100, 100]关于lengths参数我想多说一句。hmmlearn 在不同版本里对序列长度的处理方式略有变化。旧版本要求你必须传lengths新版本开始支持NaN分隔后自动推断。我建议在实际项目中两者都写上兼容性最好也方便别人阅读你的代码。如果你只有一个长序列那就不需要NaN分隔直接传二维数组即可。对于MultinomialHMM数据格式要求是二维整数数组每个元素是一个离散符号的整数索引。同样多个序列之间用-1作为分隔符注意不是NaN因为整数数组里NaN会被强制转成-1所以规范做法就是直接放-1。3.2 初始化模型的关键参数创建模型对象时最核心的参数是n_components也就是隐藏状态的数量。这个参数你不能让程序自动选必须你自己根据业务去定。怎么定三个思路一是领域知识。设备监测就知道有几种工况语音就知道有几十个音素直接拿这个数去初始化。二是数据探索。把观测值做聚类看看分成几类时轮廓系数最高用这个数作为n_components的参考。三是模型比较。分别用不同的n_components训练模型比较score的结果选对数似然提升不再明显的那个数。n_iter控制 EM 算法最大迭代次数默认是 10我在项目里通常设到 50 到 100避免迭代次数不够导致收敛不充分。tol是对数似然增益的容忍度默认 1e-2如果你的数据量很大可以调小到 1e-4让模型收敛得更精细。还有一个经常被人忽略的初始化细节hmmlearn会自动初始化转移矩阵和发射参数。如果你没有指定startprob_prior、transmat_prior之类的先验参数它用的是均匀分布加一个小的伪计数。对于大多数任务问题不大但如果你对初始状态分布有先验知识可以手动塞进去。model hmm.GaussianHMM( n_components3, covariance_typediag, n_iter50, tol1e-4, random_state42, )random_state我强烈建议固定。因为 EM 算法对随机初始化敏感不同种子可能收敛到不同的局部最优解固定种子能保证结果可复现。3.3 训练fit 之后看什么调用model.fit(X, lengths)之后训练就完成了。但这里有一个关键点hmmlearn 的fit()在 0.3.0 版本之前会自动用KMeans对观测数据聚类来初始化发射概率的均值导致训练结果很大程度上依赖初始聚类点。如果你发现两次训练结果差异很大可以先检查版本再用固定random_state压住随机性。训练完成后检查模型学到的参数是验证训练是否合理的一步。print(初始状态概率:, model.startprob_) print(转移概率矩阵:\n, model.transmat_) print(各状态均值:\n, model.means_) print(各状态协方差:\n, model.covars_)一个有业务意义的模型学出来的转移矩阵对角线元素通常会偏大因为状态本身有惯性不太会疯狂跳来跳去。如果对角线特别小说明数据里状态切换非常频繁要么状态数量设多了要么观测噪声太大。各状态的均值应该明显分得开如果两个状态的均值向量几乎一样说明其中一个状态是多余的。4. 实操一个完整的训练流程4.1 准备模拟数据为了把流程走通我构造一个简单的模拟场景三状态系统每个状态产生不同均值的二维高斯观测。两个状态有清晰区分第三个状态和第一个有轻微重叠这样可以检验模型的判别能力。import numpy as np import matplotlib.pyplot as plt from hmmlearn import hmm # 定义真实状态数和观测分布 true_means np.array([ [-2.0, 0.0], [2.0, 0.0], [0.0, 4.0], ]) covars np.array([[[0.5, 0.4], [0.4, 0.5]]] * 3) # 生成两个独立序列 lengths [200, 300] obs [] for T in lengths: states [] seq_obs [] state 0 for t in range(T): # 模拟一个简单的状态转移大部分时间留在当前状态偶尔切换 if np.random.rand() 0.03: state np.random.choice(3) states.append(state) seq_obs.append(np.random.multivariate_normal(true_means[state], covars[state])) obs.append(np.array(seq_obs)) X np.vstack([obs[0], np.full((1, 2), np.nan), obs[1]])注意这里我在两个序列之间插入了nan行这是 hmmlearn 组合多条序列的标准姿势。4.2 训练并观察收敛情况model hmm.GaussianHMM( n_components3, covariance_typefull, n_iter100, tol1e-4, random_state42, ) model.fit(X, lengths)训练结束后先看模型是否收敛。可以打印model.monitor_查看收敛信息print(收敛是否成功:, model.monitor_.converged) print(迭代次数:, model.monitor_.iter) print(历史对数似然增益:, model.monitor_.history[-5:])如果converged是False说明 100 次迭代内没有达到tol的收敛条件。这不是报错但通常意味着你需要增加n_iter或者数据本身的似然增长已经非常平缓。对于大数据集对数似然的增长曲线后期会非常缓慢tol设成默认值可能永远不收敛这时候调大tol到 1e-2 反而更实际因为继续迭代带来的提升已经可以忽略。4.3 解码预测每个时刻的状态predict用的是维特比算法也就是给定观测序列和模型参数求最可能的一条状态路径。states_pred model.predict(X)这里X是包含nan分隔行的完整数组。predict返回的是一个长度等于X.shape[0]的一维数组每个位置对应一个时刻的预测状态。需要注意的是nan分隔行对应的预测状态是需要忽略的实际处理中可画图时直接跳过。对比预测状态和真实状态如果模拟数据里有真实状态标签的时候要注意 HMM 的状态编号和真实状态编号不一定一一对应。模型学到的状态 0、1、2 和真实的状态 0、1、2 可能经过了一次置换。所以计算准确率之前先做一个状态映射或者直接算聚类一致性指标比如 ARI更科学。4.4 评估模型score 的用法和误区score返回的是对数似然log P(X | model)也就是模型解释这份观测数据有多好。注意这个值是负数而且会随着序列长度变大而变大绝对值更大所以不能拿不同长度序列的score直接比较。ll model.score(X, lengths) print(对数似然:, ll)如果你要用score来选择模型比如比较不同的n_components应该在同一份数据上做交叉验证或者留一序列评估保证比较的公平性。还有一种方式是计算 BIC 或 AIC把模型复杂度惩罚进去hmmlearn 没有内置这个功能但你可以自己算n_params model.n_features * model.n_components * 2 model.n_components * model.n_components bic -2 * ll n_params * np.log(np.sum(lengths))这里n_params的估算比较粗略严格来说还要加上协方差矩阵的参数数量但对于模型选择的参考足够了。5. 常见报错与避坑实录5.1 “Input array must be 2-dimensional”这个报错几乎是新手第一课。hmmlearn 不接受一维数组作为观测序列输入。就算你每个时间步只有一个维度的观测值X也必须是二维的形状为[n_samples, 1]。解决办法是加一个reshape(-1, 1)。X X_1d.reshape(-1, 1)5.2 出现 NaN 分隔后 predict 结果出现奇怪的负值或大整数如果用的是MultinomialHMM序列之间用-1分隔后predict结果里可能会出现一个额外状态值用来对应分隔符的位置这个值本身没有业务含义。过滤掉即可valid_idx np.all(np.isfinite(X), axis1) states_clean states_pred[valid_idx]5.3 训练时直接崩掉报协方差矩阵奇异这个问题的核心是某个状态下的协方差矩阵变成了半正定甚至奇异矩阵通常是数据量不足、特征维度太高、或者某个状态几乎没有数据点分配到。解决思路有几个第一改用covariance_typediag减少协方差矩阵的自由参数。第二检查是否有某个状态在迭代中逐渐失去了所有数据点。如果是把n_components调小或者手动调整初始均值让各状态的初始覆盖更均匀。第三如果数据确实太少可以调大min_covar参数默认是 1e-3给协方差矩阵的对角线加一个小的正则项。5.4 随机种子带来的结果不可复现我把random_state当成项目中不可省略的参数不是因为它能提高精度而是因为它能让你在调参过程中拿到可复现的结果。否则同样的代码跑两遍训练出来的模型参数不同排查问题时非常气人。提交代码或者写文档时一定要写明用的是哪个种子。5.5 状态数量选多了很多同学觉得n_components越大越好实际上不是。状态多了之后模型容易把一些本来属于同一状态的观测强行拆成几个状态导致状态之间区分度低而且转移矩阵变得稀疏甚至出现毫无意义的转移。我通常会用 BIC 来辅助判断状态数在 BIC 曲线的最低点附近就是比较合理的选择。当然最终还要看业务解释性——你分出来的状态能不能对上实际业务含义对不上的状态再多也是噪声。6. 我自己用 hmmlearn 的几个经验这套流程我在几个项目里都跑过说几个可能是文档里不会明说、但非常影响体验的细节。一个是fit的输入数据最好提前做标准化或者归一化。hmmlearn 的单高斯发射模型对特征的尺度很敏感如果某个特征的值域是 0 到 1另一个特征的值域是几千到几万协方差矩阵会很难估计score的提升也会很别扭。我一般先做StandardScaler或者MinMaxScaler这样训练稳定性和收敛速度都会好很多。另一个是多序列训练时的lengths处理。我踩过一次坑在官方文档和 GitHub issue 里翻到过好几种说法最后确认下来lengths必须是一维整数数组每个元素是每个序列的步数而且所有长度之和必须等于X中有效样本数不含nan行。如果你用的是NaN分隔那么lengths可以省略但要小心版本兼容性。还有一个经验是关于预测结果的可视化。我习惯在训练完成后把predict得到的状态序列画成阶梯图叠加在原始观测序列上。这样一眼就能看出状态切分是否符合业务直觉。很多时候模型训练指标不错但画出来才发现状态切换非常频繁、毛刺特别多这通常说明观测噪声太大或者状态数不合适比单纯看score直观得多。最后补充一句hmmlearn 对于大多数中小规模序列建模任务已经完全够用但如果你的序列特别长几十万步以上或者状态数特别多几十个甚至上百个训练速度会成为一个瓶颈。这时候可以考虑用GaussianHMM加spherical协方差跑一版快速验证再换回更精细的模型调优不要一上来就上最复杂的配置。
阅读完成 · 觉得有帮助？
