contrib/mul/mbl/mbl_stepwise_regression.cxx
Go to the documentation of this file.
00001 // This is mul/mbl/mbl_stepwise_regression.cxx
00002 #include "mbl_stepwise_regression.h"
00003 //:
00004 // \file
00005 // \brief Conduct stepwise regression
00006 // \author Martin Roberts
00007 
00008 #include <vcl_algorithm.h>
00009 #include <vcl_iterator.h>
00010 #include <vcl_utility.h>
00011 // not used? #include <vcl_cmath.h>
00012 #include <vcl_cassert.h>
00013 // not used? #include <vcl_vector.h>
00014 // not used? #include <vcl_iostream.h>
00015 #include <vnl/vnl_math.h>
00016 #include <vnl/algo/vnl_svd.h>
00017 #include <mbl/mbl_stl.h>
00018 
00019 
00020 mbl_stepwise_regression::mbl_stepwise_regression(const vnl_matrix<double>& x,
00021                                                  const vnl_vector<double>& y):
00022         x_(x),y_(y),
00023         num_examples_(y.size()),
00024         num_vars_(x.cols()),
00025         XtX_(num_vars_+1,num_vars_+1,0.0),XtY_(num_vars_+1,0.0),
00026         FthreshAdd_(2.07),FthreshRemove_(1.8),
00027         mode_(mbl_stepwise_regression::eFORWARDS)
00028 {
00029     assert(x.rows() == num_examples_);
00030     //Initialise the extended covariances
00031 
00032 #if 0 // The calculation is as follows
00033     do
00034     {
00035         // XtX += [x, -1]' * [x, -1]
00036         const vnl_vector<double> &x=inputs.current();
00037         double y = outputs[inputs.index()] ? 1.0 : -1.0;
00038         vnl_vector<double> xp(k+1);
00039         xp.update(x, 0);
00040         xp(k) = -1.0;
00041         XtX += outer_product(xp, xp);
00042         double y = outputs[inputs.index()] ? 1.0 : -1.0;
00043         XtY += y * xp;
00044     } while (inputs.next());
00045 
00046 #else// However the following version is faster
00047 
00048     rss_ = 0.0;
00049     for (unsigned irow=0;irow<num_examples_;++irow)
00050     {
00051         // XtX += [x, -1]' * [x, -1]
00052         const double* px=x_[irow];
00053         double yval = y_[irow];
00054         for (unsigned i=0; i<num_vars_; ++i)
00055         {
00056             //Just work out in triangular form first
00057             for (unsigned j=0; j<i; ++j)
00058             {
00059                 XtX_(i,j) += px[i] * px[j];
00060             }
00061             XtX_(i,i) += vnl_math_sqr(px[i]);
00062             XtX_(i,num_vars_) -= px[i];
00063             XtY_(i) += yval * px[i];
00064         }
00065         XtY_(num_vars_) += yval * -1.0;
00066         rss_ += yval * yval;
00067     }
00068     //Initialise residual sum of squares to total sum of y squares about mean
00069     double sumY = -XtY_(num_vars_);
00070     rss_ -= (sumY*sumY)/double(num_examples_);
00071     //Set other half by symmetry
00072     for (unsigned i=0; i<num_vars_; ++i)
00073     {
00074         for (unsigned j=0; j<i; ++j)
00075         {
00076             XtX_(j,i) += XtX_(i,j);
00077         }
00078         XtX_(num_vars_,i) = XtX_(i,num_vars_);
00079     }
00080     XtX_(num_vars_, num_vars_) = (double) num_examples_;
00081 
00082 #endif // !0
00083 }
00084 
00085 void mbl_stepwise_regression::operator()()
00086 {
00087     using mbl_stepwise_regression_helpers::lsfit_this_basis;
00088     if (mode_ == eFORWARDS)
00089     {
00090         basis_.clear(); //start from an empty basis
00091         basis_complement_.clear();
00092         //And fill the complement with 1..n
00093         mbl_stl_increments_n(vcl_inserter(basis_complement_,basis_complement_.end()),
00094                             num_vars_,0);
00095 
00096         //First fit a starting minimally sized basis best sum of squares first
00097         unsigned  min_basis = vcl_min(num_vars_/10,num_examples_/5);
00098         if (min_basis<1) min_basis = 1;
00099         bool forceAdd=true;
00100         for (unsigned i=0; i<min_basis;i++)
00101         {
00102             add_variable(forceAdd);
00103         }
00104         //And now carry on and do forward stepwise (with some removals)
00105         do_forward_stepwise_regression();
00106     }
00107     else // Backwards
00108     {
00109         // start from a full basis
00110         mbl_stl_increments_n(vcl_inserter(basis_,basis_.end()),
00111                             num_vars_,0);
00112         basis_complement_.clear();
00113         do_backward_stepwise_regression();
00114     }
00115     // And ensure solution is set up for the final basis
00116     lsfit_this_basis fitter(x_,y_,XtX_,XtY_);
00117     fitter.set_basis(basis_);
00118     fitter();
00119     weights_ = fitter.weights();
00120 }
00121 
00122 void mbl_stepwise_regression::do_forward_stepwise_regression()
00123 {
00124     bool carryOn=true;
00125     while (carryOn)
00126     {
00127         // Now try another addition step followed by one elimination
00128         bool addedOne = add_variable();
00129         bool removedOne = remove_variable();
00130         carryOn = addedOne || removedOne ;
00131     }
00132 }
00133 
00134 void mbl_stepwise_regression::do_backward_stepwise_regression()
00135 {
00136     bool carryOn=true;
00137     while (carryOn)
00138     {
00139         //try to remove all the unwanted variables
00140         carryOn = remove_variable();
00141     }
00142     //Having removed everything we can check if anything might now come back in that had been removed earlier
00143     do_forward_stepwise_regression();
00144 }
00145 
00146 bool mbl_stepwise_regression::add_variable(bool forceAdd)
00147 {
00148     //Loop over all variables not in the basis
00149     //Fit an extended basis with that variable in it
00150     //Store the best such basis (best residual sum of squares)
00151     //If significant improvement (or forceAdd set) then add the variable
00152 
00153     using mbl_stepwise_regression_helpers::lsfit_this_basis;
00154 
00155     if (basis_complement_.empty()) return false;
00156 
00157     lsfit_this_basis fitter(x_,y_,XtX_,XtY_);
00158     fitter.set_basis(basis_);
00159     double FratioMax=-1.0;
00160 
00161     int  knew= -1;
00162     vcl_set<unsigned>::const_iterator candIter=basis_complement_.begin();
00163     vcl_set<unsigned>::const_iterator candIterEnd=basis_complement_.end();
00164     double rssNew=rss_;
00165     while (candIter != candIterEnd)
00166     {
00167         unsigned k = *candIter;
00168         double rssPrime = fitter.add(k);
00169         double F = f_ratio(rssPrime,rss_,1);
00170         if (F>FratioMax)
00171         {
00172             FratioMax= F;
00173             knew = k;
00174             rssNew = rssPrime;
00175         }
00176         ++candIter;
00177     }
00178     bool really_added = false;
00179     bool significant = test_significance(rssNew,rss_,FthreshAdd_);
00180     if (significant || forceAdd)
00181     {
00182         if (knew>=0)
00183         {
00184             vcl_pair<vcl_set<unsigned>::iterator,bool> inserted = basis_.insert(knew);
00185             really_added = inserted.second;
00186             basis_complement_.erase(knew);
00187             rss_ = rssNew;
00188         }
00189     }
00190     return significant && really_added;
00191 }
00192 
00193 bool mbl_stepwise_regression::remove_variable()
00194 {
00195     //Loop over all variables in the basis
00196     //Fit a diminished basis with that variable in it
00197     //Store the least diminishing (to residual sum of squares) variable
00198     //If the reduction is not significant then eliminate this variable
00199 
00200     using mbl_stepwise_regression_helpers::lsfit_this_basis;
00201     lsfit_this_basis fitter(x_,y_,XtX_,XtY_);
00202     fitter.set_basis(basis_);
00203     double min_Fratio=1.0E30;
00204     double rssBest=rss_;
00205     int knew=-1;
00206     vcl_set<unsigned>::const_iterator candIter=basis_.begin();
00207     vcl_set<unsigned>::const_iterator candIterEnd=basis_.end();
00208     while (candIter != candIterEnd)
00209     {
00210         unsigned k = *candIter++;
00211         double rssNew = fitter.remove(k);
00212         double F = f_ratio(rss_,rssNew,1);
00213         if (F < min_Fratio)
00214         {
00215             min_Fratio = F;
00216             knew = k;
00217             rssBest = rssNew;
00218         }
00219     }
00220     bool significant = test_significance(rss_,rssBest,FthreshRemove_);
00221     if (!significant)
00222     {
00223         basis_.erase(knew);
00224         basis_complement_.insert(knew);
00225         rss_ = rssBest;
00226     }
00227     return !significant;
00228 }
00229 
00230 //---------------------------------------------------------------------------
00231 //------------------------------ Helpers ------------------------------------
00232 //---------------------------------------------------------------------------
00233 
00234 double mbl_stepwise_regression_helpers::lsfit_this_basis::add(unsigned k)
00235 {
00236     basis_.insert(k);
00237     double rss = (*this)();
00238     basis_.erase(k);
00239     return rss;
00240 }
00241 
00242 double mbl_stepwise_regression_helpers::lsfit_this_basis::remove(unsigned k)
00243 {
00244     basis_.erase(k);
00245     double rss = (*this)();
00246     basis_.insert(k);
00247     return rss;
00248 }
00249 
00250 double mbl_stepwise_regression_helpers::lsfit_this_basis::operator()()
00251 {
00252     // Find the solution to X w = Y;
00253     // However it is easier to find X' X w = X' Y;
00254     // because X is n_train x n_dims whereas X'X is n_dims x n_dims
00255 
00256     //Solve by SVD
00257     unsigned ndims = basis_.size();
00258     //Create working copies of mtrices containing just the subset of variables in the basis
00259     vnl_matrix<double> XtX(1+ndims,1+ndims);
00260 
00261     vcl_set<unsigned>::iterator basisVarIter=basis_.begin();
00262     vcl_set<unsigned>::iterator basisVarIterEnd=basis_.end();
00263     unsigned i = 0;
00264     vnl_vector<double> XtY(ndims+1, 0.0);
00265 
00266     while (basisVarIter != basisVarIterEnd)
00267     {
00268         unsigned k1=*basisVarIter++;
00269         vcl_set<unsigned>::iterator basisVarInnerIter=basis_.begin();
00270         unsigned j = 0;
00271         //Set half of off-diagonals
00272         while (*basisVarInnerIter < k1) //NB set is ordered
00273         {
00274             unsigned k2=*basisVarInnerIter++;
00275             XtX(i,j) = XtX_(k1,k2);
00276             ++j;
00277         }
00278         XtX(i,i) = XtX_(k1,k1); //diagonal
00279         XtX(i,ndims) = XtX_(k1,num_vars_); //extra bias column
00280         XtY(i) = XtY_(k1);
00281         ++i;
00282     }
00283     XtY(ndims) = XtY_(num_vars_);
00284     //Copy the other half by symmetry
00285     for (unsigned i=0;i<ndims;++i)
00286     {
00287         for (unsigned j=0;j<i;++j)
00288         {
00289             XtX(j,i) = XtX(i,j);
00290         }
00291         XtX(ndims,i) = XtX(i,ndims);
00292     }
00293 
00294     XtX(ndims,ndims) = double (num_examples_);
00295 
00296     vnl_svd<double> svd(XtX, 1.0e-12); // 1e-12 = zero-tolerance for singular values
00297     weights_ = svd.solve(XtY);
00298 
00299     double rss=0.0;
00300     //Now compute the residual sum of squares
00301     for (unsigned i=0;i<num_examples_;++i)
00302     {
00303         const double* pDataRow=x_[i];
00304         double ypred = 0.0;
00305         basisVarIter=basis_.begin();
00306         vnl_vector<double>::iterator weightIter = weights_.begin();
00307         while (basisVarIter != basisVarIterEnd)
00308         {
00309             ypred += pDataRow[*basisVarIter++] * (*weightIter++);
00310         }
00311         ypred -= *weightIter; //Final -1 term
00312         double dy = y_[i] - ypred;
00313         rss += dy*dy;
00314     }
00315     return rss;
00316 }