尧图网站建设 
简介一套基于MATLAB的HMM隐马尔科夫模型时间序列预测完整实现方案面向时序预测学习者、信号处理及数据科学从业者解决如何用隐马尔科夫模型对一维观测序列进行训练、解码与未来值预测的问题。资源包含3个MATLAB脚本、3张结果示意图、1份说明文档和1个数据文件共8个文件压缩包大小为238KB轻量但环节完整。脚本覆盖参数训练、前向算法、Viterbi解码与Baum-Welch参数优化并提供MSE、MAE等评估指标的可视化对比内置数据文件包含训练与测试时序数据示意图直观展示模型结构与预测效果文档梳理完整实现步骤。目前已有3414人学习下载。通过运行这些代码可快速复现完整的时序预测流程掌握状态转移矩阵与观测发射矩阵的建模细节为后续改进模型或迁移到语音、金融等序列预测任务打下基础。1. 用HMM做时间序列预测先搞清楚它在预测什么先把结论放这儿HMM隐马尔科夫模型做时间序列预测预测的不是“下一时刻的数值”而是“下一时刻最可能处于哪个隐含状态以及在这个状态下观测值的分布”。这个区别决定了你能不能把它用对。这套MATLAB实现包含完整源码和示例数据覆盖数据预处理、HMM训练、隐状态解码、逐步预测和RMSE/MAE评估适合数据量不大但存在明显模式切换的时间序列——比如用电负荷的峰谷切换、风速的静风与大风交替、金融序列的震荡与趋势切换。新手可以照脚本跑通全流程熟手可以直接替换成自己的数据调整状态数和离散化参数。2. HMM的建模逻辑五个要素、三个问题、一次映射2.1 五元组状态、观测与三个概率矩阵HMM的数学结构是五元组状态集合S{s₁,...,s_N}观测集合O{o₁,...,o_M}连续观测则对应分布状态转移矩阵AN×NA(i,j)P(s_{t1}j|s_ti)观测概率矩阵BN×MB(j,k)P(o_tk|s_tj)连续观测时是高斯的均值与方差初始状态分布π1×N。映射到时间序列预测隐状态是“看不见的模式”。比如用电负荷数据里的“工作日早高峰状态”“午间平谷状态”“晚间高峰状态”观测是“实际测到的负荷值”。HMM的核心假设有两条t1时刻的状态只依赖t时刻的状态一阶马尔可夫t时刻的观测只依赖t时刻的状态条件独立。工程上这两个假设往往不严格成立但好处是模型的参数极少、训练极快、可解释性极强。你付出的代价是模型对长程依赖和复杂非线性关系的表达能力不如LSTM这类深度模型。2.2 三个基本问题在预测里的角色HMM有三个经典问题它们分别对应预测流程的不同环节。评估问题给定模型λ(A,B,π)和观测序列O计算P(O|λ)。这其实是模型选择的核心工具。你要决定状态数取3还是5不能靠拍脑袋而是分别训练模型、比较对数似然再用AIC/BIC做权衡。解码问题给定模型和观测求最可能的隐状态序列。这一步回答“当前处于哪个模式”Viterbi算法就是干这个的。学习问题给定观测序列估计模型参数A、B、πMATLAB里的hmmtrain就是Baum-Welch算法EM的实现。基本问题算法MATLAB函数在预测中的用途评估前向算法hmmdecode输出LOGPSEQ模型选择比较不同状态数下的对数似然解码Viterbihmmviterbi获得最可能的隐状态序列用于状态解释学习Baum-WelchEMhmmtrain估计转移矩阵A和观测矩阵B预测时三个问题按顺序上场先用学习问题在训练集上学出A和B再用解码或前向滤波得到当前状态分布最后用A和B推下一时刻观测的概率分布取期望作为点预测。这就是HMM做预测的完整链路每一步都有对应的MATLAB函数这也是这套源码能写得这么短的原因。2.3 离散观测还是连续观测MATLAB工具箱的边界MATLAB的hmmtrain、hmmviterbi、hmmdecode在设计之初面向的是离散观测符号比如DNA序列的ACGT、语音识别里的码本索引。也就是说B矩阵是一个N×M的离散概率表观测序列必须是1、2、3这样的正整数索引。这就带来一个直接后果原始的时间序列是连续数值喂给hmmtrain之前必须离散化。常见做法是分桶把训练集的观测值按等频或等距分成K个桶每个值映射成一个符号。代价是丢失桶内的数值精度点预测只能回到桶中心或桶均值预测曲线呈阶梯状。另一种做法是自己写高斯发射HMM的EM更新让每个状态的观测分布是N(μ_j, σ_j²)把B矩阵换成高斯密度。这样不需要离散化预测更平滑还能输出平滑的置信区间。但这需要自己实现前向、后向和EM的M步更新代码量大约多一百行。资源包以离散版本为主推方案因为hmmtrain足够稳定不容易在数值上翻车。如果你想用连续版本第6章会给一个改造思路。2.4 同门对比HMM vs ARIMA vs LSTM做时间序列预测的人通常先想到ARIMA或LSTM。ARIMA是线性模型擅长捕捉自相关结构但遇到状态突变比如负荷从平静状态突然跳入高峰状态时反应迟钝预测曲线会出现明显的滞后。LSTM在数据量足够大时表现很好但要调的参数多、可解释性差小样本下容易过拟合——这是“lstm时间序列预测python”这类问题下最常见的坑。HMM恰好站在两者中间它对几千个点的小样本友好EM迭代通常几十步就收敛它给出显式的状态语义业务上能解释“当前在哪个模式”它本质上是概率模型能输出预测区间而不只是点预测。代价是一阶马尔可夫假设偏弱对长程依赖捕捉有限。我给的建议是数据存在明显模式切换、样本量在几百到几千、需要向业务解释预测依据时用HMM你只有海量数据且只追求精度时LSTM更合适。资源包里的示例数据就是典型的模式切换序列用HMM的预测滞后通常比ARIMA小一个量级。3. 源码拆解从数据到预测值的完整链路3.1 数据加载、差分与观测离散化第一步是加载数据并做预处理。资源包里的hmm_data.mat存了一个变量series是单变量时间序列列向量带有明显的模式切换特征。换成你自己的数据时只要保证series是列向量即可。diff做一阶差分是很多时序预测的常规操作原始负荷值非平稳HMM学“变化量”的规律比学“绝对数值水平”更稳健。差分后归一化到[0,1]再用等频分桶把连续值映射成离散符号。%% 3.1 数据加载与观测离散化 clear; clc; rng(42); % 固定随机种子保证结果可复现 load(hmm_data.mat); % 变量名 series列向量 % 差分平稳化HMM学的是变化量的规律不是绝对水平 dy diff(series); y_norm (dy - min(dy)) / (max(dy) - min(dy)); % 映射到[0,1] % 等频分桶让每个桶的样本量均衡避免某些符号样本过少 K 10; % 桶数K10是经验值第4章讲怎么调 edges prctile(y_norm, linspace(0, 100, K 1)); % K1个边界 edges(end) edges(end) eps; % 防止最大值落到桶外 [~, obs] histc(y_norm, edges); obs(obs 0) 1; % histc在等于下边界时可能返回0强制归入1号桶这段代码做了四件事差分、归一化、分桶、符号映射。关键参数K直接决定B矩阵的规模N×K也决定预测精度。K太小预测退化成几个台阶K太大每个桶样本量不足、B矩阵出现大量近似零概率、训练过拟合。我一般取10到30具体看样本量。这里有个细节值得注意edges(end) edges(end) eps这行不是可有可无的。当某个观测值恰好等于最大值时histc会把它归到最后一个桶之外返回0导致obs出现0这个非法符号。加一个eps把上边界稍微抬高保证所有值都落在桶内。3.2 训练初始化、hmmtrain与Viterbi解码训练的核心是hmmtrain它的输入是观测符号序列、初始转移矩阵、初始观测矩阵输出训练后的A和B。初始值的给法直接影响是否收敛到合理解这里先给出能跑通的版本。%% 3.2 HMM训练与状态解码 trainLen round(0.8 * length(obs)); % 前80%训练后20%测试 tr_obs obs(1:trainLen); te_obs obs(trainLen1:end); nStates 3; % 隐状态数负荷场景下解释为 平/升/降 三个模式 A0 ones(nStates) / nStates; % 均匀初始化确保没有全零行 A0 A0 0.2 * eye(nStates); % 对角线加一点权重倾向状态驻留 A0 A0 ./ sum(A0, 2); % 行归一化 B0 ones(nStates, K) / K; % 观测矩阵均匀初始化 % 训练最大迭代200次容差1e-4 [At, Bt] hmmtrain(tr_obs, A0, B0, MaxIterations, 200, Tolerance, 1e-4); % 用Viterbi解码训练集的最可能状态序列 states hmmviterbi(tr_obs, At, Bt); % 评估训练集对数似然用于模型选择 [~, logp_train] hmmdecode(tr_obs, At, Bt); fprintf(训练完成log-likelihood %.4f\n, logp_train);逻辑说明A0 A0 0.2 * eye(nStates)这一步很多人会忽略但它很关键。如果初始A是全均匀矩阵EM迭代容易陷入对称解——三个状态学出来一模一样。对角线加权重相当于先验地认为“状态倾向于维持”这符合大多数物理时序的惯性特征。参数说明MaxIterations控制最大迭代次数数据量小时几十步就收敛200步是安全上限。Tolerance控制对数似然增幅阈值1e-4够用设成1e-8反而容易在数值底部来回震荡。注意旧版MATLAB2018a之前的参数名是MaxIter新版是MaxIterations如果报错先检查这里。hmmtrain不输出初始分布π它内部按均匀分布处理所以这里不需要初始化π。3.3 预测前向概率加权而不是只看最可能状态预测阶段最常见的错误是只用Viterbi解码出的单一状态去查转移矩阵。正确做法是用滤波分布加权因为Viterbi只给你一条最可能路径但当前状态本质上是一个分布尤其在状态边界附近两个状态的概率可能各占一半。%% 3.3 时间序列预测一步预测前向概率加权 centers (edges(1:K) edges(2:K1)) / 2; % 每个桶的中心值 pred_dy zeros(length(te_obs), 1); for i 1:length(te_obs) % 只使用到 t 时刻为止的真实观测避免未来信息泄露 cur_obs [tr_obs; te_obs(1:i-1)]; % i1时就是训练集 [~, ~, forward] hmmdecode(cur_obs, At, Bt); alpha_t forward(:, end); % 滤波分布1 x nStates next_obs_dist alpha_t * At * Bt; % 下一时刻观测分布1 x K pred_dy(i) next_obs_dist * centers; % 取期望作为点预测 end % 差分还原预测的是差分序列要累积回去才是真实值 pred_raw cumsum([series(trainLen 1); pred_dy]);逻辑说明hmmdecode的第三个输出forward是前向概率矩阵forward(:, end)就是截至当前时刻的滤波分布它比Viterbi的硬解码更稳健。next_obs_dist alpha_t * At * Bt这一步的本质是先把当前状态分布用转移矩阵传播到下一时刻再对每个可能状态求观测分布最后按状态概率加权得到的就是P(o_{t1}|o_1...o_t)。点预测取期望而不是众数这样预测曲线更平滑——众数只能落在桶中心期望是加权平均等于做了软输出。cumsum还原时以series(trainLen1)为基准这个基准点选择非常关键错了整条预测曲线就平移一位。提示hmmdecode在循环里每步都在重算整个前缀序列复杂度O(t·N²)。序列超过一万点时会明显变慢可以改成在线前向更新只保留上一时刻的前向向量每一步做一次迭代前向递推效果完全一样。3.4 评估RMSE、MAE、MAPE与可视化预测做完一定要评估不能只看拟合曲线。这里同时算了差分尺度的RMSE/MAE和还原尺度的MAPE因为差分还原会累积误差只看还原后的误差无法判断模型本身在“变化量”上的预测能力。%% 3.4 评估指标 te_aligned dy(trainLen1:end); % 与pred_dy对齐的真实差分值 rmse_diff sqrt(mean((te_aligned - pred_dy).^2)); mae_diff mean(abs(te_aligned - pred_dy)); mape_raw mean(abs((series(trainLen1:end) - pred_raw) ./ ... (series(trainLen1:end) eps))) * 100; fprintf(差分RMSE: %.4f\n, rmse_diff); fprintf(差分MAE : %.4f\n, mae_diff); fprintf(还原MAPE: %.2f%%\n, mape_raw); figure; plot(series(trainLen1:end), b-, LineWidth, 1.2); hold on; plot(pred_raw, r--, LineWidth, 1.2); legend(真实值, HMM预测值); xlabel(时间步); ylabel(序列值);索引对齐是这里最容易出错的地方dy是diff(series)的结果长度比series少一个点dy(trainLen1:end)恰好与te_obs一一对应。pred_raw的还原基准是series(trainLen1)所以真实值对应series(trainLen1:end)长度也一致。我在实际项目中踩过这个坑——训练时RMSE挺好一画图整条曲线整体左移一个点就是diff和cumsum之间的索引没对齐。4. 数据准备与参数调优四个旋钮四组经验值4.1 平稳化什么时候差分什么时候不差HMM本身不要求数据平稳因为它建模的是状态转移加观测发射的联合分布理论上非平稳序列也能学。但实际效果会很差如果原始序列有明显上升趋势模型会把“上升”也当做一个稳定模式新数据一旦出现不同趋势预测立刻失效。常见做法是先做一阶差分把序列从“数值水平”变成“变化量”HMM学的是变化量的模式。我判断是否差分的方法先画ACF图如果自相关系数衰减很慢就先差分差分后再看是否在零附近波动。差分一次不够就差分两次但二次差分后观测符号分布会非常集中——大部分值落在零附近——这时候要减少K否则很多桶是空的。资源包里默认做一阶差分。如果你自己的数据本身是平稳的比如白噪声附近的残差序列把差分那行注释掉即可但归一化和分桶要基于原始值重新计算边界不能沿用差分版本的edges。4.2 状态数N怎么选对数似然与AIC/BIC状态数N是HMM里最关键的旋钮。N太小模型无法区分不同模式预测退化成均值N太大过拟合训练集似然很高测试集一塌糊涂。标准做法是遍历N从2到8对每个N训练并记录对数似然用BIC权衡。%% 4.2 状态数选择BIC nList 2:8; bic zeros(size(nList)); for idx 1:length(nList) ns nList(idx); A0 ones(ns) / ns 0.2 * eye(ns); A0 A0 ./ sum(A0, 2); B0 ones(ns, K) / K; [At_t, Bt_t] hmmtrain(obs, A0, B0, MaxIterations, 200, Tolerance, 1e-5); [~, lp] hmmdecode(obs, At_t, Bt_t); k ns * (ns - 1) ns * (K - 1); % 转移矩阵观测矩阵的自由参数 bic(idx) -2 * lp k * log(length(obs)); end figure; plot(nList, bic, o-); xlabel(状态数N); ylabel(BIC); [~, bestN] min(bic); fprintf(BIC最优状态数: %d\n, nList(bestN));逻辑说明自由参数数量按转移矩阵去掉行归一化约束后的N(N-1)个、观测矩阵去掉每行归一化约束后的N(K-1)个来计算。初始分布π的参数在样本量面前影响极小BIC里可以忽略。用BIC而不是AIC是因为HMM参数多AIC的惩罚项太轻几乎总会选到N8。实际使用中BIC给出的N只是参考。N3对应“平/升/降”三个模式业务上说得通N4如果解出来的第四个状态说不清是什么那即使BIC选了4我也会退回3。模型是拿来解释业务的不是拿来拟合似然的。4.3 桶数K与预测平滑度的关系K的选择常常被忽略因为hmmtrain只要求符号是正整数很多教程随便取一个数。但K直接影响B矩阵的规模和预测曲线的平滑度。K太大每个桶的样本量少B矩阵的估计方差大K太小预测只能在少数几个桶中心之间跳。我常用的经验公式是K 10 round(sqrt(length(series) / 1000))。2000个点取K115000个点取K15。如果你用期望值做点预测3.3节的方式K偏小带来的阶梯效应会明显减轻因为期望是加权平均不是众数。如果你改成众数预测K偏小的话预测曲线会非常难看。K和N还有一个交互效应K越大观测矩阵B的自由参数越多BIC越倾向选择较小的N。所以实际调参的顺序是先定K再遍历N。反过来先定N再调K你会发现K怎么调BIC都不稳定。4.4 初始值与随机性让训练结果可复现hmmtrain的Baum-Welch是EM算法收敛到局部最优解初始值不同结果可能不同。资源包里的脚本统一用rng(42)固定随机种子并且A0、B0按确定性方式初始化这样每次跑结果一致。遇到换一台机器结果变化很大的情况先检查是不是没固定种子。初始转移矩阵对角线加的权重0.2对收敛结果影响很大。权重越大模型越倾向状态驻留学出来的状态平均停留时间越长权重太小模型容易学出高频抖动的状态序列业务上很难解释。0.2是我在负荷类时序上的经验金融数据我一般加到0.5因为金融状态切换频率更低。有一个快速检查初始值是否合理的方法训练完打印转移矩阵At看对角线元素是否显著大于非对角线。如果学到的主对角线都在0.8以上说明你的数据确实存在强驻留模式模型结构是合理的。如果对角线全部小于0.3大概率是状态数选多了模型在拿状态拟合噪声。5. 避坑指南HMM时间序列预测的五个典型翻车点5.1 观测序列编码必须是从1开始的连续正整数现象直接拿连续数值喂给hmmtrain报错“Error using hmmtrain ... 观测值必须是从1开始的整数”或者不报错但训练结果全是NaN。原因hmmtrain把观测符号当成离散索引要求max(obs) size(B,2)且min(obs) 1。如果你把归一化后的浮点数直接当obs索引自然是非法的。解决严格按3.1节的分桶流程走一遍分桶后检查min(obs)和max(obs)。注意histc在观测值恰好等于第一个边界时返回0所以obs(obs 0) 1这行必须保留。还有一个隐蔽场景如果某个观测值在edges之外比如测试集的数值比训练集更大分桶会得到K1甚至更大的非法符号。稳妥做法是用训练集的min/max先做一次clip再把clip后的值分桶。5.2 初始转移矩阵存在全零行训练直接崩现象hmmtrain迭代一两步后报错或者返回的At里有整行NaN。原因EM的E步要计算前向概率的归一化如果A0的某一行全零前向递推中经过该状态的所有路径概率恒为0分母为0数值上出现NaN。解决初始化A0时必须行归一化且每行至少有一个正概率。ones(nStates)/nStates 0.2*eye(nStates)这个写法就是保证这一点。另一个隐蔽场景训练集里某个状态从未被访问训练后的At也可能出现整行接近0过一段时间后预测出现NaN。解决方法是把Tolerance适当调大如1e-4而不是1e-8让EM早点停别在数值底部反复震荡。如果已经出现NaN把A0和B0重新初始化再训一次通常就能避开那个坏的局部最优。5.3 差分还原的索引错位与误差累积现象还原后的预测曲线比真实值整体滞后一个点或者预测误差随时间越来越大。原因滞后一个点通常是索引错位——diff比原始序列少一个点cumsum还原时基准点选错整条曲线就平移了一位。误差越来越大则是多步预测的必然结果每次预测误差都会累积到下一步的基准里。解决索引对齐用3.4节的方式te_aligned dy(trainLen1:end)因为diff(series)的第一个元素是series(2)-series(1)而pred_dy的第一个预测值是series(trainLen2)-series(trainLen1)的估计。误差累积的解决方法是缩短预测步长或者每预测一步就用真实观测更新状态分布——3.3节的代码就是这么设计的不要让模型在自回归模式下跑太久。另外建议同时看差分尺度的RMSE和还原尺度的MAPE。差分RMSE反映模型在“变化量”上的真实预测能力还原MAPE反映业务指标。两者背离很大时不是模型坏了是累积误差在起作用你要判断的是业务上到底关心哪个尺度。5.4 状态数过多导致状态语义丢失现象训练集BIC确实下降了但解码出的状态序列在相邻时间点上频繁跳变比如1-2-1-2-1-2业务上完全说不清每个状态代表什么。原因状态数N超出数据真实的模式数量时EM会拿多余的状态去拟合噪声。学出来的状态不连续转移矩阵的驻留概率很低。BIC在样本量不大时倾向于选择更简单的模型但如果你用AIC几乎一定选到最大的N。解决不要只看BIC同时看状态的平均驻留时间。驻留时间等于1 / (1 - A(i,i))如果某个状态的平均驻留时间小于3个时间步基本可以判定这个状态是噪声拟合出来的。我会在N3、N4、N5各跑一遍挑“状态序列连续且可解释”的那个而不是BIC最小的那个。一个实用的诊断手段把解码后的状态序列和原始时间序列画在同一张图上用不同颜色标记不同状态。如果某个状态对应的曲线段没有一致的波形特征这个状态就是多余的。这一步比任何数学指标都直观。5.5 MATLAB版本差异hmmtrain可能不在基础安装里现象在新版MATLAB里调用hmmtrain报错“Undefined function hmmtrain”或者提示需要Statistics and Machine Learning Toolbox。原因hmmtrain早期属于基础MATLAB后来被移到了Statistics and Machine Learning Toolbox。R2022a之后还改过内部实现参数名从MaxIter改成了MaxIterations老脚本直接跑挂。解决先运行ver(stats)确认工具箱是否安装再用which hmmtrain看它实际指向的路径。如果工具箱没装就需要自己实现一个离散HMM的EM更新逻辑不复杂E步做前向后向M步更新A和B的行归一化。资源包里的核心代码不依赖其他工具箱函数所以即使hmmtrain不可用你也能照3.1到3.4的流程自己写一个训练函数。还有一个容易忽略的点MATLAB的hmmtrain在较新版本启动时会输出一条deprecation warning提示这个函数将来可能被移除。这个warning不影响运行不用管。但如果你的公司有MATLAB版本升级计划提前把EM训练逻辑沉淀成自己的函数是值得的——不要把自己的核心流程绑在一个可能被移除的内置函数上。6. 验证与进阶让HMM从“跑通”到“可信”模型跑通、参数调完还不够。落地之前你要验证两件事状态有没有物理意义预测区间是否可靠。验证方法很简单——把解码出的状态序列和原始序列画在同一张图上不同颜色标记不同状态人工确认每个状态对应的波形特征是否一致。状态1是否总对应低谷状态2是否总对应尖峰。如果状态序列在相邻时间点上频繁跳变这个模型即使指标好看业务上也不能用。滚动多步预测推荐用递归法先把预测值离散化追加到观测序列再做下一步预测。关键代码是把3.3节的循环改成带步长h的形式h 7; % 预测未来7步 pred_multi zeros(h, 1); cur_obs_full tr_obs; for step 1:h [~, ~, fwd] hmmdecode(cur_obs_full, At, Bt); alpha fwd(:, end); obs_dist alpha * At * Bt; pred_multi(step) obs_dist * centers; [~, pred_obs] min(abs(centers - pred_multi(step))); cur_obs_full [cur_obs_full, pred_obs]; end注意递归法会把预测误差带进下一步的状态更新所以h不宜过大。我一般控制在原始序列一个完整模式周期的一半以内。HMM是概率模型不用它输出预测区间太浪费。obs_dist本身就是观测分布直接取分位数pred_cdf cumsum(obs_dist); lo_idx find(pred_cdf 0.05, 1); hi_idx find(pred_cdf 0.95, 1); pred_lo centers(lo_idx); pred_hi centers(hi_idx);离散桶的分位数只能落在桶边界精度有限。如果换连续高斯发射HMM可以用正态分布分位数得到平滑区间。这也是我建议熟手把离散版换成连续版的主要原因——不是为了提升点预测精度而是为了拿到平滑、真实的置信区间。在状态切换点附近HMM的区间会自然变宽这是它的隐藏优势。最后说一个习惯。整理这个资源包时我踩过最深的一个坑是差分还原的索引错位——整整半天时间训练集和测试集分开打的时候RMSE挺好一画图整条曲线整体左移一个点最后发现是diff和cumsum之间少了一个点。从那以后每次写时间序列代码我都会强制走一遍“原始数据长度→差分后长度→还原后长度”的对齐检查三行代码的事能省掉最不值钱的半夜。希望这份完整源码和示例数据能让你把HMM这个老模型用得明明白白。本文还有配套的精品资源点击获取 HMM时间序列预测实战：MATLAB实现隐马尔可夫模型预测
简介一套基于MATLAB的HMM隐马尔科夫模型时间序列预测完整实现方案面向时序预测学习者、信号处理及数据科学从业者解决如何用隐马尔科夫模型对一维观测序列进行训练、解码与未来值预测的问题。资源包含3个MATLAB脚本、3张结果示意图、1份说明文档和1个数据文件共8个文件压缩包大小为238KB轻量但环节完整。脚本覆盖参数训练、前向算法、Viterbi解码与Baum-Welch参数优化并提供MSE、MAE等评估指标的可视化对比内置数据文件包含训练与测试时序数据示意图直观展示模型结构与预测效果文档梳理完整实现步骤。目前已有3414人学习下载。通过运行这些代码可快速复现完整的时序预测流程掌握状态转移矩阵与观测发射矩阵的建模细节为后续改进模型或迁移到语音、金融等序列预测任务打下基础。1. 用HMM做时间序列预测先搞清楚它在预测什么先把结论放这儿HMM隐马尔科夫模型做时间序列预测预测的不是“下一时刻的数值”而是“下一时刻最可能处于哪个隐含状态以及在这个状态下观测值的分布”。这个区别决定了你能不能把它用对。这套MATLAB实现包含完整源码和示例数据覆盖数据预处理、HMM训练、隐状态解码、逐步预测和RMSE/MAE评估适合数据量不大但存在明显模式切换的时间序列——比如用电负荷的峰谷切换、风速的静风与大风交替、金融序列的震荡与趋势切换。新手可以照脚本跑通全流程熟手可以直接替换成自己的数据调整状态数和离散化参数。2. HMM的建模逻辑五个要素、三个问题、一次映射2.1 五元组状态、观测与三个概率矩阵HMM的数学结构是五元组状态集合S{s₁,...,s_N}观测集合O{o₁,...,o_M}连续观测则对应分布状态转移矩阵AN×NA(i,j)P(s_{t1}j|s_ti)观测概率矩阵BN×MB(j,k)P(o_tk|s_tj)连续观测时是高斯的均值与方差初始状态分布π1×N。映射到时间序列预测隐状态是“看不见的模式”。比如用电负荷数据里的“工作日早高峰状态”“午间平谷状态”“晚间高峰状态”观测是“实际测到的负荷值”。HMM的核心假设有两条t1时刻的状态只依赖t时刻的状态一阶马尔可夫t时刻的观测只依赖t时刻的状态条件独立。工程上这两个假设往往不严格成立但好处是模型的参数极少、训练极快、可解释性极强。你付出的代价是模型对长程依赖和复杂非线性关系的表达能力不如LSTM这类深度模型。2.2 三个基本问题在预测里的角色HMM有三个经典问题它们分别对应预测流程的不同环节。评估问题给定模型λ(A,B,π)和观测序列O计算P(O|λ)。这其实是模型选择的核心工具。你要决定状态数取3还是5不能靠拍脑袋而是分别训练模型、比较对数似然再用AIC/BIC做权衡。解码问题给定模型和观测求最可能的隐状态序列。这一步回答“当前处于哪个模式”Viterbi算法就是干这个的。学习问题给定观测序列估计模型参数A、B、πMATLAB里的hmmtrain就是Baum-Welch算法EM的实现。基本问题算法MATLAB函数在预测中的用途评估前向算法hmmdecode输出LOGPSEQ模型选择比较不同状态数下的对数似然解码Viterbihmmviterbi获得最可能的隐状态序列用于状态解释学习Baum-WelchEMhmmtrain估计转移矩阵A和观测矩阵B预测时三个问题按顺序上场先用学习问题在训练集上学出A和B再用解码或前向滤波得到当前状态分布最后用A和B推下一时刻观测的概率分布取期望作为点预测。这就是HMM做预测的完整链路每一步都有对应的MATLAB函数这也是这套源码能写得这么短的原因。2.3 离散观测还是连续观测MATLAB工具箱的边界MATLAB的hmmtrain、hmmviterbi、hmmdecode在设计之初面向的是离散观测符号比如DNA序列的ACGT、语音识别里的码本索引。也就是说B矩阵是一个N×M的离散概率表观测序列必须是1、2、3这样的正整数索引。这就带来一个直接后果原始的时间序列是连续数值喂给hmmtrain之前必须离散化。常见做法是分桶把训练集的观测值按等频或等距分成K个桶每个值映射成一个符号。代价是丢失桶内的数值精度点预测只能回到桶中心或桶均值预测曲线呈阶梯状。另一种做法是自己写高斯发射HMM的EM更新让每个状态的观测分布是N(μ_j, σ_j²)把B矩阵换成高斯密度。这样不需要离散化预测更平滑还能输出平滑的置信区间。但这需要自己实现前向、后向和EM的M步更新代码量大约多一百行。资源包以离散版本为主推方案因为hmmtrain足够稳定不容易在数值上翻车。如果你想用连续版本第6章会给一个改造思路。2.4 同门对比HMM vs ARIMA vs LSTM做时间序列预测的人通常先想到ARIMA或LSTM。ARIMA是线性模型擅长捕捉自相关结构但遇到状态突变比如负荷从平静状态突然跳入高峰状态时反应迟钝预测曲线会出现明显的滞后。LSTM在数据量足够大时表现很好但要调的参数多、可解释性差小样本下容易过拟合——这是“lstm时间序列预测python”这类问题下最常见的坑。HMM恰好站在两者中间它对几千个点的小样本友好EM迭代通常几十步就收敛它给出显式的状态语义业务上能解释“当前在哪个模式”它本质上是概率模型能输出预测区间而不只是点预测。代价是一阶马尔可夫假设偏弱对长程依赖捕捉有限。我给的建议是数据存在明显模式切换、样本量在几百到几千、需要向业务解释预测依据时用HMM你只有海量数据且只追求精度时LSTM更合适。资源包里的示例数据就是典型的模式切换序列用HMM的预测滞后通常比ARIMA小一个量级。3. 源码拆解从数据到预测值的完整链路3.1 数据加载、差分与观测离散化第一步是加载数据并做预处理。资源包里的hmm_data.mat存了一个变量series是单变量时间序列列向量带有明显的模式切换特征。换成你自己的数据时只要保证series是列向量即可。diff做一阶差分是很多时序预测的常规操作原始负荷值非平稳HMM学“变化量”的规律比学“绝对数值水平”更稳健。差分后归一化到[0,1]再用等频分桶把连续值映射成离散符号。%% 3.1 数据加载与观测离散化 clear; clc; rng(42); % 固定随机种子保证结果可复现 load(hmm_data.mat); % 变量名 series列向量 % 差分平稳化HMM学的是变化量的规律不是绝对水平 dy diff(series); y_norm (dy - min(dy)) / (max(dy) - min(dy)); % 映射到[0,1] % 等频分桶让每个桶的样本量均衡避免某些符号样本过少 K 10; % 桶数K10是经验值第4章讲怎么调 edges prctile(y_norm, linspace(0, 100, K 1)); % K1个边界 edges(end) edges(end) eps; % 防止最大值落到桶外 [~, obs] histc(y_norm, edges); obs(obs 0) 1; % histc在等于下边界时可能返回0强制归入1号桶这段代码做了四件事差分、归一化、分桶、符号映射。关键参数K直接决定B矩阵的规模N×K也决定预测精度。K太小预测退化成几个台阶K太大每个桶样本量不足、B矩阵出现大量近似零概率、训练过拟合。我一般取10到30具体看样本量。这里有个细节值得注意edges(end) edges(end) eps这行不是可有可无的。当某个观测值恰好等于最大值时histc会把它归到最后一个桶之外返回0导致obs出现0这个非法符号。加一个eps把上边界稍微抬高保证所有值都落在桶内。3.2 训练初始化、hmmtrain与Viterbi解码训练的核心是hmmtrain它的输入是观测符号序列、初始转移矩阵、初始观测矩阵输出训练后的A和B。初始值的给法直接影响是否收敛到合理解这里先给出能跑通的版本。%% 3.2 HMM训练与状态解码 trainLen round(0.8 * length(obs)); % 前80%训练后20%测试 tr_obs obs(1:trainLen); te_obs obs(trainLen1:end); nStates 3; % 隐状态数负荷场景下解释为 平/升/降 三个模式 A0 ones(nStates) / nStates; % 均匀初始化确保没有全零行 A0 A0 0.2 * eye(nStates); % 对角线加一点权重倾向状态驻留 A0 A0 ./ sum(A0, 2); % 行归一化 B0 ones(nStates, K) / K; % 观测矩阵均匀初始化 % 训练最大迭代200次容差1e-4 [At, Bt] hmmtrain(tr_obs, A0, B0, MaxIterations, 200, Tolerance, 1e-4); % 用Viterbi解码训练集的最可能状态序列 states hmmviterbi(tr_obs, At, Bt); % 评估训练集对数似然用于模型选择 [~, logp_train] hmmdecode(tr_obs, At, Bt); fprintf(训练完成log-likelihood %.4f\n, logp_train);逻辑说明A0 A0 0.2 * eye(nStates)这一步很多人会忽略但它很关键。如果初始A是全均匀矩阵EM迭代容易陷入对称解——三个状态学出来一模一样。对角线加权重相当于先验地认为“状态倾向于维持”这符合大多数物理时序的惯性特征。参数说明MaxIterations控制最大迭代次数数据量小时几十步就收敛200步是安全上限。Tolerance控制对数似然增幅阈值1e-4够用设成1e-8反而容易在数值底部来回震荡。注意旧版MATLAB2018a之前的参数名是MaxIter新版是MaxIterations如果报错先检查这里。hmmtrain不输出初始分布π它内部按均匀分布处理所以这里不需要初始化π。3.3 预测前向概率加权而不是只看最可能状态预测阶段最常见的错误是只用Viterbi解码出的单一状态去查转移矩阵。正确做法是用滤波分布加权因为Viterbi只给你一条最可能路径但当前状态本质上是一个分布尤其在状态边界附近两个状态的概率可能各占一半。%% 3.3 时间序列预测一步预测前向概率加权 centers (edges(1:K) edges(2:K1)) / 2; % 每个桶的中心值 pred_dy zeros(length(te_obs), 1); for i 1:length(te_obs) % 只使用到 t 时刻为止的真实观测避免未来信息泄露 cur_obs [tr_obs; te_obs(1:i-1)]; % i1时就是训练集 [~, ~, forward] hmmdecode(cur_obs, At, Bt); alpha_t forward(:, end); % 滤波分布1 x nStates next_obs_dist alpha_t * At * Bt; % 下一时刻观测分布1 x K pred_dy(i) next_obs_dist * centers; % 取期望作为点预测 end % 差分还原预测的是差分序列要累积回去才是真实值 pred_raw cumsum([series(trainLen 1); pred_dy]);逻辑说明hmmdecode的第三个输出forward是前向概率矩阵forward(:, end)就是截至当前时刻的滤波分布它比Viterbi的硬解码更稳健。next_obs_dist alpha_t * At * Bt这一步的本质是先把当前状态分布用转移矩阵传播到下一时刻再对每个可能状态求观测分布最后按状态概率加权得到的就是P(o_{t1}|o_1...o_t)。点预测取期望而不是众数这样预测曲线更平滑——众数只能落在桶中心期望是加权平均等于做了软输出。cumsum还原时以series(trainLen1)为基准这个基准点选择非常关键错了整条预测曲线就平移一位。提示hmmdecode在循环里每步都在重算整个前缀序列复杂度O(t·N²)。序列超过一万点时会明显变慢可以改成在线前向更新只保留上一时刻的前向向量每一步做一次迭代前向递推效果完全一样。3.4 评估RMSE、MAE、MAPE与可视化预测做完一定要评估不能只看拟合曲线。这里同时算了差分尺度的RMSE/MAE和还原尺度的MAPE因为差分还原会累积误差只看还原后的误差无法判断模型本身在“变化量”上的预测能力。%% 3.4 评估指标 te_aligned dy(trainLen1:end); % 与pred_dy对齐的真实差分值 rmse_diff sqrt(mean((te_aligned - pred_dy).^2)); mae_diff mean(abs(te_aligned - pred_dy)); mape_raw mean(abs((series(trainLen1:end) - pred_raw) ./ ... (series(trainLen1:end) eps))) * 100; fprintf(差分RMSE: %.4f\n, rmse_diff); fprintf(差分MAE : %.4f\n, mae_diff); fprintf(还原MAPE: %.2f%%\n, mape_raw); figure; plot(series(trainLen1:end), b-, LineWidth, 1.2); hold on; plot(pred_raw, r--, LineWidth, 1.2); legend(真实值, HMM预测值); xlabel(时间步); ylabel(序列值);索引对齐是这里最容易出错的地方dy是diff(series)的结果长度比series少一个点dy(trainLen1:end)恰好与te_obs一一对应。pred_raw的还原基准是series(trainLen1)所以真实值对应series(trainLen1:end)长度也一致。我在实际项目中踩过这个坑——训练时RMSE挺好一画图整条曲线整体左移一个点就是diff和cumsum之间的索引没对齐。4. 数据准备与参数调优四个旋钮四组经验值4.1 平稳化什么时候差分什么时候不差HMM本身不要求数据平稳因为它建模的是状态转移加观测发射的联合分布理论上非平稳序列也能学。但实际效果会很差如果原始序列有明显上升趋势模型会把“上升”也当做一个稳定模式新数据一旦出现不同趋势预测立刻失效。常见做法是先做一阶差分把序列从“数值水平”变成“变化量”HMM学的是变化量的模式。我判断是否差分的方法先画ACF图如果自相关系数衰减很慢就先差分差分后再看是否在零附近波动。差分一次不够就差分两次但二次差分后观测符号分布会非常集中——大部分值落在零附近——这时候要减少K否则很多桶是空的。资源包里默认做一阶差分。如果你自己的数据本身是平稳的比如白噪声附近的残差序列把差分那行注释掉即可但归一化和分桶要基于原始值重新计算边界不能沿用差分版本的edges。4.2 状态数N怎么选对数似然与AIC/BIC状态数N是HMM里最关键的旋钮。N太小模型无法区分不同模式预测退化成均值N太大过拟合训练集似然很高测试集一塌糊涂。标准做法是遍历N从2到8对每个N训练并记录对数似然用BIC权衡。%% 4.2 状态数选择BIC nList 2:8; bic zeros(size(nList)); for idx 1:length(nList) ns nList(idx); A0 ones(ns) / ns 0.2 * eye(ns); A0 A0 ./ sum(A0, 2); B0 ones(ns, K) / K; [At_t, Bt_t] hmmtrain(obs, A0, B0, MaxIterations, 200, Tolerance, 1e-5); [~, lp] hmmdecode(obs, At_t, Bt_t); k ns * (ns - 1) ns * (K - 1); % 转移矩阵观测矩阵的自由参数 bic(idx) -2 * lp k * log(length(obs)); end figure; plot(nList, bic, o-); xlabel(状态数N); ylabel(BIC); [~, bestN] min(bic); fprintf(BIC最优状态数: %d\n, nList(bestN));逻辑说明自由参数数量按转移矩阵去掉行归一化约束后的N(N-1)个、观测矩阵去掉每行归一化约束后的N(K-1)个来计算。初始分布π的参数在样本量面前影响极小BIC里可以忽略。用BIC而不是AIC是因为HMM参数多AIC的惩罚项太轻几乎总会选到N8。实际使用中BIC给出的N只是参考。N3对应“平/升/降”三个模式业务上说得通N4如果解出来的第四个状态说不清是什么那即使BIC选了4我也会退回3。模型是拿来解释业务的不是拿来拟合似然的。4.3 桶数K与预测平滑度的关系K的选择常常被忽略因为hmmtrain只要求符号是正整数很多教程随便取一个数。但K直接影响B矩阵的规模和预测曲线的平滑度。K太大每个桶的样本量少B矩阵的估计方差大K太小预测只能在少数几个桶中心之间跳。我常用的经验公式是K 10 round(sqrt(length(series) / 1000))。2000个点取K115000个点取K15。如果你用期望值做点预测3.3节的方式K偏小带来的阶梯效应会明显减轻因为期望是加权平均不是众数。如果你改成众数预测K偏小的话预测曲线会非常难看。K和N还有一个交互效应K越大观测矩阵B的自由参数越多BIC越倾向选择较小的N。所以实际调参的顺序是先定K再遍历N。反过来先定N再调K你会发现K怎么调BIC都不稳定。4.4 初始值与随机性让训练结果可复现hmmtrain的Baum-Welch是EM算法收敛到局部最优解初始值不同结果可能不同。资源包里的脚本统一用rng(42)固定随机种子并且A0、B0按确定性方式初始化这样每次跑结果一致。遇到换一台机器结果变化很大的情况先检查是不是没固定种子。初始转移矩阵对角线加的权重0.2对收敛结果影响很大。权重越大模型越倾向状态驻留学出来的状态平均停留时间越长权重太小模型容易学出高频抖动的状态序列业务上很难解释。0.2是我在负荷类时序上的经验金融数据我一般加到0.5因为金融状态切换频率更低。有一个快速检查初始值是否合理的方法训练完打印转移矩阵At看对角线元素是否显著大于非对角线。如果学到的主对角线都在0.8以上说明你的数据确实存在强驻留模式模型结构是合理的。如果对角线全部小于0.3大概率是状态数选多了模型在拿状态拟合噪声。5. 避坑指南HMM时间序列预测的五个典型翻车点5.1 观测序列编码必须是从1开始的连续正整数现象直接拿连续数值喂给hmmtrain报错“Error using hmmtrain ... 观测值必须是从1开始的整数”或者不报错但训练结果全是NaN。原因hmmtrain把观测符号当成离散索引要求max(obs) size(B,2)且min(obs) 1。如果你把归一化后的浮点数直接当obs索引自然是非法的。解决严格按3.1节的分桶流程走一遍分桶后检查min(obs)和max(obs)。注意histc在观测值恰好等于第一个边界时返回0所以obs(obs 0) 1这行必须保留。还有一个隐蔽场景如果某个观测值在edges之外比如测试集的数值比训练集更大分桶会得到K1甚至更大的非法符号。稳妥做法是用训练集的min/max先做一次clip再把clip后的值分桶。5.2 初始转移矩阵存在全零行训练直接崩现象hmmtrain迭代一两步后报错或者返回的At里有整行NaN。原因EM的E步要计算前向概率的归一化如果A0的某一行全零前向递推中经过该状态的所有路径概率恒为0分母为0数值上出现NaN。解决初始化A0时必须行归一化且每行至少有一个正概率。ones(nStates)/nStates 0.2*eye(nStates)这个写法就是保证这一点。另一个隐蔽场景训练集里某个状态从未被访问训练后的At也可能出现整行接近0过一段时间后预测出现NaN。解决方法是把Tolerance适当调大如1e-4而不是1e-8让EM早点停别在数值底部反复震荡。如果已经出现NaN把A0和B0重新初始化再训一次通常就能避开那个坏的局部最优。5.3 差分还原的索引错位与误差累积现象还原后的预测曲线比真实值整体滞后一个点或者预测误差随时间越来越大。原因滞后一个点通常是索引错位——diff比原始序列少一个点cumsum还原时基准点选错整条曲线就平移了一位。误差越来越大则是多步预测的必然结果每次预测误差都会累积到下一步的基准里。解决索引对齐用3.4节的方式te_aligned dy(trainLen1:end)因为diff(series)的第一个元素是series(2)-series(1)而pred_dy的第一个预测值是series(trainLen2)-series(trainLen1)的估计。误差累积的解决方法是缩短预测步长或者每预测一步就用真实观测更新状态分布——3.3节的代码就是这么设计的不要让模型在自回归模式下跑太久。另外建议同时看差分尺度的RMSE和还原尺度的MAPE。差分RMSE反映模型在“变化量”上的真实预测能力还原MAPE反映业务指标。两者背离很大时不是模型坏了是累积误差在起作用你要判断的是业务上到底关心哪个尺度。5.4 状态数过多导致状态语义丢失现象训练集BIC确实下降了但解码出的状态序列在相邻时间点上频繁跳变比如1-2-1-2-1-2业务上完全说不清每个状态代表什么。原因状态数N超出数据真实的模式数量时EM会拿多余的状态去拟合噪声。学出来的状态不连续转移矩阵的驻留概率很低。BIC在样本量不大时倾向于选择更简单的模型但如果你用AIC几乎一定选到最大的N。解决不要只看BIC同时看状态的平均驻留时间。驻留时间等于1 / (1 - A(i,i))如果某个状态的平均驻留时间小于3个时间步基本可以判定这个状态是噪声拟合出来的。我会在N3、N4、N5各跑一遍挑“状态序列连续且可解释”的那个而不是BIC最小的那个。一个实用的诊断手段把解码后的状态序列和原始时间序列画在同一张图上用不同颜色标记不同状态。如果某个状态对应的曲线段没有一致的波形特征这个状态就是多余的。这一步比任何数学指标都直观。5.5 MATLAB版本差异hmmtrain可能不在基础安装里现象在新版MATLAB里调用hmmtrain报错“Undefined function hmmtrain”或者提示需要Statistics and Machine Learning Toolbox。原因hmmtrain早期属于基础MATLAB后来被移到了Statistics and Machine Learning Toolbox。R2022a之后还改过内部实现参数名从MaxIter改成了MaxIterations老脚本直接跑挂。解决先运行ver(stats)确认工具箱是否安装再用which hmmtrain看它实际指向的路径。如果工具箱没装就需要自己实现一个离散HMM的EM更新逻辑不复杂E步做前向后向M步更新A和B的行归一化。资源包里的核心代码不依赖其他工具箱函数所以即使hmmtrain不可用你也能照3.1到3.4的流程自己写一个训练函数。还有一个容易忽略的点MATLAB的hmmtrain在较新版本启动时会输出一条deprecation warning提示这个函数将来可能被移除。这个warning不影响运行不用管。但如果你的公司有MATLAB版本升级计划提前把EM训练逻辑沉淀成自己的函数是值得的——不要把自己的核心流程绑在一个可能被移除的内置函数上。6. 验证与进阶让HMM从“跑通”到“可信”模型跑通、参数调完还不够。落地之前你要验证两件事状态有没有物理意义预测区间是否可靠。验证方法很简单——把解码出的状态序列和原始序列画在同一张图上不同颜色标记不同状态人工确认每个状态对应的波形特征是否一致。状态1是否总对应低谷状态2是否总对应尖峰。如果状态序列在相邻时间点上频繁跳变这个模型即使指标好看业务上也不能用。滚动多步预测推荐用递归法先把预测值离散化追加到观测序列再做下一步预测。关键代码是把3.3节的循环改成带步长h的形式h 7; % 预测未来7步 pred_multi zeros(h, 1); cur_obs_full tr_obs; for step 1:h [~, ~, fwd] hmmdecode(cur_obs_full, At, Bt); alpha fwd(:, end); obs_dist alpha * At * Bt; pred_multi(step) obs_dist * centers; [~, pred_obs] min(abs(centers - pred_multi(step))); cur_obs_full [cur_obs_full, pred_obs]; end注意递归法会把预测误差带进下一步的状态更新所以h不宜过大。我一般控制在原始序列一个完整模式周期的一半以内。HMM是概率模型不用它输出预测区间太浪费。obs_dist本身就是观测分布直接取分位数pred_cdf cumsum(obs_dist); lo_idx find(pred_cdf 0.05, 1); hi_idx find(pred_cdf 0.95, 1); pred_lo centers(lo_idx); pred_hi centers(hi_idx);离散桶的分位数只能落在桶边界精度有限。如果换连续高斯发射HMM可以用正态分布分位数得到平滑区间。这也是我建议熟手把离散版换成连续版的主要原因——不是为了提升点预测精度而是为了拿到平滑、真实的置信区间。在状态切换点附近HMM的区间会自然变宽这是它的隐藏优势。最后说一个习惯。整理这个资源包时我踩过最深的一个坑是差分还原的索引错位——整整半天时间训练集和测试集分开打的时候RMSE挺好一画图整条曲线整体左移一个点最后发现是diff和cumsum之间少了一个点。从那以后每次写时间序列代码我都会强制走一遍“原始数据长度→差分后长度→还原后长度”的对齐检查三行代码的事能省掉最不值钱的半夜。希望这份完整源码和示例数据能让你把HMM这个老模型用得明明白白。本文还有配套的精品资源点击获取