00001
00002 #include "mbl_stats_1d.h"
00003
00004
00005
00006
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 sum_ = 0;
00030 sum_sq_ = 0;
00031 }
00032
00033 void mbl_stats_1d::obs(double v)
00034 {
00035 if (n_obs_ == 0)
00036 {
00037 min_v_ = v;
00038 max_v_ = v;
00039 sum_ = v;
00040 sum_sq_ = v * v;
00041 n_obs_++;
00042 return;
00043 }
00044
00045 if (v<min_v_) min_v_ = v;
00046 if (v>max_v_) max_v_ = v;
00047 sum_ += v;
00048 sum_sq_ += v * v;
00049 n_obs_++;
00050 }
00051
00052 int mbl_stats_1d::nObs() const
00053 {
00054 return n_obs_;
00055 }
00056
00057 double mbl_stats_1d::mean() const
00058 {
00059 if (n_obs_==0) return 0;
00060 else return sum_/n_obs_;
00061 }
00062
00063 double mbl_stats_1d::variance() const
00064 {
00065 if (n_obs_==0) return 0;
00066
00067 double mean_v = mean();
00068 return sum_sq_/n_obs_ - mean_v * mean_v;
00069 }
00070
00071
00072 double mbl_stats_1d::sd() const
00073 {
00074 if (n_obs_==0) return 0;
00075
00076 double var_v = variance();
00077 return vcl_sqrt(var_v);
00078 }
00079
00080 double mbl_stats_1d::stdError() const
00081 {
00082 if (n_obs_==0) return 0;
00083
00084 double var_v = variance();
00085 return vcl_sqrt(var_v/n_obs_);
00086 }
00087
00088 double mbl_stats_1d::min() const
00089 {
00090 if (n_obs_==0) return 0;
00091 else return min_v_;
00092 }
00093
00094 double mbl_stats_1d::max() const
00095 {
00096 if (n_obs_==0) return 0;
00097 else return max_v_;
00098 }
00099
00100 double mbl_stats_1d::sum() const
00101 {
00102 return sum_;
00103 }
00104
00105 double mbl_stats_1d::sumSq() const
00106 {
00107 return sum_sq_;
00108 }
00109
00110 double mbl_stats_1d::rms() const
00111 {
00112 return n_obs_==0 ? -1.0 : vcl_sqrt(sum_sq_/n_obs_);
00113 }
00114
00115
00116 mbl_stats_1d& mbl_stats_1d::operator+=(const mbl_stats_1d& s1)
00117 {
00118 sum_ += s1.sum();
00119 sum_sq_ += s1.sumSq();
00120 n_obs_ += s1.nObs();
00121 if (s1.min()<min_v_) min_v_ = s1.min();
00122 if (s1.max()>max_v_) max_v_ = s1.max();
00123 return *this ;
00124 }
00125
00126 const double MAX_ERROR = 1.0e-8;
00127
00128
00129 bool mbl_stats_1d::operator==(const mbl_stats_1d& s) const
00130 {
00131 return n_obs_==s.nObs() &&
00132 vcl_fabs(sum_-s.sum())<MAX_ERROR &&
00133 vcl_fabs(sum_sq_-s.sumSq())<MAX_ERROR &&
00134 vcl_fabs(min_v_-s.min())<MAX_ERROR &&
00135 vcl_fabs(max_v_-s.max())<MAX_ERROR;
00136 }
00137
00138
00139 short mbl_stats_1d::version_no() const
00140 {
00141 return 1;
00142 }
00143
00144 void mbl_stats_1d::b_write(vsl_b_ostream& bfs) const
00145 {
00146 vsl_b_write(bfs,version_no());
00147 vsl_b_write(bfs,n_obs_);
00148 if (n_obs_==0) return;
00149 vsl_b_write(bfs,min_v_); vsl_b_write(bfs,max_v_);
00150 vsl_b_write(bfs,sum_); vsl_b_write(bfs,sum_sq_);
00151 }
00152
00153 void mbl_stats_1d::b_read(vsl_b_istream& bfs)
00154 {
00155 if (!bfs) return;
00156
00157 short file_version_no;
00158 vsl_b_read(bfs,file_version_no);
00159
00160 switch (file_version_no)
00161 {
00162 case 1:
00163 vsl_b_read(bfs,n_obs_);
00164 if (n_obs_<=0) clear();
00165 else
00166 {
00167 vsl_b_read(bfs,min_v_);
00168 vsl_b_read(bfs,max_v_);
00169 vsl_b_read(bfs,sum_);
00170 vsl_b_read(bfs,sum_sq_);
00171 }
00172 break;
00173 default:
00174 vcl_cerr << "I/O ERROR: mbl_stats_1d::b_read(vsl_b_istream&)\n"
00175 << " Unknown version number "<< file_version_no << '\n';
00176 bfs.is().clear(vcl_ios::badbit);
00177 return;
00178 }
00179 }
00180
00181 void mbl_stats_1d::print_summary(vcl_ostream& os) const
00182 {
00183 os << "mbl_stats_1d: ";
00184 if (n_obs_==0)
00185 os << "No samples.";
00186 else
00187 {
00188 os << "mean: "<< mean()
00189 << " sd: "<< sd()
00190 << " ["<<min_v_<<','<<max_v_<<"] N:"<<n_obs_;
00191 }
00192 }
00193
00194 vcl_ostream& operator<<(vcl_ostream& os, const mbl_stats_1d& stats)
00195 {
00196 stats.print_summary(os);
00197 return os;
00198 }
00199
00200
00201 void vsl_print_summary(vcl_ostream& os,const mbl_stats_1d& stats)
00202 {
00203 stats.print_summary(os);
00204 }
00205
00206 mbl_stats_1d operator+(const mbl_stats_1d& s1, const mbl_stats_1d& s2)
00207 {
00208 mbl_stats_1d r = s1;
00209 r+=s2;
00210
00211 return r;
00212 }
00213
00214
00215 void vsl_b_write(vsl_b_ostream& bfs, const mbl_stats_1d& b)
00216 {
00217 b.b_write(bfs);
00218 }
00219
00220
00221 void vsl_b_read(vsl_b_istream& bfs, mbl_stats_1d& b)
00222 {
00223 b.b_read(bfs);
00224 }