ARMAX有色噪声下的RML递推极大似然参数估计

简介：围绕极大似然法与递归极大似然（RML）的系统辨识资料包，面向控制理论、自动化、信号处理方向的学习者与工程人员，重点解决线性时不变系统传递函数或状态空间模型的参数估计与在线更新问题。压缩包共19个文件，约865KB，以MATLAB脚本（.m）承载RML算法实现与PRBS、M序列激励生成，.mat保存仿真与误差数据，.jpg给出参数估计误差和输出误差对比结果，另含.pptx理论讲解及Simulink工程文件，便于从理论、代码到实验验证贯通学习。已有189人学习下载。资料可帮助读者理解极大似然估计的基本思想、递归贝叶斯更新流程，以及非高斯噪声、有色噪声条件下的估计表现；通过误差曲线与真实参数对照，还能观察算法收敛性、模型适配性与在线调整效果，适合用作课程实验、算法复现和系统辨识入门到进阶的参考素材。

1. 残差一旦有颜色，最小二乘就交出了最优性

系统辨识里的极大似然法和 RML（Recursive Maximum Likelihood，递推极大似然）常被归进"高阶方法"那一栏：能不用就不用，用了也多半是原来的最小二乘怎么调都对不上频响。这个印象有道理，但它掩盖了一个更关键的事实——当你用 ARX 结构去拟合一批带着低频漂移或量化噪声的数据时，最小二乘估出来的参数其实是有偏的，而且偏得不是一点点；加数据、减步长，偏差不会消失，它只是稳稳地停在错误的位置上。极大似然要解决的不是"怎么把残差压得更小"，而是"怎么把噪声的形状一起写进模型"，让参数回到一致的位置。搞懂这套机制，你才会明白为什么同一个阶次、同一批数据，有人拟合出来的频响就是比你的稳，也才知道 RML 那几步递推到底在逼近什么。

2. 似然函数：把"数据出现的可能性"翻过来当成参数的函数

2.1 似然不是概率，它是一张针对 θ 的评分表

概率和似然的公式长得一模一样，但看的方向相反。概率是参数已知、问数据会以多大概率出现；似然是数据已经摆在桌上、把它当成参数 θ 的函数，问"哪一组 θ 最像是能生出这批数据的人"。这个视角切换在系统辨识里特别重要，因为观测序列 u、y 是死的，能动的只有 θ。

具体到一阶的情形，假设有 N 组独立观测，每组样本的密度是 p(y(t)|θ)，那么整批数据的联合密度就是各自密度的连乘。连乘有两个麻烦：一是数值上很快下溢，二是求导时乘积法则会炸成一堆项。取对数把它压成求和，单调性保证 argmax 不变，于是最大化似然等价于最大化对数似然。这一步不是数学装饰，它是后面所有递推算法能够成立的前提——只有变成求和，才能一项一项地在线累加。

打分表还有一个容易被忽略的性质：似然值的绝对大小没有意义，比值才有意义。你在两个候选参数之间做比较，看的是 L(θ₁)/L(θ₂)，常数因子可以整个扔掉。后面推导时会反复用到这一点，很多看似复杂的归一化项最后都会被这一刀砍掉。

2.2 高斯假设下，对数似然塌缩成残差平方和

把噪声设成独立同分布的高斯白噪声，方差为 λ，单点密度是(2πλ)^(-1/2)·exp(-ε²(t,θ)/(2λ))。连乘取对数，得到

log L = -(N/2)·log(2πλ) - (1/(2λ))·Σ ε²(t,θ)

第一项里没有 θ，只要 λ 固定，它就是个常数。于是最大化对数似然这件事，在数值上只剩下最小化Σ ε²(t,θ)——这就是最小二乘。很多人到这一步会有个错觉，以为"极大似然就是最小二乘换了个名字"。差别藏在 ε 的定义里：最小二乘里的残差是"y 减去回归量的线性组合"，而似然里的 ε 是预测误差，它可以带滤波器、可以递归、可以是输入输出的任意函数。名字一样，内涵完全不同。

顺着这个式子还能看出两个推广方向。如果噪声方差随时间变化，λ变成λ(t)，目标函数就成了Σ ε²(t)/λ(t)，这就是加权最小二乘，权重的物理含义是"这个时刻的观测有多可信"。如果噪声前后相关，高斯假设下协方差阵不再是λI，直接套最小二乘会丢信息，正确的做法是把相关性参数化成一个滤波器塞进模型——这正是下一章 ARMAX 要做的事。

2.3 噪声方差可以被解析消掉：集中似然的来历

λ 是个冗余参数。它对log L的影响有闭式解：对 λ 求偏导令其为零，得到λ̂ = (1/N)·Σ ε²(t,θ)。把这个结果代回对数似然，得到集中似然（profile likelihood）：

log L_c = -(N/2)·[log(2π) + 1 + log( (1/N)·Σ ε²(t,θ) )]

括号里前两项与 θ 无关，所以最大化集中似然等价于最小化log(Σ ε²)，又等价于最小化Σ ε²。形式上绕了一圈回到原点，但工程意义很大：优化时只需要在 θ 空间搜索，λ 每一步解析给出，少了一个维度，也避免了 λ 和 θ 耦合导致的病态。

这里有个前提常被忽略：能这样干净地消掉，是因为 λ 只以标量因子的形式出现在协方差里，θ 不参与方差的参数化。如果噪声模型的增益本身依赖 θ，代回之后目标函数的曲率会变，梯度的形式也要跟着改，直接照抄集中似然会得到错误的搜索方向。

用一段代码验证"高斯下的 ML 就是 LS"这件事，成本很低：

import numpy as np from scipy.optimize import minimize rng = np.random.default_rng(0) N, sigma = 2000, 0.5 x = rng.uniform(-1, 1, N) y = 2.5 * x + 1.0 + sigma * rng.standard_normal(N) # 真值 a=2.5, b=1.0 def neg_log_lik(theta): a, b = theta eps = y - (a * x + b) lam = np.mean(eps ** 2) # 解析消去的 λ̂ return 0.5 * N * np.log(lam) # 省略与 theta 无关的常数项 res = minimize(neg_log_lik, x0=[0.0, 0.0], method="Nelder-Mead") print("ML :", res.x) print("LS :", np.polyfit(x, y, 1)[::-1])

neg_log_lik里的lam就是集中似然消掉的那个方差，0.5*N*log(lam)对应-(N/2)log(λ̂)这一项。两组输出的差距应当落在数值误差量级，这说明在同方差高斯假设下两条路殊途同归；一旦把噪声换成有色，这个等式立刻破裂，破裂的部分就是极大似然真正值钱的地方。

3. 系统辨识里的极大似然：C(q) 才是那个隐藏主角

3.1 从 ARX 到 ARMAX，多出来的多项式改变了什么

ARX 结构写成A(q)y(t) = B(q)u(t) + e(t)，e 是白噪声。这个式子里噪声直接加在方程右边，等价于假设"残差已经白了"。如果真实噪声是白噪声，最小二乘就是极大似然，两者统一，没什么好争的。麻烦在于真实系统很少这么干净：传感器有低频漂移，采样保持和量化会引入相关性，机械谐振会让残差在某个频段鼓起来。这些成分让 e(t) 变成有色序列，而 ARX 强行把它当白噪声处理，等价于在回归量里漏掉了e(t-1)、e(t-2)…这些项，遗漏变量偏误就此产生。

ARMAX 加了第三项：A(q)y(t) = B(q)u(t) + C(q)e(t)，其中 A、B、C 都是首一多项式（C 的常数项为 1），e 仍是白噪声。C(q) 的角色是一个整形滤波器，它把白噪声染成你观测到的那个颜色。参数化之后，噪声不再是被动接受的误差，而是模型的一部分，和 A、B 一起被估计。估计出 C(q) 之后，1/C(q)就是一个白化滤波器，把它作用在残差上，残差才真正变白。

这里有个反直觉的结论值得记住：多估一组 C 的参数，并不会让模型"更复杂因而更容易过拟合"这么简单。它改变的是估计量的统计性质——在噪声确实有色的情况下，忽略 C 得到的最小二乘估计是有偏且不一致的，N 再大也救不回来；正确参数化 C 之后，极大似然估计是一致且渐近有效的。代价是优化问题从凸的二次型变成了非凸的，需要迭代，需要初值，也会收敛到局部解。这是一笔很典型的交易：用优化难度换统计正确性。

3.2 预测误差框架：最大化似然等价于最小化白化残差

把 ARMAX 两边同乘1/C(q)：

ε(t,θ) = [A(q)/C(q)]·y(t) - [B(q)/C(q)]·u(t)

这个 ε 就是模型的一步预测误差。写成更一般的形式，任何y = G(q,θ)u + H(q,θ)e的模型都能导出

ε(t,θ) = H⁻¹(q,θ)·[y(t) - G(q,θ)u(t)]

似然函数只依赖 ε 的平方和。这就把"极大似然"翻译成了"预测误差最小化"（PEM），两者在 gaussian 假设下完全等价。PEM 框架的好处是它不挑模型结构：只要你能写出预测器，就能构造 ε，就能套同一套优化和收敛性分析。ARMAX、OE、Box-Jenkins 在这一点上是统一的，区别只在 G 和 H 怎么参数化、共享不共享极点。

推导里有一个容易被跳过、但面试里经常被追问的细节：从 e 到 y 的变换带有雅可比项det(∂e/∂y)，似然里按理要出现log|det|。对 ARMAX 不用加，因为 A、C 都是首一多项式，y(t) 只依赖 e(t) 及更早的 e，这个变换是下三角的，对角元全为 1，行列式恰好是 1。换成 H 非首一的参数化，log det H就必须补回去，漏掉它会让估计系统性偏移——这是个安静的错误，残差看着正常，参数就是不对。

把这个构造写成代码，用scipy.signal.lfilter做两次滤波就够了：

import numpy as np from scipy.signal import lfilter def armax_pred_error(theta, u, y, na, nb, nc): """ theta = [a1..ana, b1..bnb, c1..cnc] 返回一步预测误差 eps(t) = (1/C)(A y - B u) """ a, b, c = theta[:na], theta[na:na+nb], theta[na+nb:] A = np.r_[1.0, a] # 首一 B = np.r_[0.0, b] # 补零对齐到 t-1 起 C = np.r_[1.0, c] v = lfilter(A, [1.0], y) - lfilter(B, [1.0], u) # A(q)y - B(q)u return lfilter([1.0], C, v) # 再被 C 白化 theta = np.array([-0.8, 0.5, -0.3]) # 任意一组参数，仅演示计算流程 eps = armax_pred_error(theta, u, y, na=1, nb=1, nc=1) V = np.mean(eps ** 2) # 这就是要最小化的目标

A、B、C三个系数向量直接决定滤波器；B前面补的那个 0 是因为 ARMAX 里输入从t-1开始起作用，不补零会让整条序列错位一个采样点，这种错位在闭环数据上极难从残差曲线上看出来。V就是集中似然里被最小化的那一项，后面所有迭代算法都是在它的地形上找路。

3.3 Fisher 信息阵：参数能估多准，跑算法之前就定了

梯度向量ψ(t,θ) = -∂ε(t,θ)/∂θ是后面所有算法的核心构件。在真值附近，Fisher 信息阵可以写成J = (1/λ)·E[ψ(t)ψᵀ(t)]，参数估计的协方差下界是J⁻¹，这就是 Cramér-Rao 下界。极大似然估计渐近达到这个下界，所以它被称为渐近有效；最小二乘在噪声有色时连一致性都没有，更谈不上有效性。

J 的表达式里没有 ε 本身，只有 ψ 的自相关。这个结构说明了一件很实际的事：参数能估多准，在跑算法之前就由输入信号的频谱和模型结构决定了。如果输入在某频段没有能量，对应的 ψ 分量接近零，J 就接近奇异，那块参数怎么估都是飘的。很多人调 RML 调不出来，问题不在算法，在激励信号——给一段阶跃或者一段直流偏置的数据，任何辨识算法都无能为力。J 还解释了为什么参数之间会强相关：当两个参数的 ψ 分量高度共线时，J 接近奇异，参数只能估出它们的某种组合，单独的值没有意义。

4. RML 的递推实现：每一步只做一次微小的牛顿迭代

4.1 把批处理牛顿法拆成在线更新

批处理极大似然的标准解法是牛顿法：θ_{k+1} = θ_k + [V''(θ_k)]⁻¹·(-V'(θ_k))。目标函数V = (1/2N)Σε²的一阶导是-(1/N)Σ ψ(t)ε(t)，二阶导的精确形式包含 ε 对 θ 的二阶导，计算量大且不稳定。工程上普遍采用高斯-牛顿近似，把二阶导里的 ε 项丢掉，只保留(1/N)Σ ψ(t)ψᵀ(t)。这个近似在残差接近白噪声时非常准，因为白噪声与 ψ 不相关，被丢掉的那一项期望为零。

把这个近似写成递推形式，就是 RML 的主更新：

θ̂(t) = θ̂(t-1) + K(t)·ε(t)

其中增益K(t) = P(t)ψ(t)，而 P 是近似 Hessian 逆的递推估计。P 的更新用矩阵求逆引理把"加一个新样本后重算逆矩阵"简化成秩一修正，这就是递归最小二乘那一套标准推导。参数估计每步只做一次向量加法和一次外积，计算量是 O(n²)，n 是参数个数，完全能塞进一个采样周期。

和 RELS（递推增广最小二乘）的区别就落在 ψ 的构造上。RELS 的回归向量直接用原始的 y、u 和历史残差；RML 的 ψ 要先过一遍1/C(q)滤波。这个差别看起来只是多了一层滤波器，实际上决定了估计量是否一致：RELS 忽略了 C 对回归量的影响，相当于用错误的权重做加权，噪声模型估得越准，这个权重误差造成的影响反而越明显。RML 把这一层补上了，代价是每步多两次 IIR 滤波。

4.2 梯度 ψ 的构造：三个 1/C 滤波器

对 ARMAX 逐项求偏导，ψ 的三个分块分别是：

负号来自ψ = -∂ε/∂θ这个定义，符号搞反是 RML 实现里最常见的翻车点，表现为参数朝反方向跑、几步之内发散。三个分量都要经过同一个1/C(q)，这是 RML 和 RELS 的分水岭：RELS 那里是原始数据直接进回归向量，RML 这里进的是滤波后的版本。滤波器1/C是 IIR 的，实现上要维护三组状态：输出的滤波历史、输入的滤波历史、残差的滤波历史。每步的递推是f(t) = x(t) - c₁f(t-1) - … - c_n f(t-n)，也就是用当前估计出的 C 系数对自身的滤波输出做反馈。

这里有一个先有鸡还是先有蛋的问题：算 ψ 需要 C，算 C 又需要 ψ。RML 的处理办法是随机逼近——用上一时刻（或刚更新完）的 C 估计去滤波，误差被吸收进收敛性证明里。只要参数缓慢变化、步长合适，这个近似不影响收敛到的点，只影响收敛路径。

4.3 增益矩阵与遗忘因子在做什么

K(t) = P(t)ψ(t)的量纲是"每单位预测误差应该把参数挪多远"。P 大意味着当前对这些参数没把握，步子迈大；P 小意味着已经收敛，步子收敛到微调。P 的初值p0通常取 10³ 到 10⁶ 量级，代表"先验极弱"，前几十步参数会大幅跳动，这是正常现象，不是算法坏了。

遗忘因子 λ 控制历史数据的权重衰减，P的更新里除以 λ 会让旧信息逐渐失效。λ=1 是无限记忆，参数收敛后不再变化，适合时不变系统；λ 取 0.95 到 0.999 是有限记忆，能跟踪慢时变，但代价是稳态方差变大——λ 越小，等效数据窗口越短，估计越抖。这里没有普适最优值，只能按系统的时变速率来选，经验上让等效窗口长度是系统时间常数的三到五倍比较稳。

分母里的λ + ψᵀPψ是个归一化项，防止 ψ 特别大时增益爆掉。它的存在让算法对输入幅值的绝对大小不敏感，但挡不住量纲失衡——如果 u 的单位是伏特、y 的单位是毫米，ψ 的两个分块能差好几个数量级，P 的条件数会烂掉，这时候要么归一化数据，要么给 P 加个对角正则。

4.4 一段可运行的核心循环

下面这段是 RML 的主循环，为了展示清晰，A、B、C 取相同的阶次 n0；工程里三者阶次可以不同，把切片下标改掉即可。

import numpy as np def rml_armax(u, y, n0, lam=0.995, p0=1e3): """ ARMAX 递推极大似然 y(t)+a1 y(t-1)+... = b1 u(t-1)+... + e(t)+c1 e(t-1)+... n0 : A/B/C 公共阶次 lam : 遗忘因子，1.0 为无限记忆 p0 : 协方差阵初值，越大表示先验越弱 """ N = len(y) npar = 3 * n0 th = np.zeros(npar) # [a, b, c] P = np.eye(npar) * p0 # 近似 Hessian 的逆 yh = np.zeros(n0) # y(t-1)..y(t-n0) uh = np.zeros(n0) # u(t-1)..u(t-n0) eh = np.zeros(n0) # eps(t-1)..eps(t-n0) fy = np.zeros(n0) # (1/C)y 的滤波历史 fu = np.zeros(n0) # (1/C)u 的滤波历史 fe = np.zeros(n0) # (1/C)eps 的滤波历史 th_log = np.zeros((N, npar)) for t in range(1, N): a, b, c = th[:n0], th[n0:2*n0], th[2*n0:] # (1) 当前参数下 <p> <a href="https://download.csdn.net/download/m0_38020234/11231566" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>