Go to the documentation of this file.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 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
00097
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
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);
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
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
00253 void vsl_b_write(vsl_b_ostream& bfs, const mbl_stats_1d& b)
00254 {
00255 b.b_write(bfs);
00256 }
00257
00258
00259 void vsl_b_read(vsl_b_istream& bfs, mbl_stats_1d& b)
00260 {
00261 b.b_read(bfs);
00262 }