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在计算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.3AICsensor 2sensor 4sensor 6sensor 8sensor 10sensor 12sensor 14sensor 16 -3.35-3.4-3.45-3.5-3.55-3.6 020406080AR order100120140

AR模型预测

sensor 2. (1-step pred) 864sensor 2 (m/s2)Measuredm; fit: 81.01%20-2-4-6 11.051.11.151.21.25

AR误差计算

e@sensor 221m/s20-1-211.051.11.151.2Time 1.251.31.351.4

Correlation function of residuals. Output sensor 21.210.80.60.40.20-0.20510lag152025

附录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.