递推最小二乘法与广义最小二乘法matlab源程序

clear;
clc;
nn = normrnd(0,sqrt(0.5),1,31)';
uk=[1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.958 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 1.144 1.177 -0.390];
yk(1)=0;
yk(2)=0;
for i=1:29;
yk(i+2)=-1.642*yk(i+1)-0.715*yk(i)+0.39*uk(i+1)+0.35*uk(i)+nn(i+2)+1.642*nn(i+1)+0.715*nn(i);
end;
for i=1:29;
A(i,:)=[-yk(i+1) -yk(i) uk(i+1) uk(i)];
end
siit=inv(A'*A)*A'*(yk(3:31)+nn(2:30)')';
e(1)=yk(1);
e(2)=yk(2)+siit(1)*yk(1)-siit(3)*uk(1);
for i=3:31;
e(i)=yk(i)+siit(1)*yk(i-1)+siit(2)*yk(i-2)-siit(3)*uk(i-1)-siit(4)*uk(i-2);
end
for i=1:29;
fai(i,:)=[-e(i+1) -e(i)];
end
f=inv(fai'*fai)*fai'*e(3:31)';
for i=3:31;
yk(i)=yk(i)+f(1)*yk(i-1)+f(2)*yk(i-2);
end
yk(2)=yk(2)+f(1)*yk(1);
for i=3:30;
uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);
end
uk(2)=uk(2)+f(1)*uk(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:30
for i=1:29;
A(i,:)=[-yk(i+1) -yk(i) uk(i+1) uk(i)];
end
siit=inv(A'*A)*A'*yk(3:31)';
e(1)=yk(1);
e(2)=yk(2)+siit(1)*(yk(1))-siit(3)*uk(1);
for i=3:31;
e(i)=yk(i)+siit(1)*(yk(i-1))+siit(2)*(yk(i-2))-siit(3)*uk(i-1)-siit(4)*uk(i-2);
end
for i=1:29;
fai(i,:)=[-e(i+1) -e(i)];
end
f=inv(fai'*fai)*fai'*e(3:31)';
k1(j)=f(1);
k2(j)=f(2);
for i=3:31;
yk(i)=yk(i)+f(1)*(yk(i-1))+f(2)*(yk(i-2));
end
yk(2)=yk(2)+f(1)*yk(1);
for i=3:30
uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);
end
uk(2)=uk(2)+f(1)*uk(1);
end
siit';
if(siit(1)>0&&siit(2)>0&&siit(3)>0&&siit(4)>0)

display('方差为0.0001时，利用广义最小二乘法获得参数结果：')

siit'
end

clear;
yk=[0
0
0.47984
-1.02448
0.443946
-0.962888
1.23318
-0.584038
1.09394
-0.583966
0.564678
-0.731735
0.778365
-0.488544
0.599589
-0.878625
0.217687
-0.0144152
-0.590687
1.16106
-0.527737
0.628383
0.152147
-0.1108
0.605338
0.214697
-0.320849
0.601426
-0.000474461
0.430189
-0.446182];
yk=yk';
uk=[1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.958 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 1.144 1.177 -0.390];
for i=1:29;
faik(i,:)=[-yk(i+1) -yk(i) uk(i+1) uk(i)];
end
%faik(3,:)'*faik(3,:)
%faik(1,:)'*faik(1,:)-faik(1,:)'*faik(1,:)*faik(2,:)'/(faik(2,:)*faik(1,:)'*faik(1,:)*faik(2,:)'+1)*faik(2,:)*faik(1,:)'*faik(1,:)

%faik(1,:)'*faik(1,:)*faik(1,:)'*[yk(3)]+faik(1,:)'*faik(1,:)*faik(2,:)'/(faik(2,:)*faik(1,:)'*faik(1,:)*faik(2,:)'+1)*(yk(4)-faik(2,:)*faik(1,:)'*faik(1,:)*faik(1,:)'*[yk(3)])
A=[0,0,0.201000000000000,1.14700000000000;-0.479840000000000,0,-0.787000000000000,0.201000000000000;1.02448000000000,-0.479840000000000,-1.58900000000000,-0.787000000000000;-0.443946000000000,1.02448000000000,-1.05200000000000,-1

.58900000000000;0.962888000000000,-0.443946000000000,0.866000000000000,-1.05200000000000;-1.23318000000000,0.962888000000000,1.15200000000000,0.866000000000000;0.5840380000
00000,-1.23318000000000,1.57300000000000,1.15200000000000;-1.09394000000000,0.584038000000000,0.626000000000000,1.57300000000000;0.583966000000000,-1.09394000000000,0.433000000000000,0.626000000000000;]
B=A;
A=A'*A;
A*faik(10,:)'
A*B'
faik(10,:)*A
A=6*10000000*eye(4);

faik(10,:)*A*faik(10,:)'
kn_1=A*faik(10,:)'/(faik(10,:)*A*faik(10,:)'+1)
%pn_1=A-A*faik(1,:)'/(faik(1,:)*A*faik(1,:)'+1)*faik(1,:)*A
siit=[0.001,0.001,0.001,0.001]';
%siit=A*B'*[0.479840000000000,-1.02448000000000,0.443946000000000,-0.962888000000000,1.23318000000000,-0.584038000000000,1.09394000000000,-0.583966000000000,0.564678000000000;]'+A*faik(10,:)'/(faik(10,:)*A*faik(10,:)'+1)*(yk(12)-faik(10,:)*A*B'*[0.479840000000000,-1.02448000000000,0.443946000000000,-0.962888000000000,1.23318000000000,-0.584038000000000,1.09394000000000,-0.583966000000000,0.564678000000000;]')
for i=1:28;
pn_1=A-A*faik(1+i,:)'/(faik(1+i,:)*A*faik(1+i,:)'+1)*faik(1+i,:)*A;
siit=siit+A*faik(1+i,:)'/(faik(1+i,:)*A*faik(1+i,:)'+1)*(yk(3+i)-faik(1+i,:)*siit)
m(i,:)=siit;
A=pn_1;
end