系统模型阶选择问题

林夕小杰

832人浏览 · 2025-10-14 16:26:01

林夕小杰 · 2025-10-14 16:26:01 发布

就是决定模型有多少自由度——太少拟合不好（欠拟合），太多就抓住噪声（过拟合）。

6.3.1概念

常见使用到选阶的模型如AR(p)，MA(q)，ARMA(p,q)等模型的阶p，q在实际中，是一个待确定的参数，怎样确定最优的p值，是模型阶的确定问题。最优阶的选择也可以通过确定一定的目标准则得到最优解,Akaike较早研究了这个问题,并导出几个准则

1、最终预测误差准则(Final PredictionError,FPE),常用于 AR 模型选择，思路接近 AIC。

2、Akaike信息准则(Akaike Information Criterion,AIC),AIC=−2logL+2k，L 是模型的似然，k是参数数。倾向于稍复杂一点的模型（惩罚弱）。

3、最小描述长度准则(Minimizes The Description Length,MDL),信息论角度，类似 BIC。使准则函数最小的阶,确定为模型阶。

4、传输函数准则(Criterion Autoregressive Transfer Function,CAT)，最小化准则函数的阶选为模型的阶。

遗憾的是,这几个模型阶的准则并不能总给出相同的模型阶估计,在缺乏先验信息的情况下,可以实验用不同准则确定的阶,并对最终的谱估计结果进行分析和选择。

6.3.2 最终预测误差准则(Final PredictionError,FPE)

1、背景及介绍

FPE 是由日本统计学家 赤池弘次（Hirotugu Akaike） 在 1969 年提出的。
他当时在研究时间序列的自回归模型（AR 模型），想找到一个可以定量衡量模型预测能力的指标。

FPE 是 一种模型阶数选择方法。当你用一个数学模型（比如自回归 AR 模型）去描述信号时，会遇到一个问题：我要用多少阶？也就是，用多少个过去的样本来预测现在的样本最合适？

阶太低 → 模型太笨，抓不住规律；阶太高 → 模型太聪明，连噪声也学进去。

FPE 准则的目的就是：
在“拟合得好”和“不过度复杂”之间找到平衡点。估计“这个模型在未来预测新数据时的平均误差有多大”，选最小的那个阶。

2、基本形式（以AR模型为例）

考虑一个 p 阶自回归模型：

其中：

xt 是信号；
et 是误差项（白噪声）；
N是样本长度；
是残差的均方误差（模型拟合后留下的噪声能量）。

FPE 公式是：

1. 第一部分：

它代表模型在训练数据上的拟合误差。
越小说明模型拟合得越好。

但是！模型复杂（高阶）时，残差方差总会更小，因为模型参数多，可以“强行”贴合数据。
这就是过拟合的危险。

2. 第二部分：

这是惩罚项。
当 p 变大（模型更复杂）时，这个比值会迅速增大，对复杂模型进行惩罚。

总的来说：

若模型太简单 → 拟合误差大，FPE 大；
若模型太复杂 → 惩罚项变大，FPE 也大；
只有在最合适的阶数时，FPE 最小。

3、怎么选阶

1.设定一个可能的阶数范围，比如 p = 1 到 10

2.对每个 p：

拟合 AR(p) 模型；
计算残差方差；
套入公式；

找到 FPE 最小的 p，就是最佳阶数。

4、matlab实现

1、思路流程

生成一个已知阶数的时间序列（例如 AR(3) 模型）；
让 MATLAB 自动去尝试不同的阶；
计算每个模型的 FPE；
找到 FPE 最小的位置，就是“最佳阶数”

2、具体代码

%% -------------------------------
%  最终预测误差准则 (FPE) 自动选阶实战
%  作者：木夕水
%  -------------------------------

clc; clear; close all;

% Step 1: 生成一个已知阶的 AR(3) 信号
N = 2500;                             % 样本长度
a_true = [1, -0.75, 0.5, -0.25];     % AR(3) 系数, 注意a1,a2,a3
noise = randn(N,1)*0.5;              % 白噪声
x = filter(1, a_true, noise);        % 生成信号

% Step 2: 设置候选阶数范围
maxOrder = 10;                       % 最大测试阶数
FPE_values = zeros(maxOrder,1);      % 存放每个阶的 FPE

% Step 3: 逐阶拟合并计算 FPE
for p = 1:maxOrder
    model = ar(x, p, 'ls');          % 使用最小二乘法拟合 AR(p)
    FPE_values(p) = fpe(model);      % MATLAB自带函数 fpe(model)
end

% Step 4: 绘制结果
figure;
plot(1:maxOrder, FPE_values, 'bo-','LineWidth',1.5);
xlabel('模型阶数 p');
ylabel('最终预测误差 (FPE)');
title('FPE 准则选阶结果');
grid on;

% Step 5: 找出FPE最小点
[~, bestOrder] = min(FPE_values);
hold on;
plot(bestOrder, FPE_values(bestOrder), 'r*', 'MarkerSize',10, 'LineWidth',2);
text(bestOrder+0.3, FPE_values(bestOrder), sprintf('  最佳阶数 p = %d', bestOrder));

fprintf('基于 FPE 准则，最佳阶数为 p = %d\n', bestOrder);

3、代码讲解

Step 1：生成信号

filter(1, a_true, noise) 创建一个 AR(3) 信号。
这里的真实模型是：

（注意 MATLAB 的滤波器系数顺序与公式符号的关系。）

Step 2：定义候选阶数范围

我们打算从 1 阶到 10 阶试一遍，看哪个最合适。

Step 3：拟合模型并计算 FPE

ar(x,p,'ls') 使用最小二乘法估计 AR 模型参数。
fpe(model) 是 MATLAB 内置函数，会根据 Akaike 公式计算：

这一步完全自动完成。

Step 4：绘制结果

FPE 随着阶数的变化一般呈“先下降后上升”的趋势，最低点就是最佳阶。

Step 5：输出结果

打印并标记出最佳阶数，通常应该接近我们生成信号时的真实阶（这里是 3）。

4、运行结果

运行后，你会看到一张折线图：

横轴：模型阶数 p
纵轴：对应的 FPE 值

FPE 曲线一般会在低阶时下降（模型越复杂拟合越好），
到某个点后上升（过拟合导致泛化能力下降）(数据量较少)。
曲线最低点对应的 p，就是最佳阶。

5、优缺点

优点	缺点
思想直观，计算简单	假设模型是线性、高斯、平稳的
能很好地处理中等样本量	当样本太小（N≈p）时会不稳定（分母接近 0）
与预测性能直接相关	对非线性或非平稳序列效果有限

6.3.3 Akaike信息准则(Akaike Information Criterion,AIC)

1、背景及介绍

AIC与FPE都是通过计算“未来预测误差的估计”来选阶，只是FPE 从“预测误差”的角度看，AIC 从“信息量”角度看。FPE 关注“误差的大小”，AIC 关注“信息的损失”，结果指向一样，既不太傻也别太精。

AIC不是凭空蹦出来的“神秘公式”，而是对“模型在未来样本上的平均预测误差”的一种数学近似。Akaike 当年（1973 年）从统计推断的角度，把“预测”问题变成了一个关于信息量的度量。

Akaike 认为，一个模型的优劣在于它和真实系统之间的信息距离。如果你有一个真实的数据生成过程 f(x)，和一个你自己假设的模型 g(x∣θ)，那模型越接近真相越好。

AIC 本质上是 FPE 的统计推广版。只是 FPE 针对线性高斯模型推导，而 AIC 是更一般的形式（用信息论表述）。

2、基本形式

假设模型是线性高斯的，例如

AIC准则为：

L：模型的似然（likelihood），就是模型对数据的拟合好坏。拟合越好，L 越大。

−2ln(L)：把“拟合好坏”换成一个“误差形式”的分数，越小越好。

k：模型中要估计的参数个数（模型阶数相关）。

2k：惩罚项——模型越复杂，惩罚越重。

AIC = 拟合误差 + 模型复杂度惩罚

AIC常见公式形式

其中

N：样本点个数。

：当模型阶数为 p 时的残差方差。

p：模型的阶数，也就是参数数目。

3、如何选阶

我们试探模型该多复杂，可以从 1 阶、2 阶、3 阶……一直到10 阶，每试一个阶，就计算它的 AIC。然后选出 AIC 最小 的那个阶数。那一阶的模型，就被认为是“既不过拟合又不欠拟合”的最优模型。

4、matlab实现

1、思路流程

生成一个已知阶数的时间序列（例如 AR(3) 模型）；
让 MATLAB 自动去尝试不同的阶；
计算每个模型的 AIC；
找到 AIC 最小的位置，就是“最佳阶数”

2、完整代码

%% -------------------------------
%  AIC准则 ( AIC) 自动选阶实战
%  作者：木夕水
%  -------------------------------
clear; clc; close all;

% 1. 生成一个明显的 AR(3) 信号
N = 1000;                      % 样本点数足够长
rng(1);                        % 固定随机数种子
a_true = [1, -0.75, 0.5, -0.25]; % 真实AR(3)系数
e = randn(N,1);                % 高斯白噪声
y = filter(1, a_true, e);      % 生成信号

% 2. 设置要尝试的最大阶数
maxOrder = 15;
aic_values = zeros(maxOrder,1);

% 3. 用更鲁棒的Burg方法计算每个阶的AIC
for p = 1:maxOrder
    model = ar(y, p, 'burg');           % Burg方法估计
    noise_var = model.NoiseVariance;    % 残差方差
    N = length(y);
    aic_values(p) = N * log(noise_var) + 2 * p;  % AIC公式
end

% 4. 找出AIC最小对应的最优阶
[~, bestOrder] = min(aic_values);

% 5. 绘制AIC曲线
figure;
plot(1:maxOrder, aic_values, 'o-', 'LineWidth', 1.5, 'MarkerSize', 7);
hold on;
plot(bestOrder, aic_values(bestOrder), 'rp', 'MarkerFaceColor', 'r', 'MarkerSize', 12);
text(bestOrder+0.5, aic_values(bestOrder), ...
    ['  最优阶 = ', num2str(bestOrder)], 'FontSize', 12, 'Color', 'r');

xlabel('模型阶数 p');
ylabel('AIC 值');
title('基于 AIC 的模型自动选阶');
grid on;

3、代码讲解

第一步，生成一个AR（3）信号

第二步，设置最大尝试阶数15

第三步，用burg方法计算每个阶的AIC

第四步，找出最小AIC值，并绘制曲线表示

4、运行结果

运行后，你会看到一张折线图：

横轴：模型阶数 p
纵轴：对应的AIC值

AIC 曲线通常呈“U”型。左边模型太简单，右边太复杂，谷底位置就是最佳阶数。

5、优缺点

优点	缺点
通用性强	容易过拟合
计算简单	依赖模型假设
预测性能良好	对小样本表现不佳

6.3.4 最小描述长度准则(Minimizes The Description Length,MDL)

1、背景及介绍

MDL 其实是AIC，更严格更保守的版本，MDL更严谨、更偏向真实模型选择，而不是单纯预测性能。哲学思想是：用一个模型把数据压缩得更短，那么这个模型更接近真实。让模型 + 残差的总描述长度最小。

一个好的模型应该：

不要太复杂（否则你要花很多字节描述模型本身）；
要能解释数据（否则还得额外花字节描述误差部分）。

2、基本形式

对于AIC公式：

他有两个部分：

第一个部分：模型拟合数据的好坏。越小越好；

第二个部分2p：模型复杂度惩罚。参数越多，惩罚越大。

MDL和AIC很想，只是惩罚更严格，MDL公式为：

它和 AIC 几乎一模一样，只是把惩罚项从固定的 2p 改成了随样本量增长的。

这意味着：（1）当 N 很小时，AIC 和 MDL 差不多。（2）当 N 很大时，MDL 的惩罚会更强，更倾向于选更简单的模型。

其中：

N：样本长度
：该阶模型的残差方差
p：模型阶数

3、如何选阶

选出最优阶数：

4、matlab实现

1、思路流程

生成一个已知阶数的时间序列（例如 AR(3) 模型）；
让 MATLAB 自动去尝试不同的阶；
计算每个模型的 MDL；
找到 MDL 最小的位置，就是“最佳阶数”

2、完整代码

%% ----------------------------------------
%  基于 MDL（最小描述长度准则）的自动选阶实战
%  作者：木夕水
% ----------------------------------------
clear; clc; close all;

%% 1. 生成一个已知的 AR(3) 信号
N = 1000;                       % 信号长度（样本点数）
rng(1);                         % 固定随机数种子，保证结果可重复
a_true = [1, -0.75, 0.5, -0.25]; % 真实AR(3)系统参数
e = randn(N,1);                 % 生成白噪声（高斯随机数）
y = filter(1, a_true, e);       % 根据AR系数生成信号

%% 2. 设置搜索范围：假设最大阶数为15
maxOrder = 15;
mdl_values = zeros(maxOrder,1); % 预分配数组用于存放MDL值

%% 3. 对每个候选阶数p计算MDL值
for p = 1:maxOrder
    model = ar(y, p, 'burg');        % 用Burg算法估计AR(p)模型
    noise_var = model.NoiseVariance; % 取出残差的方差
    mdl_values(p) = N * log(noise_var) + p * log(N);  % MDL公式
end

%% 4. 找出MDL最小值对应的最优阶
[~, bestOrder] = min(mdl_values);

%% 5. 绘制MDL曲线
figure;
plot(1:maxOrder, mdl_values, 'o-', 'LineWidth', 1.5, 'MarkerSize', 7);
hold on;
plot(bestOrder, mdl_values(bestOrder), 'rp', 'MarkerFaceColor', 'r', 'MarkerSize', 12);
text(bestOrder + 0.4, mdl_values(bestOrder), ...
    ['最优阶 = ', num2str(bestOrder)], 'FontSize', 12, 'Color', 'r');

xlabel('模型阶数 p');
ylabel('MDL 值');
title('基于 MDL 准则的模型自动选阶');
grid on;

3、代码讲解

1、生成信号，这一段用 filter 函数生成一个自回归信号（AR(3)）：

2、搜索阶数范围，我们不知道信号真正的阶数是多少，所以从 1 阶、2 阶、…、15 阶都试一遍。每次都会计算一个 MDL 值，然后找最小的。

3、计算每个阶的MDL

ar(y, p, 'burg')：使用 Burg 方法 估计 AR(p) 模型，Burg 方法的优点是稳定、抗噪声好，它不会生成不稳定多项式（保证模型合理）。

noise_var：就是模型残差的方差（预测误差大小）；

mdl_values(p) 按 MDL 公式计算。

4、选出最优阶，绘图显示。

找到 MDL 曲线的最小值所对应的阶数，这个 bestOrder 就是自动选出的最佳模型阶。曲线的最低点就是最佳阶。因为 MDL 惩罚更强，更偏好简单模型。

4、结果分析

运行后，你会看到一张折线图：

横轴：模型阶数 p
纵轴：对应的MDL值

MDL 曲线通常曲线略靠左，选出的阶数略小一点，因为MDL惩罚更强，更偏好简单模型。

5、优缺点

优点	缺点
理论坚实	可能欠拟合
惩罚更合理	若涉及非线性或非高斯模型，需要复杂编码估计过程
模型稳定性好	假设噪声独立同分布、模型族正确

6.3.5 传输函数准则(Criterion Autoregressive Transfer Function,CAT)

1、背景及介绍

这个准则不像 FPE、AIC、MDL 那样常出现在入门教材里，它是 Akaike 在研究线性动态系统（尤其是ARMAX模型）时提出的一种更广义的信息准则，主要用于在传递函数模型的辨识中自动选阶。我们可以把它理解为：AIC 的“传递函数版本”。

当我们用信号去辨识系统时，往往得到的是类似这样的传递函数：

也就是：

在实际操作中，并不知道分子多项式的阶数是多少，也不知道分母的阶数是多少。

CAT 准则就是用来自动选择这些阶数的。
它的思想跟 AIC、MDL 是一脉相承的——同样在平衡“拟合误差”和“模型复杂度”之间，但它特别针对传递函数模型的结构辨识问题。

如果要描述一个系统的输入–输出关系。系统越复杂（阶数越高），你就能更准确地拟合数据，但模型参数会越来越多。应该该取多高的阶数，既能很好地拟合系统，又不会让模型过于复杂。CAT从信息论的角度出发，把传递函数的参数个数当成“描述系统所需的信息量”，并引入一个惩罚项去控制模型复杂度。

2、基本形式

CAT 准则的形式与 AIC 类似，只不过考虑了传递函数模型的特点（输入–输出、噪声模型同时存在）：

其中：

N：样本点数
：模型残差方差（拟合误差大小）
：系统分母的阶数（反馈部分）
：系统分子的阶数（输入部分）
：输入输出之间的延迟阶数（dead time）

可以看到，CAT其实是 AIC 的自然扩展版，只是把“参数个数”改成了传递函数模型的三个关键维度。

3、matlab实现

1、完整代码

%% =====================================================
%  基于 CAT（传输函数准则）的自动选阶实战示例
%  作者：木夕水
%  说明：通过直观的热力图和误差趋势展示选阶结果
% =====================================================
clear; clc; close all;

%% 1. 构造一个已知系统：ARX(2,2,1)
N = 20000;
rng(2);                                   % 固定随机数以保证复现性
u = randn(N,1);                           % 输入信号（白噪声）
a_true = [1, -0.5, 0.4];                  % 分母系数 A(z)
b_true = [0, 0.8, -0.3];                  % 分子系数 B(z)（延迟1）
e = 0.1 * randn(N,1);                     % 噪声信号
y = filter(b_true, a_true, u) + e;        % 系统输出
data = iddata(y, u, 1);                   % 转为 System ID 格式

%% 2. 设置搜索范围
na_max = 6;       % 分母阶数上限
nb_max = 6;       % 分子阶数上限
nk = 1;           % 输入延迟固定为1
CAT_val = zeros(na_max, nb_max);  % 存放 CAT 值

%% 3. 遍历每一种模型结构，计算 CAT 值
for na = 1:na_max
    for nb = 1:nb_max
        % 估计ARX模型
        model = arx(data, [na nb nk]);
        
        % 用“一步预测误差”计算残差方差（避免过拟合）
        y_pred = predict(model, data, 1);
        residual = y - y_pred.y;
        sigma2 = var(residual);
        
        % 计算 CAT 准则
        N = length(y);
        CAT_val(na, nb) = N * log(sigma2) + 2 * (na + nb + nk);
    end
end

%% 4. 找出 CAT 最小对应的模型阶数
[min_val, idx] = min(CAT_val(:));
[best_na, best_nb] = ind2sub(size(CAT_val), idx);

fprintf('最优模型阶数：na = %d, nb = %d, nk = %d\n', best_na, best_nb, nk);
fprintf('真实模型阶数：na = %d, nb = %d, nk = %d\n', ...
    length(a_true)-1, length(b_true)-1, 1);

%% 5. 可视化结果一：热力图
figure;
imagesc(1:nb_max, 1:na_max, CAT_val);
set(gca, 'YDir', 'normal');
colorbar;
xlabel('分子阶数 nb');
ylabel('分母阶数 na');
title('基于 CAT 准则的模型阶数选择');
text(best_nb, best_na, '★', 'Color','r','FontSize',20, ...
     'HorizontalAlignment','center','VerticalAlignment','middle');
grid on;

%% 6. 可视化结果二：CAT 值随阶数变化的趋势
figure;
plot(1:na_max, CAT_val(:, best_nb), 'o-', 'LineWidth', 1.5, 'MarkerSize', 7);
hold on;
plot(best_na, CAT_val(best_na, best_nb), 'rp', 'MarkerFaceColor', 'r', 'MarkerSize', 12);
xlabel('分母阶数 na');
ylabel('CAT 值');
title(['CAT随na变化趋势（nb固定为 ', num2str(best_nb), '）']);
grid on;