stepwise.m

逐步回归算法程序

   [Q,R]=qr(X(:,inmodel),0);
   b = R\(Q'*y);
   r = y - X(:,inmodel)*b;
   RI = R\eye(size(R));
   xtxid = (sum(RI'.*RI'))';
else
   r = y;
end
terms = length(inmodel);
tss = y'*y;
rss = y'*y-r'*r;
if n-terms-1 > 0
   mse = r'*r./(n-terms-1);
elseif (rss < tss * sqrt(eps))
   mse = NaN;
else
   mse = Inf;
end
rmse = sqrt(mse);
if terms == 0
   prob = 1;
   r2   = 0;
   f    = 0;
elseif n-terms-1 <= 0
   f = NaN;
   prob = NaN;
   r2 = rss ./ tss;
   seb = Inf;
   t = Inf;
elseif mse == 0
   f = Inf;
   prob = 0;
   r2 = 1;
   seb = 0;
   t = Inf;
else
   f   = (rss./terms)./mse;
   prob = 1-fcdf(f,terms,(n-terms-1));
   r2  = rss./tss;
   seb = sqrt(xtxid)*rmse;
   t = b./seb;
end

beta = zeros(p,1);
tstats = zeros(p,1);
delta = zeros(p,1);
if ~isempty(inmodel)
   delta(inmodel) = seb*crit(terms);
   beta(inmodel) = b;
   tstats(inmodel) = t;
end


for index = notin
   Xnew = [X(:,inmodel) X(:,index)];
   [Qadd, Radd] = qr(Xnew,0);
   bnew = Radd\(Qadd'*y);
   if (n <= terms+2)
     rmse1 = Inf;
     seb = Inf;
   else
     rnew = y - Xnew*bnew;
     RInew = Radd\eye(terms+1,terms+1);
     xtxidnew = (sum(RInew'.*RInew'))';
     rmse1 = sqrt(sum(rnew.*rnew)./(n-terms-2));
     seb = sqrt(xtxidnew(terms+1))*rmse1;
   end
   
   beta(index) = bnew(terms+1);
   tstats(index) = bnew(terms+1)./seb;
   delta(index) = seb*crit(terms+1);
end
delta(isnan(delta)) = 0;
stats = [rmse r2 f prob];

function u = MakeErrorbarFigure(figure_handles,p,instance,screenht)
% Local function for creating Errorbar Plot Figure and its Uicontrols.

errorbar_fig = figure_handles(1);
stepht = min(50*(p+1),480);
set(errorbar_fig,'Units','Pixels','Tag','stepfig', ...
   'Resize','off', 'Name',['Stepwise Plot #',int2str(instance)], ...
   'Interruptible', 'off',...;
    'Position',[5+30*(instance-1) screenht-stepht-30*instance-50 300 stepht]);

figpos  = get(errorbar_fig,'Position');
step_axes = axes;
set(step_axes,'NextPlot','add','DrawMode','fast','Units','Pixels',  ...
   'Tag','stepaxes','Box','on','Ydir','reverse','Ylim',[0.5 p+0.5], ...
   'Ytick',(1:p),'Position',[30 70 figpos(3)-40 figpos(4)-90], ...
   'FontSize',9,'FontName','Geneva');  

parens = ['[],',int2str(errorbar_fig),',[',int2str(figure_handles(2)),' ',int2str(figure_handles(3)),'])']; 

u.outputpopup = uicontrol(errorbar_fig,'Style','Popup','String',...
   'Export|Parameters|Confidence Intervals|Terms In|Terms Out|All',...
   'Value',1, 'Units','Pixels','Position',[190 10 100 20],... 
   'CallBack',['stepwise(''output'',', parens], 'BackgroundColor','w');

u.scalebutton = uicontrol(errorbar_fig,'Style','radio','String',...
   'Scale Inputs','Units','Pixels','Position',[100 10 80 20],...
   'CallBack',['stepwise(''scale'',',parens]);

close_button = ...
   uicontrol('Style','Pushbutton','Units','Pixels','String','Close',... 
   'Position',[30 10 60 20],'Callback',['stepwise(''close'',',parens]);

xl = xlabel('Coefficients with Error Bars','Color','w');
yl = ylabel('Model Term Number');
set(xl,'FontSize',10,'FontName','Geneva');
set(yl,'FontSize',10,'FontName','Geneva');

function UpdateCoefficientsPlot(beta,delta,inmodel,notin,handles,bhandle, ...
                                errorbar_fig)
% Local function for plotting coefficients with error bars.
p = length(beta);
for k = 1:p
   set(handles(k),'XData',[beta(k)-delta(k); beta(k)+delta(k)]);
   set(bhandle(k),'XData',beta(k));
end
xmax = max(beta+delta);
xmin = min(beta-delta);
if (~(isfinite(xmin) & isfinite(xmax)))
   xx = beta + delta;
   xmax = max(xx(isfinite(xx)));
   if (isempty(xmax)), xmax = max(beta(isfinite(beta))); end
   if (isempty(xmax)), xmax = 1; end
   xx = beta - delta;
   xmin = min(xx(isfinite(xx)));
   if (isempty(xmin)), xmin = min(beta(isfinite(beta))); end
   if (isempty(xmin)), xmin = -1; end
end
xrange = xmax - xmin;
if (xrange == 0), xrange = 1; end

set(get(errorbar_fig,'CurrentAxes'), ...
    'XLim',[xmin-0.05*xrange xmax+0.05*xrange]);

notsig = find(beta-delta<0 & beta+delta>0);
sig = 1:p;
sig(notsig) = [];
set(handles(inmodel),'Color','g');
set(bhandle(inmodel),'Color','g');

set(handles(notin),'Color','r');
set(bhandle(notin),'Color','r');

set(handles(sig),'LineStyle','-');
set(handles(notsig),'LineStyle',':');

function MakeTableFigure(figure_handles,instance,p,screenht)
% Local function for creating Parameters Table Figure and its Uicontrols.

table_fig = figure_handles(2);
fight = p*12 + 100;
figure(table_fig);
set(table_fig,'Units','Pixels', 'Color','k', 'Tag','datafig',...
    'Name',['Stepwise Table #',int2str(instance)],'Resize','off', ...
    'Interruptible', 'off',...
    'Position',[290+30*instance screenht-(fight+30*instance)-50 340 fight]);
ax = axes;
set(ax,'Units','Pixels','Position',[10 30 380 fight-30],'Color','k',...
   'Xlim',[0 340],'Ylim',[0 fight-25],'Visible','off');

parens = ['[],',int2str(figure_handles(1)),',[',int2str(table_fig),' ',int2str(figure_handles(3)),'])'];
uicontrol('Style','Pushbutton','Units','Pixels','String','Close', ...
   'Position',[10 5 100 20],'Callback',['stepwise(''close'',',parens]);

uicontrol('Style','Pushbutton','Units','Pixels','Position',[230 5 100 20],...
   'Callback',['stepwise(''help'',',parens],'String','Help');

function WriteDataTable(beta,delta,p,figure_handles,inmodel,stats)
% Local function for writing out all the entries in the Parameter Data Table.

figure(figure_handles(2));
fight = p*12 + 100;
h=text(180,fight-32,'Confidence Intervals');
set(h,'FontName','Geneva','FontSize',9,'Color','w');
colheads = str2mat('Column #','Parameter','Lower ','Upper');
colheads1 = str2mat('RMSE','R-square','F','P');
z = [10 80 175 245];
z1 = [10 80 195 265];
h = zeros(4,1);
l = zeros(4,1);
for k = 1:4
   h(k)=text(z(k),fight-45,colheads(k,:));
   set(h(k),'FontName','Geneva','FontSize',9,'Color','w');
   l(k)=text(z1(k),20,colheads1(k,:));
   set(l(k),'FontName','Geneva','FontSize',9,'Color','w');
end

x = fight-50;
h = zeros(4,p+1);
for params = 1:p
   z = [30 120 200 270];
   x = x - 12;
   row = [params beta(params) beta(params)-delta(params) beta(params)+delta(params)];
   parens = [int2str(params),',',int2str(figure_handles(1)),',[',int2str(figure_handles(2)),' ',int2str(figure_handles(3)),'])'];
   bdf = ['stepwise(''click'',',parens];
   for k = 1:4
      if k == 1,
         h(k,params) = text(z(k),x,sprintf('%4i',row(k)));
         set(h(k,params),'HorizontalAlignment','right',...
            'FontName','Geneva','FontSize',9,'Tag','text',...
            'ButtonDownFcn',bdf,'Interruptible','off');
      else
         h(k,params) = text(z(k),x,sprintf('%11.4g',row(k)));
         set(h(k,params),'HorizontalAlignment','right','Tag','text',...
            'FontName','Geneva','FontSize',9,'ButtonDownFcn',bdf,...
            'Interruptible','off');
      end
      if ~isempty(inmodel) & any(inmodel == params)
         set(h(k,params),'Color','g');
      else
         set(h(k,params),'Color','r');
      end
   end
end
z = [35 120 205 290];
for k = 1:4
   h(k,p+1) = text(z(k),6,sprintf('%11.4g',stats(k)));
   set(h(k,p+1),'HorizontalAlignment','right','Tag','text',...
      'FontName','Geneva','FontSize',9,'Color','w');
end

function MakeHistoryFigure(inmodel,n,p,instance,stats,figure_handles,...
                           histud,num_mod)
% Local function for creating Stepwise History Figure.
screen = get(0,'ScreenSize');
screenht = screen(4);
screenwidth  = screen(3);

hist_fig = figure_handles(3);
figpos = get(hist_fig,'Position');
fight = p*12 + 100;
parens = ['[]',',',int2str(figure_handles(1)),',[',int2str(figure_handles(2)),' ',int2str(hist_fig),'])'];
if screenwidth > 970
   set(hist_fig,'WindowButtonMotionFcn',['stepwise(''histmotion'',',parens],...
      'Units','Pixels','Position',[670 screenht-310 300 230]);
else 
   set(hist_fig,'WindowButtonMotionFcn',['stepwise(''histmotion'',',parens],...
      'Units','Pixels','Position',[290+30*instance ...
         max(0,screenht-(fight+30*instance)-280) 300 230]); 
end
set(hist_fig,'Tag','histfig','Interruptible', 'off', ...
    'Name',['Stepwise History #',int2str(instance)]);

trm  = n-1:-1:n-p-1;
trml = length(trm);
low  = 0.025;
hi   = 0.975;

histud.chi2crit = chi2inv([low(ones(trml,1),1) hi(ones(trml,1),1)],[trm' trm']); 
figure(figure_handles(3));

UpdateHistoryPlot(stats,inmodel,n,p,histud,figure_handles,num_mod);
xlabel('Model Number');
ylabel('RMSE');

function UpdateHistoryPlot(stats,inmodel,n,p,histud,figure_handles,num_mod)
% Local function for adding new rmse error bar to History Plot.

figure(figure_handles(3))
rmse = stats(1);
mse = rmse*rmse;
terms = length(inmodel);
rmseci  =  sqrt([(mse*sqrt((n-terms-1)./histud.chi2crit(terms+1,2))) ...
      (mse*sqrt((n-terms-1)./histud.chi2crit(terms+1,1)))]);

cirange = diff(rmseci);
ylhi = rmseci(2) + 0.1*cirange;
yllo = rmseci(1) - 0.1*cirange;

tmp = NaN;
tmp = tmp(ones(p,1),1);
tmp(1:terms) = inmodel;

if ~isfield(histud, 'instore')
   histud.instore = tmp;
else  
   in = histud.instore;
   in = [in tmp];
   histud.instore = in;
   axes(get(figure_handles(3),'CurrentAxes'));
end
histud.rmsehandle(num_mod) = plot([num_mod num_mod],rmseci);

yl = get(get(figure_handles(3),'CurrentAxes'), 'Ylim');
yhi = max(ylhi,yl(2));
ylo = min(yllo,yl(1));
set(get(figure_handles(3),'CurrentAxes'),'Xlim',[0.5 num_mod+0.5],...
   'Xtick',(1:num_mod),'Ylim',[ylo yhi], ...
   'Position',[0.15 0.15 0.8 0.8],'NextPlot','add');
axes(get(figure_handles(3),'CurrentAxes'));
rmsedothndl = plot(num_mod,rmse,'.');

parens = [int2str(num_mod),',',int2str(figure_handles(1)),',[', ...
      int2str(figure_handles(2)),' ',int2str(figure_handles(3)),'])'];
bdf = ['stepwise(''history'',',parens];
set(histud.rmsehandle(num_mod),'ButtonDownFcn',bdf);
set(rmsedothndl,'MarkerSize',20,'ButtonDownFcn',bdf);
set(figure_handles(3),'UserData',histud);

function dohelp
str{1,1} = 'Stepwise Regression';
str{1,2} = { 'Stepwise Regression Interactive Graphic User Interface'
   ' ' 
   'The interface consists of three interactively linked figure windows:'
   '   * The Stepwise Regression Plot'
   '   * The Stepwise Regression Diagnostics Table'
   '   * The Stepwise History'
   ' '   
   'All three windows have hot regions. When your mouse is above' 
   'one of these regions, the pointer changes from an arrow to a circle.'
   'Mousing down at this point initiates some activity in the interface.'
   ' '   
   'To understand how to use a given window choose one of the items'
   'from the popup menu labelled Topics above.'
};

str{2,1} = 'Stepwise Plot';
str{2,2} = {
   'Stepwise Regression Plot'
   ' '
   'Coefficients and Confidence Intervals'
   ' '
   'This plot shows the regression coefficient and confidence interval for every'
   'term (in or out of the model.) The green lines represent terms in the model'
   'while red lines indicate that the term is not currently in the model.'
   ' '
   'Statistically significant terms are solid lines. Dotted lines show that '
   'the fitted coefficient is not significantly different from zero.'
   ' '
   'Clicking on a line in this plot toggles its state. That is, a term in the'
   'model (green line) gets removed (turns red), and terms out of the model' 
   '(red line) enter the model (turn green).'
   ' '
   'The coefficient for a term out of the model is the coefficient resulting from'
   'adding that term to the current model.'
   ' '
   'Scale Inputs'
   'Pressing the button labeled Scale Inputs centers and normalizes the columns'
   'of the input matrix to have a standard deviation of one.'
   ' '
   'Export'
   'This popup menu allows you to export variables from the stepwise function to'
   'the base workspace.'
   ' '
   'Close'
   'The Close button removes all the figure windows.'
};

str{3,1} = 'Stepwise Diagnostics';
str{3,2} = {
   'Stepwise Regression Diagnostics Table'
   ' '
   'Coefficients and Confidence Intervals'
   'The table at the top of the figure shows the regression coefficient and ' 
   'confidence interval for every term (in or out of the model.) The green rows' 
   'in the table represent terms in the model while red rows indicate terms not' 
   'currently in the model.'
   ' '
   'Clicking on a row in this table toggles the state of the corresponding term.' 
   'That is, a term in the model (green row) gets removed (turns red), and terms' 
   'out of the model (red rows) enter the model (turn green).'
   ' '
   'The coefficient for a term out of the model is the coefficient resulting from'
   'adding that term to the current model.'
   ' '
   'Additional Diagnostic Statistics'
   'There are also several diagnostic statistics at the bottom of the table:' 
   '   * RMSE - the root mean squared error of the current model.'
   '   * R-square - the amount of response variability explained by the model.'
   '   * F - the overall F statistic for the regression.'
   '   * P - the associated significance probability.'
   ' ' 
   'Close Button'
   'Shuts down all windows.'
   ' '
   'Help Button'
   'Activates on-line help.'
};

str{4,1} = 'Stepwise History';
str{4,2} = {
   'Stepwise History'
   ' '
   'Description'
   'This plot shows the RMSE and a confidence interval for every model generated' 
   'in the course of the interactive use of the other windows. '
   ' '
   'Recreating a Previous Model'
   'Clicking on one of these lines will re-create the current model at that point'
   'in the analysis using a NEW set of windows. You can thus compare the two' 
   'candidate models directly.'
};

for k = 1:4,
   str{k,2} = char(str{k,2});
end
helpwin(str,'Stepwise Regression','Stepwise Regression Help')