stepwise.m

逐步回归算法程序

function stepwise(X,y,inmodel,alpha);
%STEPWISE Interactive tool for stepwise regression.
%   STEPWISE(X,Y,INMODEL,ALPHA) fits a regression model of Y on the
%   columns of X specified in the vector INMODEL. ALPHA is the significance  
%   for testing each term in the model.  The default INMODEL selects all
%   columns of X.  The default ALPHA is 1 - (1 - 0.025).^(1/p), where p
%   is the number of columns in X.  This translates to plotted 95%
%   simultaneous confidence intervals (Bonferroni) for all the coefficients.
%
%   The least squares coefficient is plotted with a green filled circle.
%   A coefficient is not significantly different from zero if its confidence
%   interval crosses the white zero line. Significant model terms are plotted 
%   using solid lines. Terms not significantly different from zero are plotted
%   with dotted lines.
%
%   Click on the confidence interval lines to toggle the state of the model
%   coefficients. If the confidence interval line is green the term is in the
%   model. If the confidence interval line is red the term is not in the
%   model. Use the pop-up menu, Export, to move variables to the base
%   workspace. 
%
%   The fitted model always includes a constant term, but the displays
%   do not show the estimated value of that term.  You can compute the
%   constant using the formula
%
%         mean(y) - mean(X(:,IN))*P(IN),
%
%   where P is the vector of parameter (coefficient) estimates and IN
%   specifies the X columns that are in the model.
%
%   Reference: Draper, Norman and Smith, Harry, Applied Regression Analysis,
%   Second Edition, John Wiley & Sons, Inc. 1981 pp. 307-312.

%   Z. You 10-29-98,  B. Jones 6-14-94
%   Copyright 1993-2002 The MathWorks, Inc. 
%   $Revision: 2.30 $  $Date: 2002/01/17 21:32:01 $
if ~isstr(X) 
   action = 'start';
else
   action = X;
end

%On recursive calls get all necessary handles and data.
num_mod = 1;
if ~strcmp(action,'start')
   flag = y;
   errorbar_fig = inmodel;
   table_fig = alpha(1);
   hist_fig = alpha(2);
   
   ch = get(0,'Children');
   good = [find(ch==table_fig) find(ch==errorbar_fig) find(ch==hist_fig)]; 
   if length(good) ~= 3
%      warndlg('Sorry, STEPWISE requires all three figure windows to run.');
      disp('Closing all three STEPWISE figure windows.');
      close(errorbar_fig); close(hist_fig); close(table_fig);
      return;
   end
   
   figure_handles = [errorbar_fig table_fig hist_fig];   
   histud = get(hist_fig,'UserData');
   if (isfield(histud, 'curmod'))
      num_mod = histud.curmod;
      if (isequal(action, 'click'))
         num_mod = num_mod + 1;
         histud.curmod = num_mod;
         set(hist_fig,'UserData',histud);
      end
   end

   u = get(errorbar_fig,'Userdata');
   
   X           = u.X;
   y           = u.y;
   inmodel     = u.inmodel;
   handles     = u.handles;
   bhandle     = u.bhandle;
   crit        = u.crit;
   outputpopup = u.outputpopup;
   scalebutton = u.scalebutton;
   instance    = u.instance;
   texthandle  = u.texthandle;
   [n, p] = size(X);
end

switch action
   
case 'start'
   [n, p]=size(X);
   colidx = 1:p;
   
   % Argument Checks
   [n1,collhs] = size(y);
   if (n1==1)
      y = y(:);
      n1 = collhs;
   end
   if n~=n1, 
      error('The number of rows in Y must equal the number of rows in X.'); 
   end 
   
   wasnan = (isnan(y) | any(isnan(X),2));
   if (any(wasnan))
      y(wasnan) = [];
      X(wasnan,:) = [];
      n = length(y);
   end

   if nargin < 3
      inmodel = colidx;
   end
   
   if nargin == 4
      sigprob = 1 - alpha/2;
   else
      sigprob = (0.975).^(1/p);
   end
   
   % Monitor Size information.
   screen = get(0,'ScreenSize');
   screenht = screen(4);
   screenwidth  = screen(3);
   
   % Do initial statistical calculations.
   mx = mean(X);
   X = (X - mx(ones(n,1),:));
   y = y - mean(y);
   sx = std(X);
   scaledx = X./sx(ones(n,1),:);
   
   df = n-2:-1:n-p-1;
   crit = tinv(sigprob,df);
   crit(df < 1) = Inf;
   [beta, delta, tstats, stats] = FitModel(X,y,inmodel,crit);
   
   % Create the three windows for the Stepwise GUI
   table_fig = figure('Visible','off');
   hist_fig = figure('Visible','off');
   errorbar_fig = figure('Visible','off');
   whitebg(hist_fig,[0 0 0]);
   whitebg(errorbar_fig,[0 0 0]);
   %instance = length(findobj(0,'Tag','stepfig'))+1;
   instance = 1;
   %drawnow;
   figure_handles = [errorbar_fig table_fig hist_fig];
   
   % Set Up for Stepwise Plot Figure
   u = MakeErrorbarFigure(figure_handles,p,instance,screenht);
   u.texthandle = [];
   u.X          = X;
   u.y          = y;
   u.inmodel    = inmodel;
   u.crit       = crit;
   u.scaledx    = scaledx;
   
   yref = get(gca,'YLim');    
   plot([0 0],yref,'w-');
   
   % Make coefficients plot with error bars
   handles = plot([(beta-delta)'; (beta+delta)'],[colidx; colidx]);
   for k = 1:length(colidx),
      bhandle(k,1) = plot(beta(k),colidx(k),'.');
      parens = [int2str(k),',',int2str(errorbar_fig),',[',int2str(table_fig),' ',int2str(hist_fig),'])'];
      set(bhandle(k,1),'ButtonDownFcn',['stepwise(''click'',',parens],...
                       'Interruptible','off','MarkerSize',20)
      set(handles(k),'ButtonDownFcn',['stepwise(''click'',',parens],...
                     'Interruptible','off','LineWidth',2)
   end
   u.handles   = handles;
   u.bhandle   = bhandle;
   u.instance  = instance;
   notin = colidx;
   notin(inmodel) = [];
   UpdateCoefficientsPlot(beta,delta,inmodel,notin,handles,bhandle, ...
                          errorbar_fig)
   
   % Set Up for Stepwise Parameter Data Table.
   MakeTableFigure(figure_handles,instance,p,screenht)
   WriteDataTable(beta,delta,p,figure_handles,inmodel,stats) 
   
   % Set Up for History Figure
   histud.curmod = 1;
   MakeHistoryFigure(inmodel,n,p,instance,stats,figure_handles,histud,num_mod);
   
   parens = ['[]',',',int2str(errorbar_fig),',[',int2str(table_fig),' ',int2str(hist_fig),'])']; 
   set(errorbar_fig,'WindowButtonMotionFcn',['stepwise(''stepmotion'',',parens])  
   set(table_fig,'WindowButtonMotionFcn',['stepwise(''datamotion'',',parens])
   
   set(hist_fig,'Visible','on','HandleVisibility','callback');
   set(table_fig,'Visible','on','HandleVisibility','callback');
   set(errorbar_fig,'Userdata',u,'Visible','on',...
                    'HandleVisibility','callback');

   set(errorbar_fig,'DeleteFcn',['stepwise(''errorkill'',',parens])
   set(table_fig,'DeleteFcn',['stepwise(''tablekill'',',parens])
   set(hist_fig,'DeleteFcn',['stepwise(''histkill'',',parens])
   % End of initialization Phase
   
   
   % Callback for motion inside the Coefficients Table figure.
 case {'datamotion', 'stepmotion', 'histmotion'}
   ax = get(gcbf,'CurrentAxes');
   cp = get(ax,'CurrentPoint');
   cx = cp(1,1);
   cy = cp(1,2);
   if strcmp(action,'datamotion')
     online = (80 < cx) & (130 > cx) & (34 < cy) & (34+p*12 > cy);
   elseif strcmp(action,'stepmotion')
     if cy > 0.9 & cy < length(handles)+0.1
       xint = get(handles(round(cy)),'XData');
     else 
       xint = zeros(1,2);
     end
     online = (abs(cy-round(cy)) < 0.1) & (cx > xint(1)) & (cx < xint(2));
   else
     online = (abs(cx-round(cx)) < 0.02) & (round(cx) > 1);
   end
   togglecursor(online,gcbf);
   
 case {'errorkill', 'tablekill', 'histkill'}
%   warndlg(['STEPWISE requires all three figure windows to' ...
%	   ' run. The other windows have also been closed.']);
      disp('Closing all three STEPWISE figure windows.');
   if ~strcmp(action,'errorkill')
     set(errorbar_fig,'DeleteFcn','');
     close(errorbar_fig); 
   end
   if ~strcmp(action,'tablekill')
     set(table_fig,'DeleteFcn','');
     close(table_fig); 
   end
   if ~strcmp(action,'histkill')
     set(hist_fig,'DeleteFcn','');
     close(hist_fig); 
   end
 case 'output',
   bmf = get(errorbar_fig,'WindowButtonMotionFcn');
   bdf = get(errorbar_fig,'WindowButtonDownFcn');
   set(errorbar_fig,'WindowButtonMotionFcn','');
   set(errorbar_fig,'WindowButtonDownFcn','');
   
   outval = get(outputpopup,'Value');
   switch outval
   case 1,
      return,         
   case 2,
      promptstr = 'Enter names for the parameters or use the supplied default.'; 
      entrystr = 'beta';
   case 3,
      promptstr = 'Enter names for the confidence intervals or use the supplied default.';
      entrystr = 'betaci';
   case 4,
      promptstr = 'Enter names for the terms in the model or use the supplied default.';
      entrystr = 'in';
   case 5,
      promptstr = 'Enter names for the terms out of the model or use the supplied default.';
      entrystr = 'out';
   case 6,
      promptstr = 'Enter names for statistical outputs or use the defaults.';
      entrystr = 'beta,betaci,in,out';
   end
   
   dlghandle = dialog('Visible','off','WindowStyle','modal');
   
   set(dlghandle,'Units','pixels','Position',[100 300 400 100],'Visible','on');
   uicontrol(dlghandle,'Style','text','Units','pixels',...
      'Position',[10 75 380 20],...
      'String',promptstr);
   
   u.texthandle = uicontrol(dlghandle,'Style','edit','Units','pixels',...
      'String',entrystr,'HorizontalAlignment','center',...
      'Position',[50 40 300 20], 'BackgroundColor','w');
   parens = ['[],',int2str(errorbar_fig),',[',int2str(table_fig),' ',int2str(hist_fig),'])'];
   ok_button = uicontrol('Style','Pushbutton','Units','pixels', ...
      'Position',[180 5 40 20], ...
      'Callback',['stepwise(''toworkspace'',' parens ';close(gcbf)'], ...
      'String','OK');  
   
   set(dlghandle,'HandleVisibility','off')
   set(errorbar_fig,'UserData',u);
   set(errorbar_fig,'WindowButtonMotionFcn',bmf);
   set(errorbar_fig,'WindowButtonDownFcn',bdf);
   
case 'toworkspace',
   if get(scalebutton,'Value') == 1
      X = u.scaledx;
   end
   
   [beta, delta, tstats, stats] = FitModel(X,y,inmodel,crit);
   notin = 1:p;
   notin(inmodel) = [];
   
   outval = get(outputpopup,'Value');
   
   s = get(texthandle,'String');
   
   switch outval
   case 1,
      return,         
   case 2,
      assignin('base',s,beta);
   case 3,
      assignin('base',s,[beta-delta beta+delta]);
   case 4,
      assignin('base',s,inmodel);
   case 5,
      assignin('base',s,notin);
   case 6,
      cpos=findstr(s,',');
      assignin('base',s(1:cpos(1)-1),beta);
      assignin('base',s((cpos(1)+1):(cpos(2)-1)),[beta-delta  beta+delta]);
      assignin('base',s((cpos(2)+1):(cpos(3)-1)),inmodel);
      assignin('base',s((cpos(3)+1):end),notin);
   end
   
   set(outputpopup,'Value',1);
   
case 'click'
   cf = gcbf;
   figure(errorbar_fig);
   if ~isempty(inmodel) & any(inmodel == flag)
      k = find(inmodel==flag);
      u.inmodel(k) = [];
   else
      u.inmodel = sort([inmodel(:); flag]);
   end
   inmodel = u.inmodel;
   if get(scalebutton,'Value') == 1
      X = u.scaledx;
   end
   
   [beta, delta, tstats, stats] = FitModel(X,y,inmodel,crit);
   notin = 1:p;
   notin(inmodel) = [];
   UpdateCoefficientsPlot(beta,delta,inmodel,notin,handles,bhandle, ...
                          errorbar_fig)
   set(errorbar_fig,'UserData',u);
   
   figure(table_fig);
   ch = findobj(table_fig,'Tag','text');
   delete(ch);
   WriteDataTable(beta,delta,p,figure_handles,inmodel,stats)
   
   figure(hist_fig);
   UpdateHistoryPlot(stats,inmodel,n,p,histud,figure_handles,num_mod)
   figure(cf);   
   
   % Callback for Scale Inputs radio button.
case 'scale',
   if get(scalebutton,'Value') == 1
      X = u.scaledx;
   end
   [beta, delta, tstats, stats] = FitModel(X,y,inmodel,crit);
   notin = 1:p;
   notin(inmodel) = [];
   UpdateCoefficientsPlot(beta,delta,inmodel,notin,handles,bhandle, ...
                          errorbar_fig)
   
   figure(table_fig);
   ch = findobj(table_fig,'Tag','text');
   delete(ch);
   WriteDataTable(beta,delta,p,figure_handles,inmodel,stats)
   
   figure(errorbar_fig);
   
case 'history',
   in = histud.instore;
   in = in(:,flag);
   in = in(find(~isnan(in)));
   set(hist_fig,'Pointer','watch');
   stepwise(X,y,in)
   set(hist_fig,'Pointer','arrow');
   
case 'close',
   set(table_fig,'DeleteFcn','')
   set(errorbar_fig,'DeleteFcn','')
   set(hist_fig,'DeleteFcn','')
   close(table_fig);
   close(errorbar_fig);
   close(hist_fig);
   
case 'help',
   dohelp;
end

function togglecursor(online, fighandle)
% toggle the cursor on and off the line

cursorstate = get(fighandle,'Pointer');
if online & strcmp(cursorstate,'arrow')
   set(fighandle,'Pointer','circle');
elseif ~online & strcmp(cursorstate,'circle'),
   set(fighandle,'Pointer','arrow');
end

function [beta, delta, tstats, stats] = FitModel(X,y,inmodel,crit);
% Local function for fitting stepwise regression model.
[n, p] = size(X);
colidx = 1:p;
notin  = colidx;
notin(inmodel) = [];

if ~isempty(inmodel)