马尔可夫预测建模实战:从核心原理到MATLAB实现
2026/8/28 7:47:40 网站建设 项目流程

1. 从“无记忆性”到预测未来:马尔可夫预测的核心思想

在数学建模,尤其是涉及时间序列预测的赛题中,我们常常会遇到一类特殊的数据:它的下一个状态,只和当前状态有关,而与过去的历史路径无关。听起来有点反直觉,对吧?现实世界如此复杂,一个系统的未来怎么可能只取决于现在,而完全“忘记”过去呢?但恰恰是这种“无记忆性”的简化假设,催生了马尔可夫预测这一强大而实用的工具。我第一次在国赛的C题里用它分析一个排队系统的状态转移时,才真正体会到,好的模型不一定是完全真实的,但一定是抓住核心矛盾且可计算的。

马尔可夫预测,本质上是一种基于概率的预测方法。它不试图去拟合复杂的曲线,也不去挖掘深层的因果关系,而是聚焦于系统状态之间“跳转”的规律。比如,你明天的情绪是“开心”还是“沮丧”,可能很大程度上取决于你今天的心情,而和你上周是否中了彩票关系不大;一个网店明天的销量是“高”还是“低”,可能更依赖于今天的促销活动和流量,而不是上个月的销售数据。马尔可夫模型就是把这些“状态”和状态间“转移的可能性”给量化出来,然后像解方程一样,推算出未来某个时刻系统最可能处在什么状态,或者各个状态的概率分布是怎样的。

它的核心价值在于处理那些具有明显“状态”特征且转移具有一定随机性的系统。在数学建模竞赛中,从资源调度(如APMCM亚太赛的AGV路径问题)、市场占有率分析、到生态种群演变、甚至是一些社会行为预测,马尔可夫模型都有一席之地。它不像深度学习算法那样是个黑箱,其过程透明,结果有明确的概率解释,这对于需要清晰建模论文的竞赛来说,是一个巨大的优势。接下来,我们就抛开复杂的公式,从实际问题出发,看看怎么把马尔可夫预测这个工具用起来。

2. 马尔可夫链的“零件拆解”:状态、转移与概率矩阵

要搭建一个马尔可夫预测模型,我们得先搞清楚它的三个基本“零件”:状态、转移和概率矩阵。这就像拼乐高,零件认清了,组合起来就顺畅了。

2.1 如何定义“状态”:模型成败的第一步

定义状态是整个建模的起点,也是最考验对问题理解深度的一步。状态划分得太粗,会丢失信息,预测不准;划分得太细,会导致状态空间爆炸,计算复杂,且转移矩阵难以从有限数据中准确估计。

实战中的状态定义技巧:

  1. 基于业务逻辑离散化:对于连续变量(如销量、温度),不要直接使用原始值。应根据业务意义划分区间。例如,预测产品日销量,可以定义为:状态1: 低销量 (0-100件)状态2: 中销量 (101-500件)状态3: 高销量 (501件以上)。区间的划分可以基于历史数据的分布(如三分位数)、业务目标(如盈亏平衡点)或自然断点。
  2. 枚举分类状态:对于本身就是分类的数据,直接作为状态。比如机器运行状态:正常预警故障;天气状态:
  3. 状态组合(谨慎使用):当系统由多个维度决定时,可以考虑组合状态。例如,研究一个地区的经济-环境系统,状态可以是(经济好, 污染轻)(经济好, 污染重)(经济差, 污染轻)(经济差, 污染重)。但要注意,这会使状态数呈乘积增长,务必确保每个组合状态都有足够的历史数据支撑。

注意:状态必须是互斥完备的。即任意时刻,系统必须且只能处于其中一个状态。

2.2 计算状态转移概率矩阵:从历史数据中学习规律

转移概率矩阵P是马尔可夫模型的心脏。它是一个方阵,元素P_{ij}表示从状态i转移到状态j的概率。P的每一行之和必须等于1(因为从状态i出发,下一时刻必然转移到所有可能状态之一)。

如何从一串状态序列计算P假设我们有一串观测到的状态序列:[A, A, B, A, B, B, A, C, B, A](状态空间为{A, B, C})。

  1. 统计频数:首先统计从每个状态出发,转移到其他状态的次数。

    • 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
  2. 计算概率:将频数矩阵的每一行除以该行的总和。

    • 行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)

    得到转移概率矩阵P

    P = [ 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(t+1); 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)就等于一步转移概率矩阵Pk次方。 即: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%继续使用A,20%转用B,10%转用C。
  • 本月使用B品牌的顾客中,下月有10%转用A,80%继续用B,10%转用C。
  • 本月使用C品牌的顾客中,下月有5%转用A,5%转用B,90%继续用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*P=S, 等价于 P'*S'=S' % 增加一个约束条件:所有分量之和为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是从非吸收态到吸收态的转移矩阵。

对于吸收链,我们关心两个核心问题:

  1. 在吸收前,平均经过多少步?这需要计算基本矩阵N = (I - Q)^(-1),其元素N_{ij}表示从非吸收态i出发,在吸收前处于非吸收态j的平均次数。
  2. 最终被各个吸收态吸收的概率是多少?这可以通过计算B = N * R得到,其中B_{ij}表示从非吸收态i出发,最终被吸收态j吸收的概率。

这类问题在设备可靠性分析、贷款风险(坏账为吸收态)、竞赛淘汰赛制等场景中非常有用。在MATLAB中,核心计算就是矩阵求逆和乘法。

% 假设一个简单系统:状态1,2为非吸收态,状态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 模型检验:马尔可夫性的验证

我们一直假设过程满足马尔可夫性(无后效性)。如何检验?一个常用的方法是卡方检验

思路:比较“观测到的转移频数”与“假设马尔可夫性成立下期望的转移频数”是否有显著差异。

  1. 根据历史数据计算经验转移概率矩阵P_emp
  2. 假设过程是一阶马尔可夫链,那么从状态i经过两步转移到状态j的期望概率,应该等于(P_emp^2)_{ij}。对应的期望频数可以通过边际频数计算。
  3. 将观测到的两步转移频数与期望频数进行卡方检验。

如果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 马尔可夫预测的局限性

认识到局限性比会用模型更重要。

  1. 无后效性假设的脆弱性:很多真实系统的未来状态依赖于更长的历史。例如,股市价格、流行病传播,显然不是只依赖前一天的状态。这时需要用高阶马尔可夫链或其它模型(如隐马尔可夫模型、时间序列模型)。
  2. 状态空间定义的任意性:状态如何划分,直接影响模型。不同的划分方式可能得到截然不同的预测结果。这需要结合领域知识进行敏感性分析。
  3. 时齐性假设:转移概率不随时间变化,这在实际中很难满足。经济周期、季节性因素都会改变转移规律。可以考虑使用时变的转移矩阵,但参数估计会非常困难。
  4. 对数据量的要求:要准确估计转移矩阵,需要足够多的历史数据,特别是对于状态数较多的情况。如果某个状态i出现的次数很少,那么P(i, :)这一行的估计就会非常不可靠。
  5. 无法预测具体数值,只能预测分布:马尔可夫链预测的是处于各个状态的概率,而不是状态的具体取值(如具体的销量数字)。它提供的是概率意义上的趋势,而非精确值。

4.3 MATLAB实战中的高效与稳定技巧

  1. 处理大矩阵乘方:当预测步数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
  2. 稀疏矩阵优化:当状态数很多(比如成百上千),但转移矩阵非常稀疏(大多数转移概率为0)时,使用稀疏矩阵存储和运算可以极大节省内存和计算时间。

    P_sparse = sparse(P); % 将稠密矩阵转换为稀疏矩阵存储 % 后续的乘法运算会自动利用稀疏算法 P_sparse_k = P_sparse^k;
  3. 模拟(蒙特卡洛)方法:对于复杂的马尔可夫链(如非时齐、状态依赖),或者想得到预测结果的分布而不仅仅是期望,可以采用模拟方法。通过模拟成千上万条可能的未来路径,来统计预测结果的分布。

    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);

    这种方法非常灵活,可以处理解析方法难以解决的复杂情况,并且能给出预测的置信区间。

  4. 与其它模型的结合:马尔可夫链常作为更复杂模型的组成部分。例如,隐马尔可夫模型假设状态是不可观测的,我们只能看到由状态生成的观测值。这在语音识别、金融序列分析中应用极广。在MATLAB中,有Statistics and Machine Learning Toolbox提供的hmmestimate,hmmdecode,hmmviterbi等函数可以用于HMM的训练和解码。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询