matlab.rar_最小二乘法

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.584038000000000,-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