00001
00002 #include "mbl_stepwise_regression.h"
00003
00004
00005
00006
00007
00008 #include <vcl_algorithm.h>
00009 #include <vcl_iterator.h>
00010 #include <vcl_utility.h>
00011
00012 #include <vcl_cassert.h>
00013
00014
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
00031
00032 #if 0 // The calculation is as follows
00033 do
00034 {
00035
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
00052 const double* px=x_[irow];
00053 double yval = y_[irow];
00054 for (unsigned i=0; i<num_vars_; ++i)
00055 {
00056
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
00069 double sumY = -XtY_(num_vars_);
00070 rss_ -= (sumY*sumY)/double(num_examples_);
00071
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();
00091 basis_complement_.clear();
00092
00093 mbl_stl_increments_n(vcl_inserter(basis_complement_,basis_complement_.end()),
00094 num_vars_,0);
00095
00096
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
00105 do_forward_stepwise_regression();
00106 }
00107 else
00108 {
00109
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
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
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
00140 carryOn = remove_variable();
00141 }
00142
00143 do_forward_stepwise_regression();
00144 }
00145
00146 bool mbl_stepwise_regression::add_variable(bool forceAdd)
00147 {
00148
00149
00150
00151
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
00196
00197
00198
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
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
00253
00254
00255
00256
00257 unsigned ndims = basis_.size();
00258
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
00272 while (*basisVarInnerIter < k1)
00273 {
00274 unsigned k2=*basisVarInnerIter++;
00275 XtX(i,j) = XtX_(k1,k2);
00276 ++j;
00277 }
00278 XtX(i,i) = XtX_(k1,k1);
00279 XtX(i,ndims) = XtX_(k1,num_vars_);
00280 XtY(i) = XtY_(k1);
00281 ++i;
00282 }
00283 XtY(ndims) = XtY_(num_vars_);
00284
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);
00297 weights_ = svd.solve(XtY);
00298
00299 double rss=0.0;
00300
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;
00312 double dy = y_[i] - ypred;
00313 rss += dy*dy;
00314 }
00315 return rss;
00316 }