irls.asv

基于M估计的迭代最小二乘算法

function [B] = irls(X,Y,type)
%IRLS   Iteratively reweighted least squares.
%   [B] = irls(X,Y,type) produce a iteratively reweighted lease square
%   regression of Y, on X, based on the influence funtion of the third
%   parameter 'type'. The returning vector B is the regression parameters.
%   Possible values for 'type' ant the influence funtion represented are
%   listed below:
%   type = 1: Huber's t function, t=2
%   type = 2: Ramsay's Ea function, a=0.3
%   type = 3: Andrews's wave funtion 
%   type = 4: Hampel's 17A funtion
[m,n]=size(X);
[m2,n2]=size(Y);
if m~=m2
     error('Column number of X and Y should be equal')
end
W0=zeros(m);
B0=pinv(X'*X)*X'Y;
e0=Y-X*B0;
s=median(abs(e0-median(e0)))/0.6745;
for i=1:m
    if type==1
        W0(i,i)=Huber_weight(e0(i)/s);
    else if type==2
            W0(i,i)=Ramsay_weight(e0(i)/s);
        else if type==3
                W0(i,i)=Andrews_weight(e0(i)/s);
            else if type==4
                    W0(i,i)=Hampel_weight(e0(i)/s);
                else 
                    error('The value specified by the last parameter is wrong')
                end
            end
        end
    end
end
iterate=true;
convergence_value=0.001;
iteration_times=10;
count_iter=0;
W=W0;
Residuals=zeros(m,iteration_times);
Weights=zeros(m,iteration_times);
while(iterate==true&&count_iter<=iteration_times)
    B=pinv(X'*W*X)*X'*W*Y;
    e=Y-X*B;
    W_pervious=W;
    for i=1:m
        if type==1
            W(i,i)=Huber_weight(e(i)/s);
        else if type==2
                W(i,i)=Ramsay_weight(e(i)/s);
            else if type==3
                    W(i,i)=Andrews_weight(e(i)/s);
                else if type==4
                        W(i,i)=Hampel_weight(e(i)/s);
                    else 
                        error('The value specified by the last parameter is wrong')
                    end
                end
            end
        end
    end
    count_iter++;
    Residuals(:,count_iter)=e;
    Weights(:,count_iter)=diag(W);
    if(norm(diag(W-W_previous))<couvergence_value)
        iterate=
end