程序

n=100;

ys=load('C:\\Users\\lenovo\\Desktop\\data.txt'); for p=1:30 for i=1:n-p y(i)=ys(p+i-1); for j=1:p

x(i,j)=ys(p+i-j); end end

fai=inv(x'*x)*(x'*y');

Delta(p)=(y'-x*fai)'*(y'-x*fai)/(n-p); criterion(p,2)=n*log(Delta)+2*p; end

请问上面程序有错误吗？matlab高手知道matlab中有criterion(p,2)这个函数吗？ 注：上述程序使用来ar模型中aic定阶的，希望高手们帮帮忙

x=[1093.2 544.3 868.9 431.2 1209.9 667.4 430.9 1484.3 647.5 370.2 606.2 671 817.6 851.1 965.8 635 1160 794.3 448.4 516.9 486.2 764.4

1

637.9 829.7 449.8 929.7 570.4 387.9 765.8 609.5 622.1 1563 519.3 1191.8 792.6 1018.8 587.8 896.2 474.3 301.1 785.3 483.1 988.9 822.8 664.3 641.7 448.5 513.2 604.5 499.3 896.3 868.1 703.1 439.7 354 824.5 597.6 769.6 465.4 617.5 569.7 668.4 450.8]; x=x'; for i=0:5

for j=0:5

2

spec=garchset('R',i,'M',j,'Display','off');%指定模型的结构 [coeffX,errorsX,LLFX]=garchfit(spec,x);%拟合参数 num=garchcount(coeffX);%计算拟合参数的个数

[aic,bic]=aicbic(LLFX,num,63);%这行matlab老提示错误，不知道错在哪里 fprintf('R=%d,M=%d,AIC=%f,BIC=%f\\n',i,j,aic,bic);%显示计算结果

end

例题：某化工生产过程每2 小时的浓度读数如表13 所示。 表13 化工生产过程浓度数据 17.0 16.6 16.3 16.1 17.1 16.9 16.8 17.4 17.1 17.0 16.7 17.4 17.2 17.4 17.4 17.0 17.3 17.2 17.4 16.8 17.1 17.4 17.4 17.5 17.4 17.6 17.4 17.3 17.0 17.8 17.5 18.1 17.5 17.4 17.4 17.1 17.6 17.7 17.4 17.8 17.6 17.5 16.5 17.8 17.3 17.3 17.1 17.4 16.9 17.3 17.6 16.9 16.7 16.8 16.8 17.2 16.8 17.6 17.2 16.6 17.1 16.9 16.6 18.0 17.2 17.3 17.0 16.9 17.3 16.8 17.3 17.4 17.7 16.8 16.9 17.0 16.9 17.0 16.6 16.7 16.8 16.7 16.4 16.5 16.4 16.6 16.5 16.7 16.4 16.4 16.2 16.4 16.3 16.4 17.0 16.9 17.1 17.1 16.7 16.9 16.5 17.2 16.4 17.0 17.0 16.7 16.2 16.6 16.9 16.5 16.6 16.6 17.0 17.1 17.1 16.7 16.8 16.3 16.6 16.8 16.9 17.1 16.8 17.0 17.2 17.3 17.2 17.3 17.2 17.2 17.5 16.9 16.9 16.9 17.0 16.5 16.7 16.8 16.7 16.7 16.6 16.5 17.0 16.7 16.7 16.9 17.4 17.1 17.0 16.8 17.2 17.2 17.4 17.2 16.9 16.8 17.0 17.4 17.2 17.2 17.1 17.1 17.1 17.4 17.2 16.9 16.9 17.0 16.7 16.9 17.3 17.8 17.8 17.6 17.5 17.0 16.9 17.1 17.2 17.4 17.5 17.9 17.0 17.0 17.0 17.2 17.3 17.4 17.4 17.0 18.0 18.2 17.6 17.8 17.7 17.2 17.4

（1）对这一生产过程建模；

（2）对这一生产过程进行10 步预测。

解 通过计算自相关函数和偏相关函数，确定取d = 1。利用AIC 准则定阶，取ARIMA(3,1,3)模型。计算得

（2）经计算，其10 步预报值见表14。 表 14 10 步预报值

步数 1 2 3 5

预报值17.3611 17.4446 17.5811 17.6441 17.5821

3

4 步数 6 7 8 9 10 预报值17.4644 17.4079 17.4644 17.5759 17.6364 计算的 Matlab 程序如下： clc,clear

a=textread('hua.txt'); %把原始数据按照原来的排列格式存放在纯文本文件hua.txt a=nonzeros(a')'; %按照原来数据的顺序去掉零元素 r11=autocorr(a) %计算自相关函数 r12=parcorr(a) %计算偏相关函数 da=diff(a); %计算1 阶差分

r21=autocorr(da) %计算自相关函数 r22=parcorr(da) %计算偏相关函数

n=length(da); %计算差分后的数据个数 for i=0:3 for j=0:3

spec= garchset('R',i,'M',j,'Display','off'); %指定模型的结构 [coeffX,errorsX,LLFX] = garchfit(spec,da); %拟合参数 num=garchcount(coeffX); %计算拟合参数的个数 %compute Akaike and Bayesian Information Criteria [aic,bic]=aicbic(LLFX,num,n);

fprintf('R=%d,M=%d,AIC=%f,BIC=%f\\n',i,j,aic,bic); %显示计算结果 end end

r=input('输入阶数R＝');m=input('输入阶数M＝');

spec2= garchset('R',r,'M',m,'Display','off'); %指定模型的结构 [coeffX,errorsX,LLFX] = garchfit(spec2,da) %拟合参数

[sigmaForecast,w_Forecast] = garchpred(coeffX,da,10) %计算10 步预报值 x_pred=a(end)+cumsum(w_Forecast) %计算原始数据的10 步预测值

我完全按照原代码运行，使用的是matlabR2012a，结果为 R=0,M=0,AIC=214.909040,BIC=221.465270 R=0,M=1,AIC=142.101300,BIC=151.935644 R=0,M=2,AIC=140.537062,BIC=153.649521 R=0,M=3,AIC=142.061727,BIC=158.452301 R=1,M=0,AIC=177.357950,BIC=187.192294 R=1,M=1,AIC=140.116716,BIC=153.229175 R=1,M=2,AIC=142.888600,BIC=159.279174 R=1,M=3,AIC=143.725801,BIC=163.394489 R=2,M=0,AIC=163.520336,BIC=176.632794 R=2,M=1,AIC=141.783127,BIC=158.173701 R=2,M=2,AIC=143.739875,BIC=163.408563 R=2,M=3,AIC=145.350387,BIC=168.297189 R=3,M=0,AIC=161.615677,BIC=178.006250

4

R=3,M=1,AIC=143.597322,BIC=163.266010 R=3,M=2,AIC=144.969054,BIC=167.915857 R=3,M=3,AIC=142.044630,BIC=168.269548

若按照AIC.BIC准则定阶，应选取AIC.BIC 的最小值，由以上结果看R=1,M=1,AIC=140.116716,BIC=153.229175 才是最小值，并不是答案给出的R=3,M=3啊，这是怎么回事？？

另外，我自己的运行结果是 coeffX =

Comment: 'Mean: ARMAX(3,3,0); Variance: GARCH(0,0) ' Distribution: 'Gaussian' R: 3 M: 3 C: 0.0035

AR: [-0.0101 -0.9030 0.1848] MA: [-0.6932 0.8359 -0.8797] 和答案给出的参数值也不一样，怎么回事？？

VarianceModel: 'GARCH' K: 0.1114 Display: 'off'

errorsX =

Comment: 'Mean: ARMAX(3,3,0); Variance: GARCH(0,0) ' Distribution: 'Gaussian' R: 3 M: 3 C: 0.0070

AR: [0.0834 0.0353 0.0836] MA: [0.0494 0.0239 0.0489] VarianceModel: 'GARCH' K: 0.0124

LLFX =

-63.0223

sigmaForecast =

0.3337 0.3337 0.3337 0.3337 0.3337 0.3337 0.3337

5

0.3337 0.3337 0.3337

w_Forecast =

0.5151 -0.0198 -0.0116 0.1167 0.0091 -0.1041 0.0179 0.0990 -0.0329 -0.0823

x_pred =

17.3151 17.2954 17.2838 17.4005 17.4096 17.3055 17.3234 17.4224 17.3895 17.3072

预测值和给出的答案也不一样，怎么回事呢？ 请高手指点，急求答案！！谢谢

很早之前发过一个帖子，关于AR模型阶次估计的，《MATLAB中AR模型功率谱估计中AR阶次估计的实现》，后来很多人问我要这几个程序，其实帖子里面已经将程序上传了，里面的程序是对原始的工具箱函数进行了修改，加入了定阶准则这一步，然后返回模型参数以及模型阶次。

定阶准则定出的阶次一般也只有参考意义，对于数据较短的序列，定阶准则给出的阶次比较小（胡光书《数字信号处理-理论、算法与实现》）。MATLAB的信号处理工具箱中只有参数求解的函数，没有给出模型阶次求解的功能。高阶谱分析工具箱（HOSA）里面的arorder函数可以进行定阶，但是给出的阶次很小（我是用过几次，每次大概2-3），相比定阶准则给出的阶次来说，和实际使用的AR模型相差更大，使用诸如FPE等定阶准则给出的阶次还是有一定的参考意义。

这里我们不再使用通用的signal工具箱，直接将burg算法编写出来，中间加入定阶准则的内容，可以直接用来进行参数求解以及功率谱的求解，效率比使用工具箱函数要高一些（因为工具箱函数为了适应各种参数要求，里面加入了太多的判断内容）。

关于burg算法求解AR参数的过程，这里不再详细给出了，大家可以参考随机信号处理的相关教材。步骤就是一个递推的过程，直到某一个阶次p为止。

游客，如果您要查看本帖隐藏内容请回复

6

上面的程序为burg算法求解AR模型阶次，以及使用FPE和AIC定阶，输入参数为3个，前两个固定为信号和采样率，后一个如果为数值型，那么就是模型阶次，也就是不用定阶，直接使用该阶次求解模型参数。如果最后一个参数为字串，FPE或者AIC，那么就是定阶准则，需要在程序中定阶。模型参数求解过程参考了arburg函数，工具箱函数有很大的参考意义，毕竟是各个领域的大师级人物写的，效率发挥到极致。

下面为一个测试： 01.

02.randn('state', 1);

03.x = randn(100, 1);

04.y = filter(1, [1 1/2 1/3 1/4 1/5], x);

05.pburg(y, 4, [], 1000);

06.[psd, f, p] = myBurg(y, 1000, 'FPE');

07.figure;

08.plot(f, 10*log10(psd)); 09.

复制代码

使用FPE定阶给出的阶次为5，功率谱曲线如下（前一个为pburg工具箱函数，后面为自己编写的）

近在做一个matlab gui 程序，是用时间序列法预测数据的，有一个m文件一个fig文件，在matlab中运行没有问题，但是生成exe文件后，可执行文件可以运行，但是会在DOS页面上报错，算法部分不能实现（就是不能返回预测值），m文件中主要用到armax函数等。

我的qq：980246573 欢迎讨论！

算法部分主要代码如下： function sjxlxiugai(wenjianlujing)

A=xlsread(wenjianlujing,'sheet1','A1：A720'); %用于ARIMA建模的彭水电站3月1日-3月30日每小时的入库流量数据a=diff(A);

7

%计算原始序列的一阶差分%%模型的建立 X1=iddata(a); %将a转化为matlab接受的格式 a1=[a;zeros(24,1)];

x1=iddata(a1);test = [];

for p = 1:10 %自回归对应偏自相关函数（PACF） for q = 1:10 %移动平均对应自相关函数（ACF）

m=armax(X1,[p,q]); %利用不同的p、q建立arma模型 AIC = aic(m); %求得建立的arma模型的AIC值 test = [test;p q AIC]; %求得所有p、q组合的AIC值 endendfor k = 1:size(test,1)

if test(k,3) == min(test(:,3)) %选择AIC值最小对应的p、q值 p_test = test(k,1); q_test = test(k,2); break; end end

%拟合过程

m1 = armax(X1,[p_test q_test]) ;

%用AIC最小对应的p,q建立ARMA模型%%模型的预测与分析 p=predict(m1,x1,1);

p = get(p,'outputData'); %将iddata的数据类型转换为double类型的数据%差分还原 y1=[A;zeros(24,1)];y=zeros(24,1); for i=721:744

y1(i)= p(i-1)+y1(i-1); y(i-720)=y1(i);end %figure(1)

%matlab中绘出y的图像%plot(y,'-*b') %xlabel('时间/1h');

%ylabel('入库流量/m3/s');

%输出未来一天每小时的入库流量

xlswrite('C:\\Documents and Settings\\Administrator\\桌面\\预测结果.xls',y,'时间序列法预测结果');

end

各位大牛们，希望能帮个忙啊，最近在看chi2gof函数这块，可是不知道为什么返回的检验值P总为NaN，希望大牛们能帮我解决下，谢谢！

以下是我的程序： clc,clear

a=[17.98 16.92 16.78 16.85 16.85 16.97 16.92 17.54 16.49 17.08 17.26 17.16 16.90]； a=nonzeros(a')'; %按照原来数据的顺序去掉零元素 da=diff(a);%消除趋势性的差分运算

8

daa=diff(da);%消除趋势性的差分运算 w=length(daa); %计算差分后的数据个数 n=20;%预报的数据个数 for i=0:3

for j=0:3

spec= garchset('R',i,'M',j,'Display','off'); %指定模型的结构 [coeffX,errorsX,LLFX] = garchfit(spec,daa); %拟合参数 num=garchcount(coeffX); %计算拟合参数的个数 %compute Akaike and Bayesian Information Criteria [aic,bic]=aicbic(LLFX,num,w);

fprintf('R=%d,M=%d,AIC=%f,BIC=%f\\n',i,j,aic,bic); %显示计算结果 end end

r=input('输入阶数R＝');m=input('输入阶数M＝');

spec2= garchset('R',r,'M',m,'Display','off'); %指定模型的结构 [coeffX,errorsX,LLFX] = garchfit(spec2,daa) %拟合参数

e(1)=0; %残差赋初值 for s=2:13

e(s)=a(s)-coeffX.AR*a(s-1); end

[h,p,stats] = chi2gof(e,'nbins',1,'cdf',@(z)normcdf(z,0,std(e)),'nparams',1)

返回的结果如下：coeffX =

Comment: 'Mean: ARMAX(1,0,0); Variance: GARCH(0,0) ' 'Gaussian' R: 1 C: 0.1308 VarianceModel: 'GARCH' K: 0.3994 Display: 'off'

errorsX =

Comment: 'Mean: ARMAX(1,0,0); Variance: GARCH(0,0) ' 'Gaussian' R: 1 C: 0.2310 VarianceModel: 'GARCH' K: 0.2325

LLFX =

-10.5608 h = 0 p =

NaN

stats =

chi2stat: 0 df: 0 edges: [4.9407e-324 28.0419] 9

Distribution: AR: -0.6186 Distribution: AR: 0.2954 O: 13

E: 13

y=xlsread('E:\\PM2.5 1.xls','c1:c225')

Data=y; %共300个数据 SourceData=Data(1:225,1); %前250个训练集 step=50; %后50个测试 TempData=SourceData;

TempData=detrend(TempData);%去趋势线 TrendData=SourceData-TempData;%趋势函数 %--------差分，平稳化时间序列--------- H=adftest(TempData); difftime=0;

SaveDiffData=[]; while ~H

SaveDiffData=[SaveDiffData,TempData(1,1)];

TempData=diff(TempData);%差分，平稳化时间序列 difftime=difftime+1;%差分次数

H=adftest(TempData);-f检验，判断时间序列是否平稳化 end

%---------模型定阶或识别-------------- u = iddata(TempData); test = [];

for p = 1:5 %自回归对应PACF,给定滞后长度上限p和q，一般取为T/10、ln(T)或T^(1/2),这里取T/10=12

for q = 1:5 %移动平均对应ACF m = armax(u,[p q]);

AIC = aic(m); %armax(p,q),计算AIC test = [test;p q AIC]; end end

for k = 1:size(test,1)

if test(k,3) == min(test(:,3)) %选择AIC值最小的模型 p_test = test(k,1);

10

q_test = test(k,2); break; end end

%------1阶预测-----------------

TempData=[TempData;zeros(step,1)]; n=iddata(TempData);

m = armax(u,[p_test q_test]);

%m

= armax(u(1:ls),[p_test

q_test]); %armax(p,q),[p_test q_test]对应AIC值最小，自动回归滑动平均模型

P1=predict(m,n,1); PreR=P1.OutputData; PreR=PreR';

%----------还原差分----------------- if size(SaveDiffData,2)~=0

for index=size(SaveDiffData,2):-1:1

PreR=cumsum([SaveDiffData(index),PreR]); end end

%-------------------预测趋势并返回结果---------------- mp1=polyfit([1:size(TrendData',2)],TrendData',1); xt=[];

for j=1:step

xt=[xt,size(TrendData',2)+j]; end

TrendResult=polyval(mp1,xt);

PreData=TrendResult+PreR(size(SourceData',2)+1:size(PreR,2)); tempx=[TrendData',TrendResult]+PreR; % tempx为预测结果 plot(tempx,'r'); hold on

plot(Data,'b');

function arma=armaa(x) clc

GG=load('d:/gr.txt'); KK=GG; GG=dtrend(GG);

aic=autocorr(GG);%自相关 paic=parcorr(GG);%偏相关

11

mm=-aic; nn=-paic; paic(1)=[]; aic(1)=[]; nn(1)=1; for i=2:13, for j=1:i, dd(j)=nn(j); end

c=filter(dd,1,GG); for op=1:i-1, c(op)=[]; end

canca(i-1)=var(c); end for i=1:12,

BIC(i)=log(canca(i))+1/length(KK)*log(length(KK)); AIC(i)=log(canca(i))+2*i/length(KK); end aic paic

canca %残差平方和 AIC BIC

输出结果：

aic =

12

-0.5037 0.1220 -0.2108 0.0796 0.0187 0.1169 -0.2162 0.2517 -0.1932 0.0559 -0.1036 0.0123 0.2180 -0.0591 -0.1431 0.0697 -0.0505 0.0115 0.0312 -0.0766 paic = -0.5041 -0.1781 -0.3551 -0.2831 -0.1772 0.0360

13

-0.2209 0.1458 0.0688 0.0009 -0.1723 -0.2912 0.1872 0.1437 -0.2386 0.1059 0.2155 -0.3839 -0.2884 -0.2151 canca = 1.0e+003 *

Columns 1 through 11

2.5172 2.4685 2.3110 8 2.1040 Column 12 2.1688 AIC =

Columns 1 through 11

14

2.2903 2.3482

2.3340 2.2696 2.2007 1.9902 2.022 7.8660 7.8815 7.8507 7.8768 7.9368 7.9659 7.9730 7.9772 7.9118 7.9631 8.0376 Column 12 8.1030 BIC =

Columns 1 through 11

7.9019 7.8823 7.8164 7.8074 7.8323 2 7.7225 Column 12 7.7529

对于acf和pacf，跟EVIEWS输出结果一样。

15

7.8263 7.7983 7.7675 7.6669 7.683

=IF(INDEX(圆通!$A$3:$L$308,MATCH(价格计算!A2,圆通!$A$3:$A$308,0),IF(B2=\不限\圆通!$A$3:$L$308,MATCH(价格计算!A2,圆通!$A$3:$A$308,0),IF(B2=\不限\

=IF(G2<=1,(F2*G2+13)/G2,IF(H2<=12,(F2*G2+11*CEILING(H2,1)+2)/G2,IF(G2<=25,(F2*G2+132.72)/G2,(F2*G2+CEILING(H2,1)*1.84+16.72+CEILING(INT(H2/25),1)*70)/G2)))

=IF(0

??+??????×

??

????

+????+????.????,??>25???? ?????? ?????? ≤?????? ??????????

??+??????× ???? +????×?????? ???? +????.????,??>25???? ?????? ??????≤

?????? ???? ≤???????? ??+??????× ???? +??????.????+????.????,??>25???? ??????????????≤ ?????? ???? ≤????????

16

????

??

??=

????

??+??????× +????+????.????,??>25???? ?????? ?????? ≤??????

????????

=IF(INDEX(顺丰!$A$3:$L$261,MATCH(价格计算!A2,顺丰!$A$3:$A$261,0),IF(B2=\不限\顺丰!$A$3:$L$261,MATCH(价格计算!A2,顺丰!$A$3:$A$261,0),IF(B2=\不限\

??+????,??≤1??????+??????+??,1????≤??≤12??????+??????.????,????????≤??≤????????

17

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