Matlab_AR模型阶数确定

自回归（AR）模型

理论模型

自回归（AutoRegressive, AR）模型又称为时间序列模型，数学表达式为

AR:y(t) a1y(t 1) ... anay(t na) e(t)

其中，e(t)为均值为0，方差为某值的白噪声信号。

Matlab Toolbox

研究表明，采用Yule-Walker方法可得到优化的AR模型[1]，故采用aryule程序估计模型参数。

[m,refl] = ar(y,n,approach,window)

模型阶数的确定

有几种方法来确定。如Shin提出基于SVD的方法，而AIC和FPE方法是目前应用最广泛的方法。若计算出的AIC较小，例如小于-20，则该误差可能对应于损失函数的10-10级别，则这时阶次可以看成是系统合适的阶次。

am = aic(model1,model2,...)

fp = fpe(Model1,Model2,Model3,...)

AR预测

yp = predict(m,y,k)

m表示预测模型；y为实际输出；k预测区间；yp为预测输出。 y(1),y(2),...,y(t k 1),y(t k),...,y(t 2),y(t 1),y(t)

当k<Inf时，yp(t)为模型m与y(1,2,…t-k)的预测值；当k=Inf时，yp(t)为模型m的纯仿真值；默认情况下，k=1。

在计算AR模型预测时，k应取1，原因参照AR模型理论公式。

compare(y,m,k)

[yh,fit,x0] = compare(y,m,k)

Compare的预测原理与predict相同，但其对预测进行了比较。

||y yh||

fit 100 1 ||y ||

AR误差

e = pe(m,data)

pe误差计算。采用yh=predict(m,data,1)进行预测，然后计算误差e=data-yh;

[e,r]= resid(m,data,mode,lags); resid(r)

resid计算并检验误差。采用pe计算误差；在无输出的情况下，绘出误差图，误差曲线应足够小，黄色区域为99%的置信区间，误差曲线在该区域内表明通过检验。

Matlab练习

确定模型阶数

采用ASCE benchmark模型120DOF，选取y方向的响应，共8个。首先，对响应数据进行标准化处理；其次，将标准化处理后的数据建立AR模型；最后，确定合适的模型阶次，通过选取一系列阶数，分别计算对应的AIC值，从图中可以看出，阶次80以后的AIC值变化不大，因此，合适的阶次选择为80。

-3.1-3.15-3.2-3.25-3.3AIC

-3.35-3.4-3.45-3.5-3.55-3.60

20

40

6080AR order

100

120

140

AR模型预测

sensor 2. (1-step pred)

864sensor 2 (m/s2)

20-2-4-6

1

1.05

1.1

1.15

1.2

1.25

AR误差计算

e@sensor 2

2

1

m/s2

-1

-2

11.051.11.15

1.2Time

1.251.31.351.4

Correlation function of residuals. Output sensor 2

0510

lag

152025

附录MATLAB 代码

（1）

% AR model order clc;clear;

addpath(genpath(pwd),1);

data = load('Case2Damage0_1_0.01_0.001_40_10_150_123_1.mat'); dofy = [2,4,6,8,10,12,14,16];

x = zscore(data.acc(:,dofy)); order = [10:10:130]; for i = 1:length(order) for j =1:size(x,2)

m = ar(x(:,j),order(i),'yw'); am(i,j) = aic(m); end end

plot(order,am,'LineWidth',2,'Marker','*') xlabel('AR order'),ylabel('AIC')

legend('sensor 2','sensor 4','sensor 6','sensor 8',... 'sensor 10','sensor 12','sensor 14','sensor 16')

（2）

clc;clear;

addpath(genpath(pwd),1);

data = load('Case2Damage0_1_0.01_0.001_40_10_150_123_1.mat'); dofy = [2;4;6;8;10;12;14;16]; order = 80;

ts = zscore(data.acc(:,dofy)); k = 1;

name = {'sensor 2','sensor 4','sensor 6','sensor 8',... 'sensor 10','sensor 12','sensor 14','sensor 16'}; Ounit = repmat({'m/s^2'},length(dofy),1);

ts = iddata(zscore(data.acc(:,dofy)),[],data.dt,... 'OutputName',name,'OutputUnit',Ounit,... 'Name','ASCE-benchmark 120DOF');

m = ar(ts(:,k),order,'yw');

compare(ts(:,k),m,'r-.',1,1000:1300);

（3）

resid(m,ts(:,k),'corr',25)

[e,r] = resid(m,ts(:,k),'corr',40);

resid(r)

plot(e(1000:1300))

参考文献

Matlab help

[1] DA SILVA S, DIAS J NIOR M, LOPES JUNIOR V. Damage detection in a benchmark structure using AR-ARX models and statistical pattern recognition [J]. J Braz Soc Mech Sci Eng, 2007, 29(2): 174-84.

本文来源：https://www.bwwdw.com/article/fir4.html