Go to the documentation of this file.00001 #include "mbl_sample_stats_1d.h"
00002
00003
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
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
00137 assert(q>=0.0 && q<=1.0);
00138 assert(n>0);
00139
00140
00141 double float_index = (n-1)*q;
00142
00143
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
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
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
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
00243
00244
00245
00246
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
00272
00273
00274
00275
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
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
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);
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
00414 void mbl_sample_stats_1d::print_all(vcl_ostream& os,
00415 const vcl_string& delim) 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
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
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
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
00469 void vsl_b_read(vsl_b_istream& bfs, mbl_sample_stats_1d& b)
00470 {
00471 b.b_read(bfs);
00472 }
00473