gls.m

最小二乘法中的递推最小二乘法

function[n_die,theta,ef]=GLS(N,n,In,Out)
p=2;
%n_die为迭代的次数
%ef迭代过程中的误差
ee=0.000001;
x=rand(N,2*n+1);
for i=1:N
    for j=1:n
        if(i-j>0)
        x(i,j)=-Out(i-j);
        else x(i,j)=0;
        end
    end
    for k=(n+1):(2*n+1)
        if(i-k+n+1>0)
        x(i,k)=In(i-k+n+1);
        else x(i,k)=0;
        end
    end
end
theta=inv(x'*x)*x'*Out';
%广义二乘
for m=1:N
    %生成误差序列
    for i=n+1:N
        e(i-n)=Out(i);
        for j=1:n
            e(i-n)=e(i-n)+theta(j)*Out(i-j);
        end
        for k=(n+1):(2*n+1)
            e(i-n)=e(i-n)-theta(k)*In(i-k+n+1);
        end
    end
    %对误差序列进行评估
    for i=1:N-n;%求E
        for j=1:p
            if(i-j>0)
            E(i,j)=-e(i-j);
            else E(i,j)=0;
            end           
        end
    end
    f=inv(E'*E)*E'*e';
    for k=n+1:N
        for j=1:p
            Out(k)=Out(k)+f(j)*Out(k-j);
            In(k)=In(k)+f(j)*In(k-j);
        end
    end
    %最小二乘法
    x=rand(N,2*n+1);
for i=1:N
    for j=1:n
        if(i-j>0)
        x(i,j)=-Out(i-j);
        else x(i,j)=0;
        end
    end
    for k=(n+1):(2*n+1)
        if(i-k+n+1>0)
        x(i,k)=In(i-k+n+1);
        else x(i,k)=0;
        end
    end
end
thetaN=theta;
theta=inv(x'*x)*x'*Out';
    ef(m)=norm(theta-thetaN)/norm(thetaN);
    if(ef(m)<ee)
        n_die=m;
        break;
    end
    n_die=N;
end