diff options
Diffstat (limited to 'src/vchisq.hpp')
| -rw-r--r-- | src/vchisq.hpp | 137 |
1 files changed, 137 insertions, 0 deletions
diff --git a/src/vchisq.hpp b/src/vchisq.hpp new file mode 100644 index 0000000..8cd8481 --- /dev/null +++ b/src/vchisq.hpp @@ -0,0 +1,137 @@ +/**
+ \file vchisq.hpp
+ \brief chi-square statistic
+ \author Junhua Gu
+ */
+
+#ifndef VCHI_SQ_HPP
+#define VCHI_SQ_HPP
+
+#define OPT_HEADER
+
+#include <core/fitter.hpp>
+#include <iostream>
+#include <vector>
+#include <misc/optvec.hpp>
+#include <cmath>
+
+using std::cerr;
+using std::endl;
+
+namespace opt_utilities
+{
+ template<typename T>
+ class vchisq
+ :public statistic<std::vector<T>,std::vector<T>,std::vector<T>,T,std::string>
+ {
+ private:
+ bool verb;
+ int n;
+ bool limit_bound;
+ typedef std::vector<T> Tp;
+
+ vchisq<T>* do_clone()const
+ {
+ return new vchisq<T>(*this);
+ }
+
+ const char* do_get_type_name()const
+ {
+ return "chi^2 statistic";
+ }
+
+ public:
+ void verbose(bool v)
+ {
+ verb=v;
+ }
+
+ void set_limit()
+ {
+ limit_bound=true;
+ }
+
+ void clear_limit()
+ {
+ limit_bound=false;
+ }
+ public:
+ vchisq()
+ :verb(false),limit_bound(false)
+ {}
+
+ T do_eval(const std::vector<T>& p)
+ {
+ if(limit_bound)
+ {
+ if(!this->get_fitter().get_model().meets_constraint(p))
+ {
+ return 1e99;
+ }
+ }
+ T result(0);
+ std::vector<float> vx;
+ std::vector<float> vy;
+ std::vector<float> vye;
+ std::vector<float> my;
+ float x1=1e99,x2=-1e99,y1=1e99,y2=-1e99;
+ if(verb)
+ {
+ n++;
+ if(n%100==0)
+ {
+ vx.resize(this->get_data_set().get_data(0).get_y().size());
+ vy.resize(this->get_data_set().get_data(0).get_y().size());
+ vye.resize(this->get_data_set().get_data(0).get_y().size());
+ my.resize(this->get_data_set().get_data(0).get_y().size());
+ }
+ }
+ for(int i=(this->get_data_set()).size()-1;i>=0;--i)
+ {
+ const std::vector<double> y_model(this->eval_model(this->get_data_set().get_data(i).get_x(),p));
+ const std::vector<double>& y=this->get_data_set().get_data(i).get_y();
+ const std::vector<double>& ye=this->get_data_set().get_data(i).get_y_lower_err();
+ for(size_t j=0;j<y.size();++j)
+ {
+ double chi=(y_model[j]-y[j])/ye[j];
+ result+=chi*chi;
+ }
+
+ if(verb&&n%100==0)
+ {
+ for(size_t j=0;j<y.size();++j)
+ {
+ vx.at(j)=((this->get_data_set().get_data(i).get_x().at(j)+this->get_data_set().get_data(i).get_x().at(j+1))/2.);
+ vy.at(j)=(y[j]);
+ vye.at(j)=ye[j];
+ my.at(j)=(y_model[j]);
+ x1=std::min(vx.at(j),x1);
+ y1=std::min(vy.at(j),y1);
+ x2=std::max(vx.at(j),x2);
+ y2=std::max(vy.at(j),y2);
+ vye[j]=log10(vy[j]+vye[j])-log10(vy[j]);
+ vx[j]=log10(vx[j]);
+ vy[j]=log10(vy[j]);
+ my[j]=log10(my[j]);
+ }
+ }
+ }
+ if(verb)
+ {
+ if(n%100==0)
+ {
+ cerr<<result<<"\t";
+ for(size_t i=0;i<get_size(p);++i)
+ {
+ cerr<<get_element(p,i)<<",";
+ }
+ cerr<<endl;
+ }
+ }
+
+ return result;
+ }
+ };
+
+}
+#endif
|