作者shyfang (老師你有空？)
看板NCTU-STAT94G
標題Matlab stepwisefit
時間Wed Feb 25 21:46:49 2009
內建的好像有點問題，改寫如下
(目前我寫的是沒看到問題，有興趣可以試試看)
後面有testing data
pin=0.06
pout=0.05
keep=[]
intercept=true
會選入 x4 x1 x2 -x4 。 即，最後剩下 x1 x2
在STATS.inmodel中 row 代表第幾步，column代表第幾個變數，0:沒有選到 1:選到
% 沒有內設的參數，所以，每一個arguments都要填滿。
%
% STATS=stepwf(X,Y,pin,pout,keep,intercpt)
%
% Algorithm
% 1. k=1,加入最小RSS的變數 (至少有一個變數比 null model好)
% 2. k=k+1
% 3. 在(p-k+1)個變數中，找有最小的RSS的變數
% 4. 若其RSS的p-value<pin 加入。
% 否則，goto 8
% 5. 只要檢查前面k-1個選入的變數 (注意keep住的變數不可在其中)
% 找刪除後有最小partial F 的變數
%
% 6. 若刪除該變數後，其p-value>pout 則刪除該變數
% 7. 回到 2
% 8. STOP
%
% Arguments
% X:predictors
% Y:response
% pin:entering probability
% pout: removing probability
% keep: keep some predictors ALWAYS in model ('column numbers' or [])
% intercept: add intercept term or not (true or false)
%
% Return Values :
% STATS.inmodel=inmodel;
% STATS.SSE;
% STATS.dfE;
% STATS.SSR;
% STATS.dfR;
% STATS.sst;
% STATS.pf (partial F statistic)
%
function STATS=stepwf(X,Y,pin,pout,keep,intercept)
[nrow ncol]=size(X);
% first step
k=1;
sse=10^10;
SSE=zeros(ncol,1);
dfE=zeros(ncol,1);
SSR=zeros(ncol,1);
dfR=zeros(ncol,1);
pf=zeros(ncol,1);
index=setdiff(1:ncol,keep); %candidate variables
idx0=zeros(ncol,1);
if ~isempty(keep)
idx0(keep,1)=1;
end
for i=index
idx1=sort([find(idx0==1);i]);
if(intercept)
x=[ones(nrow,1),X(:,idx1)];
else
x=X(:,idx1);
end
beta=inv(x'*x)*x'*Y;
yhat=x*beta;
e=Y-yhat;
r=yhat-mean(Y);
ssr1=r'*r;
sse1=e'*e;
dfr=k;
dfe=nrow-k-intercept;
PF=ssr1/sse1/dfr*dfe;
pval=1-fcdf(PF,dfr,dfe);
% pval小表示和null model 有差異所以要選非null model
if sse1<sse
sse=sse1;
idx=idx1;
end
end
inmodel=zeros(1,ncol);
inmodel(idx)=1;
SSE(1)=sse;
dfR(1)=k+(~isempty(keep))*length(keep);
dfE(1)=nrow-intercept-dfR(1);
SSR(1)=ssr1;
pf(1)=PF;
sse0=sse;
t=Y-mean(Y);
sst=t'*t;
% setpwise
%先加一個再看看可不可以丟一個
pval0=0;
index=setdiff(1:ncol,keep);
while pval0<pin
[k ncol]=size(inmodel);
k=k+1;
sse0=SSE(k-1);
sse=10^10;
in0=sign(sum(inmodel,1));
for i=find(in0==0)
in=in0;
in(i)=1;
if(intercept)
x=[ones(nrow,1),X(:,in==1)];
else
x=X(:,in==1);
end
beta=inv(x'*x)*x'*Y;
yhat=x*beta;
e=Y-yhat; %residual
sse1=e'*e;
r=yhat-mean(Y); %regression
ssr1=r'*r;
if sse1<sse
flag=i;
sse=sse1;
ssr=ssr1;
end
end
dfr=sum(in);
dfe=nrow-dfr-intercept;
PF=(sse0-sse)/sse*dfe;
pval=1-fcdf(PF,1,dfe);
pval0=pval;
%若沒有增加變數，表示上一步就是最終結果，若有增加變數就要看前面有沒有可以
丟。
if pval<pin
inmodel=[inmodel;zeros(1,ncol)];
inmodel(k,flag)=1;
if ~isempty(keep)
inmodel(k,keep)=1;
end
SSE(k)=sse;dfE(k)=dfe;SSR(k)=ssr;dfR(k)=dfr;pf(k)=PF;
sse0=sse;
sse=0; %
% remove variable
clear flag;
out0=sign(sum(inmodel,1));
size(out0)
PF=10^10;
dfr=sum(out0)-1;
dfe=nrow-dfr-intercept;
indd=find(sign(sum(inmodel(1:(k-1),:),1))==1);%
%把keep的位置移掉
if ~isempty(keep)
indd(indd==keep)=[];
end
for i=indd
i
out=zeros(1,ncol);
out(i)=1
out=out0-out
if(intercept)
x=[ones(nrow,1),X(:,out==1)]; %intercept=true
else
x=X(:,out==1); %intercept=false
end
beta=inv(x'*x)*x'*Y;
yhat=x*beta;
e=Y-yhat;
sse1=e'*e;
r=yhat-mean(Y);
ssr1=r'*r;
PF1=(sse1-sse0)/sse0*(dfe-1);
if PF1<PF %挑一個最不顯著改善殘值的變數，看是否可刪除該變數
flag=i;
sse=sse1;
ssr=ssr1;
PF=PF1;
end
end
% disp('dfe');dfe;
pval=1-fcdf(PF,1,dfe-1);
%pval小表示二個model有差異，故不可刪除該變數
if pval>pout
[a b]=find(inmodel==1);
inmodel=[inmodel;zeros(1,ncol)];
k=k+1;
disp('k=');k
inmodel(k,flag)=-1;
SSE(k)=sse;dfE(k)=dfe;SSR(k)=ssr;dfR(k)=dfr;pf(k)=PF;
end
else
k=k-1;
break;
end
end
inmod=sign(cumsum(STATS.inmodel,1));
STATS.inmodel=inmod;
STATS.SSE=SSE(1:k);
STATS.dfE=dfE(1:k);
STATS.SSR=SSR(1:k);
STATS.dfR=dfR(1:k);
STATS.pf=pf(1:k);
STATS.sst=sst;
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.114.204.50
※ 編輯: shyfang 來自: 140.114.204.50 (02/25 21:49)
※ 編輯: shyfang 來自: 140.114.36.148 (03/04 11:03)
※ 編輯: shyfang 來自: 140.114.36.148 (04/12 13:52)