contrib/mul/mbl/mbl_stats_1d.cxx
Go to the documentation of this file.
00001 // This is mul/mbl/mbl_stats_1d.cxx
00002 #include "mbl_stats_1d.h"
00003 //:
00004 // \file
00005 // \brief Simple statistics on a 1D variable.
00006 // \author Tim Cootes
00007 
00008 #include <vcl_cmath.h>
00009 #include <vcl_iostream.h>
00010 
00011 mbl_stats_1d::mbl_stats_1d()
00012 {
00013   clear();
00014 }
00015 
00016 mbl_stats_1d::mbl_stats_1d(const vcl_vector<double>& observations)
00017 {
00018   clear();
00019   vcl_vector<double>::const_iterator it;
00020   for (it=observations.begin(); it != observations.end(); ++it)
00021   {
00022     obs(*it);
00023   }
00024 }
00025 
00026 void mbl_stats_1d::clear()
00027 {
00028   n_obs_ = 0;
00029   w_obs_ = 0;
00030   sum_ = 0;
00031   sum_sq_ = 0;
00032 }
00033 
00034 void mbl_stats_1d::obs(double v)
00035 {
00036   if (n_obs_ == 0)
00037   {
00038     min_v_ = v;
00039     max_v_ = v;
00040     sum_ = v;
00041     sum_sq_ = v * v;
00042     w_obs_=1.0;
00043     n_obs_=1;
00044     return;
00045   }
00046 
00047   if (v<min_v_) min_v_ = v;
00048   if (v>max_v_) max_v_ = v;
00049   sum_ += v;
00050   sum_sq_ += v * v;
00051   n_obs_++;
00052   w_obs_++;
00053 }
00054 
00055 void mbl_stats_1d::obs(double v, double weight)
00056 {
00057   if (n_obs_ == 0)
00058   {
00059     min_v_ = v;
00060     max_v_ = v;
00061     sum_ = v * weight;
00062     sum_sq_ = v * v * weight;
00063     w_obs_=weight;
00064     n_obs_=1;
00065     return;
00066   }
00067 
00068   if (v<min_v_) min_v_ = v;
00069   if (v>max_v_) max_v_ = v;
00070   sum_ += v * weight;
00071   sum_sq_ += v * v * weight;
00072   w_obs_+=weight;
00073   n_obs_++;
00074 }
00075 
00076 double mbl_stats_1d::mean() const
00077 {
00078   if (n_obs_==0) return 0;
00079   else return sum_/w_obs_;
00080 }
00081 
00082 double mbl_stats_1d::variance() const
00083 {
00084   if (n_obs_==0) return 0;
00085 
00086   double mean_v = mean();
00087   return sum_sq_/w_obs_ - mean_v * mean_v;
00088 }
00089 
00090 
00091 double mbl_stats_1d::sd() const
00092 {
00093   if (n_obs_==0) return 0;
00094 
00095   double var_v = variance();
00096   // Use of numerically dodgy Sum{x^2} - {Sum x}^2
00097   // can return negative numbers.
00098   if (var_v<0) var_v=0;
00099   return vcl_sqrt(var_v);
00100 }
00101 
00102 double mbl_stats_1d::stdError() const
00103 {
00104   if (n_obs_==0) return 0;
00105 
00106   double var_v = variance();
00107   return vcl_sqrt(var_v/w_obs_);
00108 }
00109 
00110 double mbl_stats_1d::min() const
00111 {
00112   if (n_obs_==0) return 0;
00113   else return min_v_;
00114 }
00115 
00116 double mbl_stats_1d::max() const
00117 {
00118   if (n_obs_==0) return 0;
00119   else return max_v_;
00120 }
00121 
00122 double mbl_stats_1d::sum() const
00123 {
00124   return sum_;
00125 }
00126 
00127 double mbl_stats_1d::sumSq() const
00128 {
00129   return sum_sq_;
00130 }
00131 
00132 double mbl_stats_1d::rms() const
00133 {
00134   return n_obs_==0 ? -1.0 : vcl_sqrt(sum_sq_/w_obs_);
00135 }
00136 
00137 
00138 mbl_stats_1d& mbl_stats_1d::operator+=(const mbl_stats_1d& s1)
00139 {
00140   sum_ += s1.sum_;
00141   sum_sq_ += s1.sum_sq_;
00142   n_obs_ += s1.n_obs_;
00143   w_obs_ += s1.w_obs_;
00144   if (s1.min()<min_v_) min_v_ = s1.min_v_;
00145   if (s1.max()>max_v_) max_v_ = s1.max_v_;
00146   return *this ;
00147 }
00148 
00149 const double MAX_ERROR = 1.0e-8;
00150 
00151 //: Test for equality
00152 bool mbl_stats_1d::operator==(const mbl_stats_1d& s) const
00153 {
00154   return n_obs_==s.n_obs_ &&
00155          vcl_fabs(w_obs_-s.w_obs_)<MAX_ERROR &&
00156          vcl_fabs(sum_-s.sum_)<MAX_ERROR &&
00157          vcl_fabs(sum_sq_-s.sum_sq_)<MAX_ERROR &&
00158          vcl_fabs(min_v_-s.min_v_)<MAX_ERROR &&
00159          vcl_fabs(max_v_-s.max_v_)<MAX_ERROR;
00160 }
00161 
00162 
00163 void mbl_stats_1d::b_write(vsl_b_ostream& bfs) const
00164 {
00165   const short version = 2;
00166   vsl_b_write(bfs,version);
00167   vsl_b_write(bfs,n_obs_);
00168   if (n_obs_==0) return;
00169   vsl_b_write(bfs,min_v_); vsl_b_write(bfs,max_v_);
00170   vsl_b_write(bfs,sum_); vsl_b_write(bfs,sum_sq_);
00171   vsl_b_write(bfs,w_obs_);
00172 }
00173 
00174 void mbl_stats_1d::b_read(vsl_b_istream& bfs)
00175 {
00176   if (!bfs) return;
00177 
00178   short file_version_no;
00179   vsl_b_read(bfs,file_version_no);
00180 
00181   switch (file_version_no)
00182   {
00183    case 1:
00184     {
00185       int tmp;
00186       vsl_b_read(bfs, tmp);
00187       n_obs_ = static_cast<unsigned>(tmp);
00188     }
00189     if (n_obs_<=0) clear();
00190     else
00191     {
00192       vsl_b_read(bfs,min_v_);
00193       vsl_b_read(bfs,max_v_);
00194       vsl_b_read(bfs,sum_);
00195       vsl_b_read(bfs,sum_sq_);
00196     }
00197     w_obs_ = n_obs_;
00198     break;
00199    case 2:
00200     vsl_b_read(bfs, n_obs_);
00201     if (n_obs_<=0) clear();
00202     else
00203     {
00204       vsl_b_read(bfs,min_v_);
00205       vsl_b_read(bfs,max_v_);
00206       vsl_b_read(bfs,sum_);
00207       vsl_b_read(bfs,sum_sq_);
00208       vsl_b_read(bfs,w_obs_);
00209     }
00210     break;
00211    default:
00212     vcl_cerr << "I/O ERROR: mbl_stats_1d::b_read(vsl_b_istream&)\n"
00213              << "           Unknown version number "<< file_version_no << '\n';
00214     bfs.is().clear(vcl_ios::badbit); // Set an unrecoverable IO error on stream
00215     return;
00216   }
00217 }
00218 
00219 void mbl_stats_1d::print_summary(vcl_ostream& os) const
00220 {
00221   os << "mbl_stats_1d: ";
00222   if (n_obs_==0)
00223     os << "No samples.";
00224   else
00225   {
00226     os << "mean: "<< mean()
00227        << " sd: "<< sd()
00228        << " ["<<min_v_<<','<<max_v_<<"] N:"<<n_obs_;
00229   }
00230 }
00231 
00232 vcl_ostream& operator<<(vcl_ostream& os, const mbl_stats_1d& stats)
00233 {
00234   stats.print_summary(os);
00235   return os;
00236 }
00237 
00238   //: Stream output operator for class reference
00239 void vsl_print_summary(vcl_ostream& os,const mbl_stats_1d& stats)
00240 {
00241   stats.print_summary(os);
00242 }
00243 
00244 mbl_stats_1d operator+(const mbl_stats_1d& s1, const mbl_stats_1d& s2)
00245 {
00246   mbl_stats_1d r = s1;
00247   r+=s2;
00248 
00249   return r;
00250 }
00251 
00252   //: Binary file stream output operator for class reference
00253 void vsl_b_write(vsl_b_ostream& bfs, const mbl_stats_1d& b)
00254 {
00255   b.b_write(bfs);
00256 }
00257 
00258   //: Binary file stream input operator for class reference
00259 void vsl_b_read(vsl_b_istream& bfs, mbl_stats_1d& b)
00260 {
00261   b.b_read(bfs);
00262 }