ce6_6.m

生成一组带有高斯噪声的样本,分别用一阶,二阶,三阶的最小二乘估计方法进行拟合,然后分别用AIC,MDL,FPE,CAT四种评测模型对其性能进行比较,得到最优的拟合模型.

clear all
t = -9:9;
z = zeros(1,19);
v = randn(1,19);
for i = 1:19
    z(i) = 10 + t(i) + 0.1*t(i)^2 + v(i);
end
H1 = zeros(19,2);        %第i行表示H(ti)
H2 = zeros(19,3);
H3 = zeros(19,4);
H1(:,1) = 1;
H1(:,2) = t;
H2(:,1) = 1;
H2(:,2) = t;
H2(:,3) = (t.*t/2)';
H3(:,1) = 1;
H3(:,2) = t;
H3(:,3) = (t.*t/2)';
H3(:,4) = (t.*t.*t)';

temp1 = zeros(2,2); temp2 = zeros(2,1);
for i = 1:19
    temp1 = temp1 + H1(i,:)'*H1(i,:);
    temp2 = temp2 + H1(i,:)'*z(i);
end
X1 = (temp1)^(-1)*temp2;  % X1 = [a0,a1]

temp1 = zeros(3,3); temp2 = zeros(3,1);
for i = 1:19
    temp1 = temp1 + H2(i,:)'*H2(i,:);
    temp2 = temp2 + H2(i,:)'*z(i);
end
X2 = (temp1)^(-1)*temp2;   % X2 = [a0,a1,a2]

temp1 = zeros(4,4); temp2 = zeros(4,1);
for i = 1:19
    temp1 = temp1 + H3(i,:)'*H3(i,:);
    temp2 = temp2 + H3(i,:)'*z(i);
end
X3 = (temp1)^(-1)*temp2;   % X3 = [a0,a1,a2,a3]
plot(t,z,'*R')
hold
i = -9:0.01:9;
z1 = X1(1)+X1(2)*i;
z2 = X2(1)+X2(2)*i+X2(3)*(i.*i)/2
z3 = X3(1)+X3(2)*i+X3(3)*(i.*i)/2+X3(4)*(i.*i.*i)
plot(i,z1)
figure
plot(t,z,'*R');hold
plot(i,z2)
figure
plot(t,z,'*R');hold
plot(i,z3)

error1 = 0;        %calculate prediction error power for m=1
for i = 1:19
    error1 = error1 + (z(i)-H1(1,:)*X1)'*(z(i)-H1(1,:)*X1);
end
error1 = error1/19;

error2 = 0;        %calculate prediction error power for m=2
for i = 1:19
    error2 = error2 + (z(i)-H2(1,:)*X2)'*(z(i)-H2(1,:)*X2);
end
error2 = error2/19;

error3 = 0;        %calculate prediction error power for m=3
for i = 1:19
    error3 = error3 + (z(i)-H3(1,:)*X3)'*(z(i)-H3(1,:)*X3);
end
error3 = error3/19;

AIC1 = 19*log(error1) + 2*1;
AIC2 = 19*log(error2) + 2*2;
AIC3 = 19*log(error3) + 2*3;

MDL1 = 19*log(error1) + 1*log(19);
MDL2 = 19*log(error2) + 2*log(19);
MDL3 = 19*log(error3) + 3*log(19);

FPE1 = error1*(19+1)/(19-1);
FPE2 = error2*(19+2)/(19-2);
FPE3 = error3*(19+3)/(19-3);

CAT1 = (error1*19/(19-1))^(-1) /19 - (error1*19/(19-1))^(-1);
CAT2 = ((error1*19/(19-1))^(-1)+(error2*19/(19-2))^(-1))/19 - (error2*19/(19-2))^(-1);
CAT3 = ((error1*19/(19-1))^(-1)+(error2*19/(19-2))^(-1)+(error3*19/(19-3))^(-1))/19 - (error3*19/(19-3))^(-1);