前幾天，偶然在網上看到，建立多元線性迴歸時須對其中的異常資料進行篩選剔除，那這裡的異常點指的是什麼呢？這裡的異常點指的是人為採集資料誤差或者某些異常的個例等等一些不太準確的資料。

例子說明

例如，對於下表中的一些資料我們來進行資料篩選。
在醫學上，糖尿病人的血糖量y(mmol/L)$y(mmol/L)$與總膽固醇 1 X (mmol / L)，甘油三脂X2(mmol/L)$X_2(mmol/L)$胰島素X3(μU/ml)$X_3(\mu U/ml)$糖化血紅蛋白X4($X_4($等有關。附表中是某醫院的原始資料。我們將對其中的異常資料進行篩選。

具體做法如下：
我們先通過Matlab$Matlab$對原始資料進行檢驗，對殘差進行分析，得到了殘差分析圖，剔除其中的異常點。

從圖上可以看出，第13$13$個點和第27$27$個點是異常點，這樣在做資料處理時就可以將其剔除。
下面是我遇到的一些問題：

• 至今沒有找到這個演算法的數學原理，沒看到書上有殘差向量的相關介紹。
• Matlab$Matlab$程式是改進網上的開源程式而來的，從程式中看不出他為什麼要這麼做。
• 這種篩選資料的方法是否只適應於線性迴歸，非線性迴歸是否也可以借鑑。

Matlab程式碼

clear all;
clc;

%輸入資料
z=xlsread('data.xls');
z1=z;
y=z(:,1);
X=[ones(size(y)),z(:,2:5)];
alpha=0.05;
b=inv((X'*X))*X'*y;  %迴歸係數
[n,ncolX] = size(X);

%剔除異常資料
wasnan=(isnan(y)|any(isnan(X),2));
havenans=any(wasnan);
[Q,R,perm]=qr(X,0);
p=ncolX;
RI=R\eye(p);
nu=max(0,n-p);
yhat=X*b;
r=y-yhat;  %殘差
normr=norm(r);
rmse=normr/sqrt(nu);    % Root mean square error.
tval=tinv((1-alpha/2),nu);
s2=rmse^2;
hatdiag=sum(abs(Q).^2,2);
ok=((1-hatdiag)>sqrt(eps(class(hatdiag))));
denom=(nu-1).*(1-hatdiag);
sigmai=zeros(length(denom),1);
sigmai(ok)=sqrt(max(0,(nu*s2/(nu-1))-(r(ok).^2./denom(ok))));
ser=sqrt(1-hatdiag).*sigmai;
rint=[(r-tval*ser) (r+tval*ser)];
rcoplot(r,rint);
kk=[];
for i=1:n
    if (rint(i,1)>0&&rint(i,2)>0)||(rint(i,1)<0&&rint(i,2)<0)
       kk=[kk,i]; 
    end
end
X(kk,:)=[];
y(kk)=[];

%剔除異常點後，求解迴歸係數
beta=inv((X'*X))*X'*y;  %迴歸係數

%迴歸分析
X1=X;
X1(:,1)=[];
n=size(y,1);    %觀察單位數
m=size(X,2);    
p1=m-1;          %自變數個數

alpha=0.05;  
yhat=X*beta;

%方差分析表（F檢驗）
SSR=(yhat-mean(y))'*(yhat-mean(y));  %迴歸平方和

SSE=(yhat-y)'*(yhat-y);        %殘差平方和

SST=(y-mean(y))'*(y-mean(y));    %總平方和

Fb=(SSR/(m-1))/(SSE/(n-m));       %顯著性檢驗的統計量

Falpha=2*(1-fcdf(abs(Fb),m-1,n-m));

table1=cell(4,6);    %建立元胞
table1(1,:)={'模型','偏差平方和','自由度','均方','F值','F.Sig'};
table1(2,1:6)={'迴歸',SSR,m-1,SSR/(m-1),Fb,Falpha};
table1(3,1:6)={'殘差',SSE,n-m,SSE/(n-m),[],[]};
table1(4,1:3)={'總和',SST,n-1};

C=diag(inv((X'*X)));
bj2=beta.*beta;     %迴歸係數平方
SSj=bj2(2:end)./C(2:end);     %偏回歸係數平方和

%決定係數檢驗
R2=SSR/SST;   %決定係數
R=sqrt(R2);   %複相關係數
RC=1-(1-R2)*(n-1)/(n-1-p1);
Sy=sqrt(SSE/(n-m));   %剩餘標準差
table2=cell(2,5);    %建立元胞
table2(1,:)={'模型','R','R平方','校正決定係數','剩餘標準差'};
table2(2,1)={1};
table2(2,2)={R};
table2(2,3)={R2};
table2(2,4)={RC};
table2(2,5)={Sy};

%t檢驗
s=zeros(m,1);       %迴歸係數的標準誤差
for i=1:m
    s(i,1)=sqrt(C(i))*sqrt(SSE/(n-m));
end

mnX=mean(X1);
MNX=repmat(mnX,n,1);   %複製mnX到矩陣MNX中
Ljj=diag((X1-MNX)'*(X1-MNX));     %Ljj的對角線元素為(X1-MNX)'*(X1-MNX)
Pj=beta(2:end).*sqrt(Ljj/SST);   %標準偏回歸係數；

t=zeros(m,1);
for i=1:m
    t(i,1)=beta(i,1)/s(i,1);
end

p2=zeros(m,1);
for i=1:m
    p2(i,1)=2*(1-tcdf(abs(t(i,1)),n-m));
end

table3=cell(m+1,6);    %建立元胞
table3(1,:)={'模型','偏回歸係數','迴歸係數的標準誤差','標準偏回歸係數','t值','P值'};
table3(2,:)={'常量',beta(1,1),s(1,1),[],t(1,1),p2(1,1)};

for kk=1:m-1
    table3(kk+2,:)={['x',num2str(kk)],beta(kk+1,1),s(kk+1,1),Pj(kk,1),t(kk+1,1),p2(kk+1,1)};
end

disp('係數分析表');
disp(table3);
disp('模型彙總');
disp(table2);
disp('方差分析表');
disp(table1);