contrib/mul/mbl/mbl_sample_stats_1d.cxx
Go to the documentation of this file.
00001 #include "mbl_sample_stats_1d.h"
00002 //:
00003 // \file
00004 #include <vsl/vsl_vector_io.h>
00005 #include <vcl_cassert.h>
00006 #include <vcl_cmath.h>
00007 #include <vcl_limits.h>
00008 #include <vcl_algorithm.h>
00009 
00010 
00011 //=========================================================================
00012 mbl_sample_stats_1d::mbl_sample_stats_1d(const vcl_vector<double> &samples)
00013 {
00014   clear();
00015   for (unsigned i=0, n=samples.size(); i<n; ++i)
00016   {
00017     add_sample(samples[i]);
00018   }
00019 }
00020 
00021 
00022 //=========================================================================
00023 mbl_sample_stats_1d::mbl_sample_stats_1d(const vnl_vector<double> &samples)
00024 {
00025   clear();
00026   for (unsigned i=0, n=samples.size(); i<n; ++i)
00027   {
00028     add_sample(samples[i]);
00029   }
00030 }
00031 
00032 
00033 //=========================================================================
00034 mbl_sample_stats_1d::mbl_sample_stats_1d()
00035 {
00036   clear();
00037 }
00038 
00039 
00040 //=========================================================================
00041 mbl_sample_stats_1d::~mbl_sample_stats_1d()
00042 {
00043 }
00044 
00045 
00046 //=========================================================================
00047 void mbl_sample_stats_1d::clear()
00048 {
00049   samples_.resize(0);
00050   stats_1d_.clear();
00051   use_mvue_=true;
00052 }
00053 
00054 
00055 //=========================================================================
00056 void mbl_sample_stats_1d::add_sample(double v)
00057 {
00058   stats_1d_.obs(v);
00059   samples_.push_back(v);
00060   return;
00061 }
00062 
00063 
00064 //=========================================================================
00065 unsigned mbl_sample_stats_1d::n_samples() const
00066 {
00067   return samples_.size();
00068 }
00069 
00070 
00071 //=========================================================================
00072 double mbl_sample_stats_1d::mean() const
00073 {
00074   return stats_1d_.mean();
00075 }
00076 
00077 
00078 //=========================================================================
00079 double mbl_sample_stats_1d::mean_of_absolutes() const
00080 {
00081   double abs_sum = 0;
00082   for (unsigned i=0, n=samples_.size(); i<n; ++i)
00083     abs_sum+=vcl_fabs(samples_[i]);
00084   return abs_sum/samples_.size();
00085 }
00086 
00087 
00088 //=========================================================================
00089 double mbl_sample_stats_1d::median() const
00090 {
00091   double ret;
00092 
00093   if (samples_.size()>0)
00094   {
00095     if ( samples_.size() % 2 == 0 )
00096     {
00097       unsigned index = samples_.size() / 2 - 1;
00098 
00099       vcl_vector<double> tmp=samples_;
00100 
00101       vcl_vector<double>::iterator index_it0 = tmp.begin() + index;
00102       vcl_nth_element(tmp.begin(),index_it0,tmp.end(),vcl_less<double>());
00103       double v0 = *index_it0;
00104 
00105       vcl_vector<double>::iterator index_it1 = tmp.begin() + index + 1;
00106       vcl_nth_element(tmp.begin(),index_it1,tmp.end(),vcl_less<double>());
00107       double v1 = *index_it1;
00108       ret = v0 + v1;
00109       ret /= 2.0;
00110     }
00111     else
00112     {
00113       unsigned index = (samples_.size() - 1) / 2;
00114 
00115       vcl_vector<double> tmp=samples_;
00116 
00117       vcl_vector<double>::iterator index_it = tmp.begin() + index;
00118       vcl_nth_element(tmp.begin(),index_it,tmp.end(),vcl_less<double>());
00119 
00120       ret = *index_it;
00121     }
00122   }
00123   else // crazy value if  no samples
00124   {
00125     ret = vcl_numeric_limits<double>::max();
00126   }
00127   return ret;
00128 }
00129 
00130 
00131 //=========================================================================
00132 double mbl_sample_stats_1d::quantile(double q) const
00133 {
00134   const unsigned n = samples_.size();
00135 
00136   // These checks are only asserts because client code is responsible for avoiding these errors.
00137   assert(q>=0.0 && q<=1.0);
00138   assert(n>0);
00139 
00140   // Map the specified quantile to a real-valued "index", i.e. a float lying between 2 integer indices
00141   double float_index = (n-1)*q;
00142 
00143   // Get the integer index immediately below (and enforce the bounds)
00144   double f0 = vcl_floor(float_index);
00145   f0 = f0<0.0 ? 0.0 : f0>n-1.0 ? n-1.0 : f0;
00146   unsigned i0 = static_cast<unsigned>(f0);
00147 
00148   // Get the integer index immediately above (and enforce the bounds)
00149   double f1 = vcl_ceil(float_index);
00150   f1 = f1<0.0 ? 0.0 : f1>n-1.0 ? n-1.0 : f1;
00151   unsigned i1 = static_cast<unsigned>(f1);
00152 
00153   // Get the 2 values bracketing the specified quantile position
00154   vcl_vector<double> tmp = samples_;
00155 
00156   vcl_vector<double>::iterator index_it0 = tmp.begin() + i0;
00157   vcl_nth_element(tmp.begin(), index_it0, tmp.end(), vcl_less<double>());
00158   double v0 = *index_it0;
00159 
00160   vcl_vector<double>::iterator index_it1 = tmp.begin() + i1;
00161   vcl_nth_element(tmp.begin(), index_it1, tmp.end(), vcl_less<double>());
00162   double v1 = *index_it1;
00163 
00164   // Linearly interpolate between the 2 values
00165   double f = float_index - f0;
00166   double ret = ((1.0-f)*v0) + (f*v1);
00167   return ret;
00168 }
00169 
00170 
00171 //=========================================================================
00172 double mbl_sample_stats_1d::nth_percentile(int n) const
00173 {
00174   if (samples_.size()==0)
00175     return vcl_numeric_limits<double>::max();
00176 
00177   double fact = double(n)/100.0;
00178   int index=int(fact*(samples_.size()-1));
00179   vcl_vector<double> tmp=samples_;
00180 
00181   vcl_vector<double>::iterator index_it = tmp.begin() + index;
00182   vcl_nth_element(tmp.begin(),index_it,tmp.end(),vcl_less<double>());
00183   double ret = *index_it;
00184   return ret;
00185 }
00186 
00187 
00188 //=========================================================================
00189 double mbl_sample_stats_1d::variance() const
00190 {
00191   double v=0;
00192 
00193   if (samples_.size()>1)
00194   {
00195     double mean_v = mean();
00196     double sum_sq = sum_squares();
00197     v = sum_sq - samples_.size()*(mean_v * mean_v);
00198 
00199     if (use_mvue_)
00200     {
00201       v /= (samples_.size()-1);
00202     }
00203     else
00204     {
00205       v /= samples_.size();
00206     }
00207   }
00208 
00209   return v;
00210 }
00211 
00212 
00213 //=========================================================================
00214 double mbl_sample_stats_1d::sd() const
00215 {
00216   return vcl_sqrt(variance());
00217 }
00218 
00219 
00220 //=========================================================================
00221 double mbl_sample_stats_1d::stdError() const
00222 {
00223   double se = variance();
00224   if (use_mvue_)
00225   {
00226     se /= samples_.size()-1;
00227   }
00228   else
00229   {
00230     se /= samples_.size();
00231   }
00232 
00233   return vcl_sqrt(se);
00234 }
00235 
00236 
00237 //=========================================================================
00238 double mbl_sample_stats_1d::skewness() const
00239 {
00240   double skew = 0;
00241 
00242   // skew
00243   // calculated as
00244   // ( Sum_i (Y_i-MEAN)^3 ) / ((N-1)*sigma^3)
00245   // where N is the number of samples
00246   // sigma is the standard deviation
00247 
00248   if (samples_.size()>1)
00249   {
00250     double s=sd();
00251     double m=mean();
00252 
00253     for (unsigned i=0, n=samples_.size(); i<n; ++i)
00254     {
00255       double tmp=samples_[i]-m;
00256       skew += (tmp*tmp*tmp) ;
00257     }
00258 
00259     skew /= ( (samples_.size()-1) * s * s * s );
00260   }
00261 
00262   return skew;
00263 }
00264 
00265 
00266 //=========================================================================
00267 double mbl_sample_stats_1d::kurtosis() const
00268 {
00269   double kurt = 0;
00270 
00271   // kurtosis
00272   // calculated as
00273   // -3 + ( Sum_i (Y_i-MEAN)^4 ) / ((N-1)*sigma^4)
00274   // where N is the number of samples
00275   // sigma is the standard deviation
00276 
00277   if (samples_.size()>1)
00278   {
00279     double s=sd();
00280     double m=mean();
00281 
00282     for (unsigned i=0, n=samples_.size(); i<n; ++i)
00283     {
00284       double tmp=samples_[i]-m;
00285       kurt += (tmp*tmp*tmp*tmp) ;
00286     }
00287 
00288     kurt /= ( (samples_.size()-1) * s * s * s *s);
00289     kurt -= 3;
00290   }
00291   return kurt;
00292 }
00293 
00294 
00295 //=========================================================================
00296 double mbl_sample_stats_1d::min() const
00297 {
00298   if (samples_.size()==0) return vcl_numeric_limits<double>::max();
00299   else return stats_1d_.min();
00300 }
00301 
00302 
00303 //=========================================================================
00304 double mbl_sample_stats_1d::max() const
00305 {
00306   if (samples_.size()==0) return vcl_numeric_limits<double>::min();
00307   else return stats_1d_.max();
00308 }
00309 
00310 
00311 //=========================================================================
00312 double mbl_sample_stats_1d::sum() const
00313 {
00314   return stats_1d_.sum();
00315 }
00316 
00317 
00318 //=========================================================================
00319 double mbl_sample_stats_1d::sum_squares() const
00320 {
00321   return stats_1d_.sumSq();
00322 }
00323 
00324 
00325 //=========================================================================
00326 double mbl_sample_stats_1d::rms() const
00327 {
00328   double ms=sum_squares()/stats_1d_.nObs();
00329   return vcl_sqrt( ms );
00330 }
00331 
00332 
00333 //=========================================================================
00334 mbl_sample_stats_1d& mbl_sample_stats_1d::operator+=(const mbl_sample_stats_1d& s1)
00335 {
00336   // add new samples
00337   for (unsigned i=0;i<s1.samples().size();++i)
00338   {
00339     add_sample(s1.samples()[i]);
00340   }
00341 
00342   return *this ;
00343 }
00344 
00345 
00346 //=========================================================================
00347 // Test for equality
00348 bool mbl_sample_stats_1d::operator==(const mbl_sample_stats_1d& s) const
00349 {
00350   return samples_==s.samples_ && use_mvue_==s.use_mvue_;
00351 }
00352 
00353 // =============================================
00354 short mbl_sample_stats_1d::version_no() const
00355 {
00356   return 1;
00357 }
00358 
00359 
00360 //=========================================================================
00361 void mbl_sample_stats_1d::b_write(vsl_b_ostream& bfs) const
00362 {
00363   vsl_b_write(bfs,version_no());
00364   vsl_b_write(bfs,samples_);
00365   vsl_b_write(bfs,stats_1d_);
00366   vsl_b_write(bfs,use_mvue_);
00367 }
00368 
00369 
00370 //=========================================================================
00371 void mbl_sample_stats_1d::b_read(vsl_b_istream& bfs)
00372 {
00373   if (!bfs) return;
00374 
00375   short file_version_no;
00376   vsl_b_read(bfs,file_version_no);
00377 
00378   switch (file_version_no)
00379   {
00380   case 1:
00381     vsl_b_read(bfs,samples_);
00382     vsl_b_read(bfs,stats_1d_);
00383     vsl_b_read(bfs,use_mvue_);
00384     break;
00385   default :
00386     vcl_cerr << "I/O ERROR: mbl_sample_stats_1d::b_read(vsl_b_istream&)\n"
00387              << "           Unknown version number "<< file_version_no << '\n';
00388     bfs.is().clear(vcl_ios::badbit); // Set an unrecoverable IO error on stream
00389     return;
00390   }
00391 }
00392 
00393 
00394 //=========================================================================
00395 void mbl_sample_stats_1d::print_summary(vcl_ostream& os) const
00396 {
00397   os << "mbl_sample_stats_1d: ";
00398   if (samples_.size()==0)
00399   {
00400     os << "No samples.";
00401   }
00402   else
00403   {
00404     os << "mean: "<< mean()
00405        << " use MVUE: "<< use_mvue_
00406        << " sd: "<< sd()
00407        << " ["<<stats_1d_.min()<<','<<stats_1d_.max()<<"] N:"<<samples_.size();
00408   }
00409 }
00410 
00411 
00412 //=========================================================================
00413 // Print all data
00414 void mbl_sample_stats_1d::print_all(vcl_ostream& os,
00415                                     const vcl_string& delim/*="\n"*/) const
00416 {
00417   unsigned nSamples = samples_.size();
00418   for (unsigned i=0; i<nSamples; ++i)
00419   {
00420     os << samples_[i] << delim;
00421   }
00422 }
00423 
00424 
00425 //=========================================================================
00426 vcl_ostream& operator<<(vcl_ostream& os, const mbl_sample_stats_1d& stats)
00427 {
00428   stats.print_summary(os);
00429   return os;
00430 }
00431 
00432 
00433 //=========================================================================
00434 // Stream output operator for class reference
00435 void vsl_print_summary(vcl_ostream& os,const mbl_sample_stats_1d& stats)
00436 {
00437   stats.print_summary(os);
00438 }
00439 
00440 
00441 //=========================================================================
00442 // Print all data
00443 void vsl_print_all(vcl_ostream& os, const mbl_sample_stats_1d& stats)
00444 {
00445   stats.print_all(os);
00446 }
00447 
00448 
00449 //=========================================================================
00450 mbl_sample_stats_1d operator+(const mbl_sample_stats_1d& s1, const mbl_sample_stats_1d& s2)
00451 {
00452   mbl_sample_stats_1d r = s1;
00453   r+=s2;
00454 
00455   return r;
00456 }
00457 
00458 
00459 //=========================================================================
00460 // Binary file stream output operator for class reference
00461 void vsl_b_write(vsl_b_ostream& bfs, const mbl_sample_stats_1d& b)
00462 {
00463   b.b_write(bfs);
00464 }
00465 
00466 
00467 //=========================================================================
00468 // Binary file stream input operator for class reference
00469 void vsl_b_read(vsl_b_istream& bfs, mbl_sample_stats_1d& b)
00470 {
00471   b.b_read(bfs);
00472 }
00473