简介：面向系统辨识与参数估计学习者，《极大似然法.zip》围绕极大似然法与递归极大似然（RML）展开，适用于希望理解LTI系统传递函数或状态空间模型参数估计的读者，也可作为控制、信号处理课程的实验补充。压缩包共19个文件，约865KB，包含6个MATLAB脚本（m）用于算法实现与仿真、2个mat数据文件保存实验数据、9张jpg结果图展示参数估计误差与输出对比，另有1份pptx理论讲义和1个tmw工程文件，便于从理论、代码到结果对照学习。已有189人学习下载。资料重点放在RML的递归更新流程、白噪声与有色噪声下的估计表现，以及PRBS激励和误差分析，读者可据此复现实验、观察参数收敛过程，并对照图表检查算法实现是否正确。整体体量轻便，适合具备概率统计与系统辨识基础的学习者按模块查阅。
1. 从一组带噪伺服数据说起：极大似然法.zip 解决的是有偏估计问题
最小二乘跑出拟合度 0.92，残差一看还有明显自相关，这种结果在系统辨识里最坑人。极大似然法.zip 是一份围绕极大似然法做系统辨识的资源包，核心落在 RML 递推极大似然上，用来处理 ARMAX 模型在有色噪声下的参数估计，而不是只给一页公式。它适合做课程设计、算法复现、工程预研的人：手里有输入输出数据，想在线更新参数，又发现普通最小二乘估计有偏。反过来，数据本身激励不足、系统强非线性、噪声脉冲很重，先别指望 RML 救场。拆包顺序我一般这样走：先找入口脚本跑通示例数据，再对照第 3 章换自己的数据，最后才动阶次和遗忘因子。
2. 极大似然与 RML 的底层账本：目标函数、噪声模型和递推边界
2.1 为什么系统辨识里要用极大似然而不是只做最小二乘
系统辨识里最常用的 ARMAX 结构可以写成 A(q)y(t)=B(q)u(t-d)+C(q)e(t)，其中 e(t) 是白噪声，A、B、C 都是滞后算子多项式。最小二乘在 ARX 结构、噪声为白的情况下和极大似然是等价的，但一旦噪声有色，回归向量里的 y(t-k) 和残差相关，最小二乘估计就不再一致。极大似然法的思路不是直接拟合输出，而是把噪声模型 C(q) 也当成参数的一部分，用 C(q) 去白化残差，让预报误差 ε(t)=A(q)/C(q)y(t)-B(q)/C(q)u(t-d) 尽量接近白噪声。负对数似然在 σ² 集中掉之后，等价于最小化预报误差平方和，这也是 RML 递推能落到预测误差框架里的原因。
RML 的“递推”两个字很关键。批处理极大似然要把所有数据攒起来做非线性优化，数据一长、模型一复杂，计算量和内存都上去了；RML 每来一个采样点就更新一次参数和协方差，适合在线辨识、慢时变系统跟踪，也适合嵌入式或工控现场那种不想存整段数据的场景。代价是它对初值、遗忘因子和数据尺度更敏感，参数跑飞的情况并不少见。常见做法是先用最小二乘或 ARX 给一个初值，再切到 RML 细调，这样收敛过程会稳很多。
从选型上看，RML 适合单输入单输出、离散线性、噪声可以用有理谱描述的场合。它比 OE 模型多了一个噪声多项式 C(q)，所以对噪声结构有一定建模能力；它比 BJ 模型简单，计算量小，工程上更容易落地。边界也很清楚：强非线性、死区、摩擦、大延迟跳变，这些不是 ARMAX 能兜住的；数据里如果全是阶跃或者常值段，持续激励不够，RML 会把噪声当信息学进去，参数看着收敛，实际预测一塌糊涂。系统辨识里判断一个模型能不能用，不能只看拟合度，残差白噪声检验和交叉验证至少要做一项。
把这条链路拆成可执行步骤，我一般按下面顺序走：
- 明确模型结构，先定 na、nb、nc 和延迟 d，不要一上来就加阶次。
- 对输入输出去均值，必要时做尺度归一化，避免协方差矩阵数量级差异过大。
- 用 ARX 或最小二乘给一组初值，把 RML 的起点放在合理区域。
- 构造回归向量和预报误差梯度，进入递推循环。
- 每步记录参数、残差和协方差迹，收敛后做残差自相关与交叉验证。
提示：RML 对数据尺度很敏感，输入输出量纲差两个数量级以上时，先归一化再递推，否则协方差矩阵容易病态。
2.2 RML 的递推形式和包内实现边界
RML 的递推形式和递归最小二乘很像，差别在梯度向量 ψ(t) 不是简单的回归向量，而是预报误差对参数的偏导。参数更新写成 θ(t)=θ(t-1)+K(t)ε(t)，增益 K(t)=P(t-1)ψ(t)/(λ+ψᵀ(t)P(t-1)ψ(t))，协方差更新 P(t)=(P(t-1)-K(t)ψᵀ(t)P(t-1))/λ。λ 是遗忘因子，λ=1 时记忆全部历史，λ 小于 1 时更看重近期数据。ψ(t) 的计算依赖当前参数下的 A(q)、B(q)、C(q) 滤波，所以每步都要先算预报误差，再算梯度，再更新参数和协方差。这个顺序不能乱，否则递推出来的参数会漂。
包内实现边界要提前看清：它大概率是离散时间、线性时不变、SISO 的 ARMAX 递推，输入是 u、y 两列向量，输出是参数估计和残差序列。常见做法是把 na、nb、nc 作为入口参数，延迟 d 单独给，遗忘因子 λ 和协方差初值 P0 放在配置区。它不会自动帮你判断数据有没有持续激励，也不会自动选阶次，这些得自己补。若包内是 MATLAB 实现，矩阵运算和滤波函数会比较直接；若拿到的是 Python 版，递推公式和参数含义不变，只是把矩阵运算换成 numpy，滤波用 scipy.signal.lfilter 也能写。
| 符号 | 含义 | 调参时的实际影响 |
|---|---|---|
| na | 输出自回归阶次 | 太小残差白不了，太大容易过拟合 |
| nb | 输入多项式阶次 | 和延迟 d 一起决定输入通道记忆 |
| nc | 噪声多项式阶次 | 不够则残差自相关，过大则参数方差变大 |
| d | 输入延迟 | 错一拍，B 参数会整体偏移 |
| λ | 遗忘因子 | 越小跟踪快，但稳态波动大 |
| P0 | 协方差初值 | 太大前期猛冲，太小收敛慢 |
RML 的收敛判断不能只看参数变化量。我一般同时看三个量：参数增量的范数、预报误差的滑动均方、协方差迹。参数增量小但残差还有强自相关，说明模型结构不对，不是没收敛。协方差迹掉到很小但拟合没改善，可能是初值把算法带进了局部平坦区。这个时候别硬调 λ，先回去检查 nc 和延迟，或者换一段激励更充分的数据。系统辨识的“玄学”多半来自数据，不是算法本身。
3. 把 RML 跑起来：回归向量、递推循环与收敛落盘
3.1 数据分帧与回归向量构造
拿到包后先不要改核心函数，先看入口脚本怎么读数据、怎么分帧。我一般把数据切成辨识段和验证段，辨识段至少覆盖几个完整动态过程，验证段留一段没参与递推的数据。下面这段 MATLAB 代码构造 ARMAX 的回归向量和预报误差梯度所需的滞后项，重点是去均值、延迟对齐和矩阵预分配。参数 na、nb、nc 从配置区传入，d 是输入延迟，N 是样本数。
function [phi, ylag, ulag, elag] = build_regressors(u, y, e, na, nb, nc, d) % 输入去均值，避免直流分量混进参数 u = u(:) - mean(u(:)); y = y(:) - mean(y(:)); N = length(y); % 预分配回归矩阵，列顺序：y 滞后、u 滞后、e 滞后 phi = zeros(N, na + nb + nc); ylag = zeros(N, na); ulag = zeros(N, nb); elag = zeros(N, nc); for t = max([na, nb + d, nc]) + 1 : N % 输出滞后项 y(t-1)...y(t-na) ylag(t, :) = y(t-1:-1:t-na)'; % 输入滞后项 u(t-d)...u(t-d-nb+1) ulag(t, :) = u(t-d:-1:t-d-nb+1)'; % 噪声滞后项 e(t-1)...e(t-nc) elag(t, :) = e(t-1:-1:t-nc)'; phi(t, :) = [ylag(t, :), ulag(t, :), elag(t, :)]; end end这段代码的逻辑是把 ARMAX 的回归项按时间对齐，y 和 u 的滞后起点不同，d 错一位后面 B 参数就会整体偏移。参数说明：na 控制输出自回归项数量，nb 控制输入项数量，nc 控制噪声项数量，d 是纯延迟采样点数。预分配矩阵是为了避免递推过程中反复扩容，数据长的时候这个细节很影响耗时。去均值放在函数入口，避免把直流工作点当成动态参数辨识出来。
注意：e 的滞后项在第一次递推时还没有历史残差，前 max(na,nb+d,nc) 个点要跳过，或者用零初始化。不要拿验证段的残差回灌到辨识段，否则交叉验证会虚高。
3.2 RML 递推核心循环
核心循环里每步要做四件事：算预报误差、算梯度、更新参数、更新协方差。下面这段函数把 RML 的递推骨架写出来，θ 的顺序和 phi 的列顺序一致，P 是协方差矩阵，lambda 是遗忘因子。代码里保留了残差记录和参数轨迹，方便后面画收敛曲线。
function [theta_hist, e_hist, P_hist] = rml_core(u, y, na, nb, nc, d, lambda, P0, theta0) N = length(y); npar = na + nb + nc; theta = theta0(:); P = P0 * eye(npar); e = zeros(N, 1); theta_hist = zeros(N, npar); P_hist = zeros(N, 1); for t = max([na, nb + d, nc]) + 1 : N % 当前参数下的回归向量，e 用上一步残差 [phi_t, ~, ~, ~] = build_regressors(u(1:t), y(1:t), e(1:t), na, nb, nc, d); phi = phi_t(t, :)'; % 一步预报误差 e(t) = y(t) - phi' * theta; % 近似梯度用回归向量，严格 RML 需按 C(q) 滤波 psi = phi; % 增益与协方差更新 K = P * psi / (lambda + psi' * P * psi); theta = theta + K * e(t); P = (P - K * psi' * P) / lambda; % 记录轨迹 theta_hist(t, :) = theta'; e_hist(t) = e(t); P_hist(t) = trace(P); end end这里的逻辑说明要分两层看：phi 是当前时刻的回归向量，e(t) 是用上一步参数算出来的一步预报误差，K 是增益，P 是参数协方差。参数说明：lambda 越接近 1 越稳，越小于 1 跟踪越快但稳态噪声越大；P0 是协方差初值，数据归一化后一般取 1 到 100 之间，太大前期参数会猛冲；theta0 是初值，建议先用 ARX 或最小二乘给。代码里 psi 直接用了 phi，这是简化写法，严格 RML 要把 phi 按 C(q) 滤波后再当梯度用，包内实现如果做了这一步，以包内函数为准。
提示：递推前先把 u、y 归一化到相近量级，P0 不要拍脑袋给 1e6。协方差爆炸多数不是算法错，是尺度和初值一起翻车。
3.3 收敛判断与结果落盘
收敛判断我一般看滑动窗口内的参数增量和残差均方，而不是单步值。下面这段脚本在递推结束后计算最后 200 步的参数变化范数、残差均方和协方差迹，并保存结果。参数 window 是滑动窗口长度，tol_theta 是参数增量阈值，tol_p 是协方差迹阈值。
window = 200; dtheta = vecnorm(diff(theta_hist(end-window:end, :), 1, 1), 2, 2); resid_mse = mean(e_hist(end-window:end).^2); p_trace = P_hist(end); fprintf('dtheta=%.3e, resid_mse=%.3e, traceP=%.3e\n', ... mean(dtheta), resid_mse, p_trace); save('rml_result.mat', 'theta_hist', 'e_hist', 'P_hist');这段脚本的输出不是用来“证明收敛”的，而是用来判断有没有异常。dtheta 很小、resid_mse 也小、traceP 平稳，才算比较健康；dtheta 很小但 resid_mse 还很大，多半是模型结构不对；traceP 突然掉到零，可能是 P 更新除了零或者数据长期无激励。保存 theta_hist 和 e_hist 是为了后面画参数轨迹和残差自相关，不要只存最后一行参数，否则出了问题连回看的机会都没有。
4. 参数怎么设：na/nb/nc、遗忘因子与初值的三组对照
4.1 模型阶次：从 AIC 和残差自相关两头夹
阶次选择是 RML 最容易拍脑袋的地方。na 太小，输出自回归不够，残差会留下明显自相关；na 太大，参数方差上去，验证段预测反而变差。nb 要和延迟 d 一起看，d 错一拍，nb 再大也补不回来。nc 是噪声模型阶次，ARMAX 里它负责白化残差，nc 不够时残差自相关图会在低阶滞后上翘，nc 过大则参数估计变慢，数据短的时候尤其明显。我一般先跑一组小范围组合，记录 AIC、验证段拟合度和残差自相关最大值，再挑一个不是“拟合最高”而是“残差最白、验证不塌”的组合。
| na | nb | nc | AIC | 验证段 fit (%) | 残差自相关最大值 |
|---|---|---|---|---|---|
| 2 | 2 | 1 | -1.82e3 | 86.4 | 0.21 |
| 2 | 2 | 2 | -1.91e3 | 88.1 | 0.09 |
| 3 | 2 | 2 | -1.90e3 | 87.6 | 0.12 |
| 3 | 3 | 3 | -1.88e3 | 85.9 | 0.15 |
这张表是示意对照，不是固定答案。看表时不要只盯 AIC 最小，验证段 fit 掉太多说明过拟合；残差自相关最大值超过 0.2 时，先怀疑 nc 和延迟，而不是继续加 na。系统辨识里“阶次越高越好”是典型的坑，尤其数据只有几千点时，高阶模型会把噪声也拟合进去，RML 递推看着收敛，换一段数据就原形毕露。
4.2 遗忘因子与协方差初值：RML 的“记忆长度”
遗忘因子 λ 决定 RML 记住多久的历史。λ=1 时等价于无限记忆，适合时不变系统；λ=0.98 到 0.995 适合慢时变，但稳态参数会有一点点波动；λ 太小，比如 0.9，参数会跟着噪声跑，看起来跟踪很快，实际预测发散。协方差初值 P0 决定前期步长，P0 大前期更新猛，P0 小收敛慢。数据归一化后，我一般从 P0=1 或 P0=10 起步，先跑 λ=1 看能不能收敛，再根据时变程度往下调。
| λ | 适用场景 | 风险 |
|---|---|---|
| 1.0 | 时不变、数据量中等 | 时变时跟踪慢 |
| 0.995 | 慢时变、工控采样 | 稳态略有波动 |
| 0.98 | 较快时变 | 噪声被当信息，参数抖 |
| 0.95 以下 | 短时跟踪 | 稳态方差大，慎用 |
lambda = 0.995; % 慢时变系统先取 0.995 P0 = 10; % 归一化后协方差初值 theta0 = theta_ls; % 先用最小二乘给初值 [theta_hist, e_hist, P_hist] = rml_core(u_id, y_id, na, nb, nc, d, lambda, P0, theta0);这段设置的意思是：先用一个不太激进的遗忘因子把递推跑通，确认残差和参数轨迹正常，再调 λ。P0 不要和 λ 同时大改，否则出了问题分不清是初值还是记忆长度导致的。参数说明：lambda 越接近 1 越稳，P0 越大前期步长越大，theta0 是起点。初值给得好，RML 前几百步就不会乱冲，后面残差检验才有意义。
4.3 初值策略：两步法先 LS 后 RML
RML 是局部优化，初值差的时候很容易收敛到奇怪的地方。我一般用两步法：先跑 ARX 最小二乘拿 θ_LS，把 A、B 部分填进 RML 的初值，C 部分先用 1 或者小系数；然后固定 λ=1 跑一遍 RML，看残差有没有白化；最后再开遗忘因子做时变跟踪。下面这段代码演示怎么从 ARX 结果拼初值，phi_ls 是 ARX 的回归矩阵，y_id 是辨识段输出。
theta_ls = phi_ls \ y_id; % ARX 最小二乘初值 theta0 = zeros(na + nb + nc, 1); theta0(1:na+nb) = theta_ls(1:na+nb); % A、B 部分直接搬 theta0(na+nb+1:end) = 0.1; % C 部分给小系数，避免初始残差过大这段代码的逻辑是只把 A、B 的可靠部分搬过来，C 不给大初值，因为噪声多项式初值太大时预报误差会剧烈跳变。参数说明：theta_ls 是 ARX 估计结果，theta0 是 RML 起点，C 部分给小系数是工程上的保守做法。若包内已经有初始化函数，优先用包内的，但要知道它大概在做什么，否则换一组数据初值失效时无从排查。
5. 避坑与排查：RML 跑飞、残差不白、换数据就变的四个现场
5.1 现象：递推几步后参数爆炸
现象很直观：theta_hist 前几十步还正常，后面某个参数突然到 1e6 量级，残差也跟着炸。原因通常不是 RML 公式写错，而是协方差初值过大、数据没归一化、回归矩阵病态三者叠加。输入输出量纲差太多时，psiᵀPpsi 数量级失衡，增益 K 变得很大，一步就把参数推飞。解决顺序是先把 u、y 去均值并缩放到相近方差，再把 P0 降到 1 到 10，最后检查回归矩阵条件数。若条件数仍然很大，减少 na 或 nc，或者换一段激励更均匀的数据。血泪经验是：先怀疑尺度和初值，再怀疑算法。
5.2 现象：参数收敛但残差自相关明显
参数轨迹看起来很稳，dtheta 很小，但残差自相关图在滞后 1 到 3 上还有 0.3 以上。原因多半是 nc 不够、延迟 d 错位、或者系统里存在 ARMAX 描述不了的动态。先加 nc 看残差自相关是否下降，再扫 d 从 0 到 3 看验证段 fit 哪个最高，最后检查输入输出是不是有饱和或死区。若加 nc 后参数方差明显变大但残差仍不白，别继续加阶次，换 OE 或 BJ 结构试试。系统辨识里残差白不了就是模型没抓住动态，硬调 RML 参数只是把问题往后拖。
5.3 现象：换一段数据结果完全变
辨识段拟合很好，验证段一跑就崩，参数符号甚至反了。原因通常是激励不足、遗忘因子过小、或者阶次过高过拟合。检查输入段有没有持续变化，阶跃太少时先做预滤波或换数据；λ 如果小于 0.98，先回到 1 跑时不变版本；阶次组合如果 AIC 很小但验证 fit 很差，说明模型在背噪声。解决方法是留出独立验证段，用一步预测误差和交叉验证说话，不要用辨识段残差均方给自己发合格证。
5.4 现象：RML 比最小二乘还差
有人跑完发现 RML 的验证 fit 还不如 ARX，就觉得极大似然法没用。常见原因是噪声模型初值太差、C(q) 不稳定、或者数据里噪声根本不是有理谱。先把 C 部分初值调小，固定 λ=1 跑；再检查 C(q) 的根有没有跑到单位圆外；最后对比 ARX 和 ARMAX 的残差自相关，如果 ARX 残差本来就接近白噪声，RML 的优势不会明显。RML 不是万能升级包，噪声结构不对时它确实可能更差。判断标准是残差白噪声和验证误差，不是方法名字。
6. 进阶验证：交叉验证与频域残差给 RML 结果上保险
6.1 交叉验证与一步预测误差
递推结束后别急着收工，拿验证段做一步预测，算 fit 和残差均方。下面这段代码用最终参数在验证段上滚动预测，y_hat 是一步预报输出，fit 按 100*(1-NRMSE) 算。参数说明：theta_final 是递推最后一步参数，u_val、y_val 是验证段，na、nb、nc、d 与辨识段保持一致。
y_hat = zeros(size(y_val)); for t = max([na, nb+d, nc]) + 1 : length(y_val) [phi_t, ~, ~, ~] = build_regressors(u_val(1:t), y_val(1:t), e_val(1:t), na, nb, nc, d); y_hat(t) = phi_t(t, :) * theta_final; end fit = 100 * (1 - norm(y_val - y_hat) / norm(y_val - mean(y_val)));这段代码的逻辑是用历史真实输出和残差构造回归向量，只预测下一步，不做多步滚动，避免误差累积把问题放大。fit 只是入口指标，还要看残差自相关和频谱。若验证段 fit 比辨识段低超过 10 个百分点，优先怀疑过拟合或数据段动态不一致。
6.2 频域残差与置信区间
时域残差看着小，不代表频域干净。我一般再看残差的自相关和功率谱，确认低频段没有残留峰。下面这张表是我常用的验证组合，三项都过再认为模型可用。
| 验证项 | 通过标准 | 不通过时先查 |
|---|---|---|
| 验证段 fit | 与辨识段差距小于 10% | 阶次、激励、数据段一致性 |
| 残差自相关 | 95% 置信带内 | nc、延迟 d |
| 残差频谱 | 无明显低频峰 | 噪声模型、预滤波 |
从那以后我每次跑 RML 都强制先做数据分帧和残差白噪声检验，参数轨迹再漂亮也不跳过这两步。希望帮到你。
本文还有配套的精品资源，点击获取
