1. 从“无记忆性”到预测未来马尔可夫预测的核心思想在数学建模尤其是涉及时间序列预测的赛题中我们常常会遇到一类特殊的数据它的下一个状态只和当前状态有关而与过去的历史路径无关。听起来有点反直觉对吧现实世界如此复杂一个系统的未来怎么可能只取决于现在而完全“忘记”过去呢但恰恰是这种“无记忆性”的简化假设催生了马尔可夫预测这一强大而实用的工具。我第一次在国赛的C题里用它分析一个排队系统的状态转移时才真正体会到好的模型不一定是完全真实的但一定是抓住核心矛盾且可计算的。马尔可夫预测本质上是一种基于概率的预测方法。它不试图去拟合复杂的曲线也不去挖掘深层的因果关系而是聚焦于系统状态之间“跳转”的规律。比如你明天的情绪是“开心”还是“沮丧”可能很大程度上取决于你今天的心情而和你上周是否中了彩票关系不大一个网店明天的销量是“高”还是“低”可能更依赖于今天的促销活动和流量而不是上个月的销售数据。马尔可夫模型就是把这些“状态”和状态间“转移的可能性”给量化出来然后像解方程一样推算出未来某个时刻系统最可能处在什么状态或者各个状态的概率分布是怎样的。它的核心价值在于处理那些具有明显“状态”特征且转移具有一定随机性的系统。在数学建模竞赛中从资源调度如APMCM亚太赛的AGV路径问题、市场占有率分析、到生态种群演变、甚至是一些社会行为预测马尔可夫模型都有一席之地。它不像深度学习算法那样是个黑箱其过程透明结果有明确的概率解释这对于需要清晰建模论文的竞赛来说是一个巨大的优势。接下来我们就抛开复杂的公式从实际问题出发看看怎么把马尔可夫预测这个工具用起来。2. 马尔可夫链的“零件拆解”状态、转移与概率矩阵要搭建一个马尔可夫预测模型我们得先搞清楚它的三个基本“零件”状态、转移和概率矩阵。这就像拼乐高零件认清了组合起来就顺畅了。2.1 如何定义“状态”模型成败的第一步定义状态是整个建模的起点也是最考验对问题理解深度的一步。状态划分得太粗会丢失信息预测不准划分得太细会导致状态空间爆炸计算复杂且转移矩阵难以从有限数据中准确估计。实战中的状态定义技巧基于业务逻辑离散化对于连续变量如销量、温度不要直接使用原始值。应根据业务意义划分区间。例如预测产品日销量可以定义为状态1: 低销量 (0-100件)状态2: 中销量 (101-500件)状态3: 高销量 (501件以上)。区间的划分可以基于历史数据的分布如三分位数、业务目标如盈亏平衡点或自然断点。枚举分类状态对于本身就是分类的数据直接作为状态。比如机器运行状态正常、预警、故障天气状态晴、阴、雨、雪。状态组合谨慎使用当系统由多个维度决定时可以考虑组合状态。例如研究一个地区的经济-环境系统状态可以是(经济好 污染轻)、(经济好 污染重)、(经济差 污染轻)、(经济差 污染重)。但要注意这会使状态数呈乘积增长务必确保每个组合状态都有足够的历史数据支撑。注意状态必须是互斥且完备的。即任意时刻系统必须且只能处于其中一个状态。2.2 计算状态转移概率矩阵从历史数据中学习规律转移概率矩阵P是马尔可夫模型的心脏。它是一个方阵元素P_{ij}表示从状态i转移到状态j的概率。P的每一行之和必须等于1因为从状态i出发下一时刻必然转移到所有可能状态之一。如何从一串状态序列计算P假设我们有一串观测到的状态序列[A, A, B, A, B, B, A, C, B, A]状态空间为{A, B, C}。统计频数首先统计从每个状态出发转移到其他状态的次数。从A出发序列中A出现了4次作为起点。看每个A后面是什么A-A: 第1个A后是A 第4个A后是B 第7个A后是C。所以A-A出现1次A-B出现1次A-C出现1次。最后一个A是序列末尾无后续不计入。因此从A出发的转移总次数为3次。从B出发B出现了3次作为起点。B-A(第3个B后是A)B-B(第5个B后是B)B-?(第9个B后是A但序列结束这里序列是...B, A]所以第9个B后是A应计入)。所以B-A出现2次B-B出现1次。总次数3。从C出发C只出现1次C-B出现1次。总次数1。我们可以列出频数矩阵F:到A 到B 到C 从A 1 1 1 从B 2 1 0 从C 0 1 0计算概率将频数矩阵的每一行除以该行的总和。行A: (1,1,1) / 3 (0.333, 0.333, 0.333)行B: (2,1,0) / 3 (0.667, 0.333, 0.000)行C: (0,1,0) / 1 (0.000, 1.000, 0.000)得到转移概率矩阵PP [ 0.333 0.333 0.333 0.667 0.333 0.000 0.000 1.000 0.000 ]在MATLAB中的实现对于更长的序列手动计算不现实。我们可以用MATLAB向量化操作快速计算。假设states是一个包含状态索引如1,2,3的向量。% 假设 states 是状态索引序列例如 [1, 1, 2, 1, 2, 2, 1, 3, 2, 1] num_states max(states); % 状态总数 P zeros(num_states); for t 1:length(states)-1 from states(t); to states(t1); P(from, to) P(from, to) 1; end % 将频数转换为概率 row_sums sum(P, 2); % 避免除以0对和为0的行即某些状态在历史中从未作为起点出现进行处理 row_sums(row_sums 0) 1; P P ./ row_sums; % 利用广播机制每行除以对应的和 disp(转移概率矩阵 P:); disp(P);2.3 一步与多步转移预测的时空延伸我们得到的P是一步转移概率矩阵它描述了经过一个时间单位如一天、一月后的状态变化。那如果要预测两步、三步甚至更远呢这里就用到了切普曼-柯尔莫哥洛夫方程。简单来说k步转移概率矩阵P(k)就等于一步转移概率矩阵P的k次方。 即P(k) P^k为什么可以直观理解要从状态i经过两步到状态j第一步必须先跳到某个中间状态r然后从r跳到j。对所有可能的中间状态r求和正好就是矩阵乘法的定义。因此P(2) P * P P^2 依此类推。在MATLAB中计算多步预测% 已知初始状态分布 S0 例如 S0 [0.2, 0.5, 0.3] 表示初始时刻处于状态1,2,3的概率 % 预测 k 步后的状态分布 Sk k 5; % 预测5步后 P_k P^k; % 计算k步转移矩阵 Sk S0 * P_k; % 初始分布右乘k步转移矩阵 disp([预测 , num2str(k), 步后的状态概率分布:]); disp(Sk);一个关键假设时齐性。我们上面计算的前提是转移概率矩阵P不随时间改变。这在短期预测或系统相对稳定时近似成立。如果系统存在季节性或趋势则需要使用非时齐马尔可夫链或其它模型这大大增加了复杂性。在数学建模中通常先假设时齐性然后在模型检验部分讨论其局限性。3. 预测实战以市场占有率分析为例理论讲起来总是抽象的我们用一个经典的数学建模案例——市场占有率预测——来走一遍完整的流程。假设市场上有A、B、C三个品牌每月初消费者可能因为广告、口碑等原因更换品牌。我们通过市场调查得到了上个月的消费者流动情况。3.1 问题构建与数据准备已知数据转移频数本月使用A品牌的顾客中下月仍有70%继续使用A20%转用B10%转用C。本月使用B品牌的顾客中下月有10%转用A80%继续用B10%转用C。本月使用C品牌的顾客中下月有5%转用A5%转用B90%继续用C。当前市场占有率初始状态分布A: 30%, B: 45%, C: 25%。问题1预测下个月、三个月后各品牌的市场占有率。问题2长期来看市场会趋于一个稳定的分布吗这个稳定分布是什么3.2 模型建立与计算过程首先根据数据写出一步转移概率矩阵P。注意行表示“从”哪个状态列表示“到”哪个状态。我们按A, B, C的顺序。P [ 0.70 0.20 0.10 0.10 0.80 0.10 0.05 0.05 0.90 ]初始状态向量S0 [0.30, 0.45, 0.25]。在MATLAB中求解% 定义转移矩阵和初始状态 P [0.70, 0.20, 0.10; 0.10, 0.80, 0.10; 0.05, 0.05, 0.90]; S0 [0.30, 0.45, 0.25]; % 预测下个月一步的市场占有率 S1 S0 * P; disp(下个月市场占有率预测:); fprintf(A: %.2f%%, B: %.2f%%, C: %.2f%%\n, S1*100); % 预测三个月后三步的市场占有率 P_3 P^3; % 计算三步转移矩阵 S3 S0 * P_3; disp(三个月后市场占有率预测:); fprintf(A: %.2f%%, B: %.2f%%, C: %.2f%%\n, S3*100); % 尝试计算长期稳定分布极限分布 % 方法求解方程 S * P S 且 S各分量之和为1 % 即 S * (P - I) 0 其中I是单位阵 % 转化为求解 (P - I) 的零空间并归一化 I eye(3); A (P - I); % 转置是因为我们要解 S*PS 等价于 P*SS % 增加一个约束条件所有分量之和为1即 sum(S) 1 A [A; ones(1,3)]; b [zeros(3,1); 1]; % 前三个方程是齐次的最后一个方程和为1 % 使用左除求解最小二乘解因为方程可能超定或欠定这是稳健的做法 S_stable (A \ b); disp(长期稳定市场占有率极限分布:); fprintf(A: %.4f, B: %.4f, C: %.4f\n, S_stable);运行这段代码你会得到类似以下结果下个月预测A ≈ 26.5% B ≈ 43.3% C ≈ 30.3%。三个月后预测A ≈ 22.2% B ≈ 40.1% C ≈ 37.7%。长期稳定分布A ≈ 18.2% B ≈ 36.4% C ≈ 45.5%。结果分析从预测可以看出品牌C虽然初始占有率最低但由于其客户忠诚度最高90%留存率并且能从A、B品牌吸引少量客户其市场份额在未来会持续增长。品牌A的客户流失相对严重份额下降最快。长期来看市场将稳定在A:18.2%, B:36.4%, C:45.5%的格局。这个稳定分布是转移矩阵P的内在属性与初始分布S0无关只要P满足一定正则条件。这意味着无论市场初期格局如何在当前的客户流动规律下最终都会收敛到这个比例。3.3 模型扩展带吸收态的马尔可夫链在某些问题中存在一些“吸收态”一旦进入就无法离开比如机器“故障”状态、游戏“结束”状态。这类链被称为吸收马尔可夫链。它的转移矩阵可以写成标准形式P [ Q R 0 I ]其中I是单位矩阵对应吸收态Q是非吸收态之间的转移矩阵R是从非吸收态到吸收态的转移矩阵。对于吸收链我们关心两个核心问题在吸收前平均经过多少步这需要计算基本矩阵N (I - Q)^(-1)其元素N_{ij}表示从非吸收态i出发在吸收前处于非吸收态j的平均次数。最终被各个吸收态吸收的概率是多少这可以通过计算B N * R得到其中B_{ij}表示从非吸收态i出发最终被吸收态j吸收的概率。这类问题在设备可靠性分析、贷款风险坏账为吸收态、竞赛淘汰赛制等场景中非常有用。在MATLAB中核心计算就是矩阵求逆和乘法。% 假设一个简单系统状态12为非吸收态状态3为吸收态 % P [0.5, 0.3, 0.2; % 从状态1出发 % 0.2, 0.6, 0.2; % 从状态2出发 % 0.0, 0.0, 1.0]; % 从状态3吸收态出发 P [0.5, 0.3, 0.2; 0.2, 0.6, 0.2; 0.0, 0.0, 1.0]; Q P(1:2, 1:2); % 非吸收态部分 R P(1:2, 3); % 到吸收态的部分 I eye(size(Q)); N inv(I - Q); % 基本矩阵 % 计算从每个非吸收态出发在吸收前经历的平均步数包括自身 % 即基本矩阵N的每一行之和 avg_steps_before_absorption sum(N, 2); disp(从各非吸收态出发被吸收前的平均步数:); disp(avg_steps_before_absorption); % 计算最终被吸收态吸收的概率 B N * R; disp(从各非吸收态出发最终被吸收态吸收的概率:); disp(B);4. 马尔可夫模型的检验、局限与MATLAB实战技巧建立一个马尔可夫模型后我们不能直接拿着结果就去写论文。必须对模型进行检验并清醒地认识其局限性。同时在MATLAB实现中也有一些效率与稳定性的技巧。4.1 模型检验马尔可夫性的验证我们一直假设过程满足马尔可夫性无后效性。如何检验一个常用的方法是卡方检验。思路比较“观测到的转移频数”与“假设马尔可夫性成立下期望的转移频数”是否有显著差异。根据历史数据计算经验转移概率矩阵P_emp。假设过程是一阶马尔可夫链那么从状态i经过两步转移到状态j的期望概率应该等于(P_emp^2)_{ij}。对应的期望频数可以通过边际频数计算。将观测到的两步转移频数与期望频数进行卡方检验。如果p值大于显著性水平如0.05则没有足够证据拒绝原假设即可以认为数据满足马尔可夫性。MATLAB实现简化版思路由于完整的卡方检验涉及构建复杂的列联表在竞赛中我们可以采用一种更直观的“可视化”或“近似”检验计算自相关系数对于一个时间序列如果它是一阶马尔可夫的那么它的二阶及以上的自相关系数应该很小理论上对于马尔可夫过程自相关系数呈指数衰减。我们可以计算序列滞后1、2、3阶的自相关系数观察衰减情况。% 将状态序列转换为数值序列例如 states [1,2,1,3,2,...] [acf, lags] autocorr(states, NumLags, 5); figure; stem(lags, acf); xlabel(滞后阶数); ylabel(自相关系数); title(状态序列自相关图); grid on; % 观察acf(2), acf(3)...是否迅速接近0分状态检验对于每个状态i收集所有从i出发的后续状态。然后检验这些后续状态是否与再前一个状态独立。这可以通过比较条件分布来实现但操作较复杂。在时间有限的数学建模竞赛中更常见的做法是在模型假设部分明确提出马尔可夫性假设并在模型评价与灵敏度分析部分讨论该假设若不成立对结果可能产生的影响。这是一种务实的策略。4.2 马尔可夫预测的局限性认识到局限性比会用模型更重要。无后效性假设的脆弱性很多真实系统的未来状态依赖于更长的历史。例如股市价格、流行病传播显然不是只依赖前一天的状态。这时需要用高阶马尔可夫链或其它模型如隐马尔可夫模型、时间序列模型。状态空间定义的任意性状态如何划分直接影响模型。不同的划分方式可能得到截然不同的预测结果。这需要结合领域知识进行敏感性分析。时齐性假设转移概率不随时间变化这在实际中很难满足。经济周期、季节性因素都会改变转移规律。可以考虑使用时变的转移矩阵但参数估计会非常困难。对数据量的要求要准确估计转移矩阵需要足够多的历史数据特别是对于状态数较多的情况。如果某个状态i出现的次数很少那么P(i, :)这一行的估计就会非常不可靠。无法预测具体数值只能预测分布马尔可夫链预测的是处于各个状态的概率而不是状态的具体取值如具体的销量数字。它提供的是概率意义上的趋势而非精确值。4.3 MATLAB实战中的高效与稳定技巧处理大矩阵乘方当预测步数k很大时直接计算P^k可能导致数值不稳定元素趋于0或1。对于求极限分布更稳健的方法是求解特征向量。因为稳态分布π满足πP π即π是P的转置矩阵对应于特征值1的左特征向量。[V, D] eig(P); % 计算P转置的特征值和特征向量 % 找到特征值接近1的特征向量 idx find(abs(diag(D) - 1) 1e-10); if ~isempty(idx) pi_vec V(:, idx(1)); % 取对应的左特征向量行向量 pi_vec abs(pi_vec); % 取绝对值概率非负 pi_vec pi_vec / sum(pi_vec); % 归一化 disp(通过特征向量法求得的稳态分布:); disp(pi_vec); end稀疏矩阵优化当状态数很多比如成百上千但转移矩阵非常稀疏大多数转移概率为0时使用稀疏矩阵存储和运算可以极大节省内存和计算时间。P_sparse sparse(P); % 将稠密矩阵转换为稀疏矩阵存储 % 后续的乘法运算会自动利用稀疏算法 P_sparse_k P_sparse^k;模拟蒙特卡洛方法对于复杂的马尔可夫链如非时齐、状态依赖或者想得到预测结果的分布而不仅仅是期望可以采用模拟方法。通过模拟成千上万条可能的未来路径来统计预测结果的分布。num_simulations 10000; num_steps 20; current_state 1; % 假设当前处于状态1 future_states zeros(num_simulations, num_steps); for sim 1:num_simulations state current_state; for step 1:num_steps % 根据当前状态state的转移概率分布P(state, :)随机选择下一个状态 next_state randsample(1:size(P,2), 1, true, P(state, :)); future_states(sim, step) next_state; state next_state; end end % 分析第10步后状态的分布 step_to_analyze 10; dist histcounts(future_states(:, step_to_analyze), 1:size(P,2)1) / num_simulations; disp([通过模拟得到的第, num2str(step_to_analyze), 步状态分布:]); disp(dist);这种方法非常灵活可以处理解析方法难以解决的复杂情况并且能给出预测的置信区间。与其它模型的结合马尔可夫链常作为更复杂模型的组成部分。例如隐马尔可夫模型假设状态是不可观测的我们只能看到由状态生成的观测值。这在语音识别、金融序列分析中应用极广。在MATLAB中有Statistics and Machine Learning Toolbox提供的hmmestimate,hmmdecode,hmmviterbi等函数可以用于HMM的训练和解码。