clc
clear all
filein = fopen('gaocheng.txt');
indat = fscanf(filein,'%f',[36,inf]); %fscanf(fopen打开文件的返回值，格式，大小[m n]的向量)
emdat = indat;
cindat = sym('x', size(emdat));
for p =1:5 %开始循环 p（AR的阶数）从1到7
% % while p>0
[A,b,Ca,Cb]= composeAb(emdat,cindat,p); %调用自定义函数得到 A是系数矩阵 b是观测矩阵 Ca和Cb分别是字符串形式的A和b
[rw,cu]=size(A); %返回矩阵A的尺寸信息， 并存储在rw、cu中。其中rw中存储的是行数，cu中存储的是列数。
DX = eye(rw); %把rw*rw的单位矩阵 赋值给 DX
EDX =DX; %把矩阵DX 赋值给 EDX
EPX =inv(DX); %把矩阵DX的逆矩阵 赋值给 EPX相当于Qy的逆矩阵
lsp = (A'*inv(EPX)*A)\(A'*inv(EPX)*b); %就是迭代算法里的计算参数x(0)
lsv = A*lsp-b;
lsSSR = lsv'*inv(EDX)*lsv; %残差平方和，用来判断阶数
lsSSR = lsv'*inv(EDX)*lsv;
AIC(p) = 2*p+(rw+p)*log(lsSSR/(rw+p));
end
AIC
p=4;
[A,b,Ca,Cb]= composeAb(emdat,cindat,p);
[rw,cu]=size(A); %返回矩阵A的尺寸信息， 并存储在rw、cu中。其中rw中存储的是行数，cu中存储的是列数。
DX = eye(rw); %把rw*rw的单位矩阵 赋值给 DX
EDX =DX; %把矩阵DX 赋值给 EDX
EPX =inv(DX); %把矩阵DX的逆矩阵 赋值给 EPX相当于Qy的逆矩阵
lsp = (A'*inv(EPX)*A)\(A'*inv(EPX)*b); %就是迭代算法里的计算参数x(0)
lsv = A*lsp-b;
H = double(jacobian(reshape(Ca,rw*cu,1),cindat));
Q=H*emdat*H'
