作者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)