aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--misc/bootstrap.hpp80
-rw-r--r--misc/optvec.hpp32
-rw-r--r--statistics/chisq.hpp18
-rw-r--r--test/test_fitter.cpp4
-rw-r--r--vmodels/beta1d.hpp64
-rw-r--r--vmodels/constant.hpp54
-rw-r--r--vmodels/poly1d.hpp78
7 files changed, 302 insertions, 28 deletions
diff --git a/misc/bootstrap.hpp b/misc/bootstrap.hpp
index f2ba525..88332ed 100644
--- a/misc/bootstrap.hpp
+++ b/misc/bootstrap.hpp
@@ -3,11 +3,12 @@
*/
-#ifndef BOOT_STRIP
-#define BOOT_STRIP
+#ifndef BOOT_STRAP
+#define BOOT_STRAP
#define OPT_HEADER
#include <core/fitter.hpp>
#include <vector>
+#include <list>
#include <cstdlib>
#include <iostream>
#include <utility>
@@ -16,6 +17,18 @@
using std::cerr;
namespace opt_utilities
{
+ template <typename Ty>
+ inline Ty rand_norm(Ty y0,Ty y_err)
+ {
+ Ty y;
+ do
+ {
+ y=(rand()/(Ty)RAND_MAX-(Ty).5)*(10*y_err)+y0;
+ }
+ while(rand()/(Ty)RAND_MAX>exp(-(y-y0)*(y-y0)/(y_err*y_err)));
+ return y;
+ }
+
/**
\brief using bootstrap method to estimate confidence interval
@@ -29,23 +42,16 @@ namespace opt_utilities
class bootstrap
{
private:
- Ty rand_norm(Ty y0,Ty y_err)const
- {
- Ty y;
- do
- {
- y=(rand()/(Ty)RAND_MAX-(Ty).5)*(10*y_err)+y0;
- }
- while(rand()/(Ty)RAND_MAX>exp(-(y-y0)*(y-y0)/(y_err*y_err)));
- return y;
-
- }
public:
- std::vector<Tp> param_pool;
+ typedef typename std::list<Tp>::const_iterator param_iterator;
+ public:
+ std::list<Tp> param_pool;
default_data_set<Ty,Tx> current_data_set;
default_data_set<Ty,Tx> origin_data_set;
fitter<Ty,Tx,Tp,Ts,Tstr>* p_fitter;
Tp origin_param;
+ mutable bool bstop;
+
public:
/**
@@ -83,11 +89,12 @@ namespace opt_utilities
*/
void reset()
{
- if(p_fitter!=NULL)
+ if(p_fitter!=0)
{
p_fitter->load_data(origin_data_set);
p_fitter->set_param_value(origin_param);
}
+ p_fitter=0;
}
@@ -97,9 +104,10 @@ namespace opt_utilities
*/
void sample(int n)
{
+ bstop=false;
if(p_fitter!=NULL)
{
- for(int i=0;i<n;++i)
+ for(int i=0;i<n&&!bstop;++i)
{
sample();
}
@@ -113,15 +121,45 @@ namespace opt_utilities
}
}
+ void stop()
+ {
+ bstop=true;
+ }
+
/**
get the param of i-th sampling
\param i the order of parameter
\return the parameter
*/
- const Tp& get_param(int i)const
+ const Tp& get_param(int n)const
{
- return param_pool.at(i);
+ //return param_pool.at(i);
+ int cnt=0;
+ for(typename std::list<Tp>::const_iterator i=param_pool.begin();
+ i!=param_pool.end();++i)
+ {
+ if(cnt==n)
+ {
+ return *i;
+ }
+ ++cnt;
+ }
+ throw opt_exception("excesses param_pool size");
+ return *param_pool.begin();
+ }
+
+ typename std::list<Tp>::const_iterator
+ param_begin()const
+ {
+ return param_pool.begin();
+ }
+
+ typename std::list<Tp>::const_iterator
+ param_end()const
+ {
+ return param_pool.end();
}
+
int get_param_pool_size()const
{
@@ -140,7 +178,7 @@ namespace opt_utilities
{
data<Ty,Tx> d;
d=origin_data_set.get_data(i);
- d.set_y(rand_norm(d.get_y(),(d.get_y_upper_err()+d.get_y_lower_err())/2));
+ d.set_y(rand_norm(d.get_y(),(d.get_y_upper_err()+d.get_y_lower_err())/2.));
current_data_set.add_data(d);
}
p_fitter->load_data(current_data_set);
@@ -176,7 +214,7 @@ namespace opt_utilities
//sample();
std::vector<typename element_type_trait<Tp>::element_type> _tmp;
int order=p_fitter->get_param_order(param_name);
- for(typename std::vector<Tp>::iterator i=param_pool.begin();
+ for(typename std::list<Tp>::iterator i=param_pool.begin();
i!=param_pool.end();++i)
{
_tmp.push_back((*i)[order]);
@@ -186,7 +224,7 @@ namespace opt_utilities
typename std::vector<typename element_type_trait<Tp>::element_type>::iterator>
itv=equal_range(_tmp.begin(),_tmp.end(),origin_param[order]);
int current_param_position=itv.second-_tmp.begin();
- std::cerr<<_tmp.size()<<std::endl;
+ //std::cerr<<_tmp.size()<<std::endl;
return std::pair<typename element_type_trait<Tp>::element_type,
typename element_type_trait<Tp>::element_type>(
_tmp.at((int)((1-level)*current_param_position)),
diff --git a/misc/optvec.hpp b/misc/optvec.hpp
index cc85ccf..b259029 100644
--- a/misc/optvec.hpp
+++ b/misc/optvec.hpp
@@ -25,8 +25,8 @@ namespace opt_utilities
optvec& operator=(const optvec& rhs)
{
- //dynamic_cast<std::vector<T>&>(*this).operator=(rhs);
- std::vector<T>::operator=(rhs);
+ static_cast<std::vector<T>&>(*this).operator=(rhs);
+ //std::cerr<<"rhs.size="<<rhs.size()<<std::endl;
return *this;
}
public:
@@ -208,8 +208,8 @@ namespace opt_utilities
optvec<T> operator/(const T& x1,const optvec<T>& x2)
{
- optvec<T> result(x1.size());
- for(size_t i=0;i!=x1.size();++i)
+ optvec<T> result(x2.size());
+ for(size_t i=0;i!=x2.size();++i)
{
result[i]=x1/x2[i];
}
@@ -260,6 +260,30 @@ namespace opt_utilities
return false;
}
+ template <typename Ty>
+ inline Ty randn(Ty y0,Ty y_err)
+ {
+ Ty y;
+ do
+ {
+ y=(rand()/(Ty)RAND_MAX-(Ty).5)*(10*y_err)+y0;
+ }
+ while(rand()/(Ty)RAND_MAX>exp(-(y-y0)*(y-y0)/(y_err*y_err)));
+ return y;
+ }
+
+
+ template <typename T>
+ optvec<T> rand_norm(const optvec<T>& y0, const optvec<T>& y_err)
+ {
+ optvec<T> result(y0.size());
+ for(int i=0;i<y0.size();++i)
+ {
+ result[i]=randn(y0[i],y_err[i]);
+ }
+ return result;
+ }
+
}
diff --git a/statistics/chisq.hpp b/statistics/chisq.hpp
index 339ab6f..c460394 100644
--- a/statistics/chisq.hpp
+++ b/statistics/chisq.hpp
@@ -247,14 +247,26 @@ namespace opt_utilities
Ts result(0);
for(int i=(this->get_data_set()).size()-1;i>=0;--i)
{
- Ty chi=(this->get_data_set().get_data(i).get_y()-eval_model(this->get_data_set().get_data(i).get_x(),p))/this->get_data_set().get_data(i).get_y_upper_err();
+ Ty chi(this->get_data_set().get_data(0).get_y().size());
+ for(int j=0;j<chi.size();++j)
+ {
+ Ty model_y(eval_model(this->get_data_set().get_data(i).get_x(),p));
+ if(model_y[j]>this->get_data_set().get_data(i).get_y()[j])
+ {
+ chi[j]=(this->get_data_set().get_data(i).get_y()[j]-model_y[j])/this->get_data_set().get_data(i).get_y_upper_err()[j];
+ }
+ else
+ {
+ chi[j]=(this->get_data_set().get_data(i).get_y()[j]-model_y[j])/this->get_data_set().get_data(i).get_y_lower_err()[j];
+ }
+ }
result+=sum(chi*chi);
-
+
}
return result;
}
};
-
+
}
diff --git a/test/test_fitter.cpp b/test/test_fitter.cpp
index 5e5b139..226d4ce 100644
--- a/test/test_fitter.cpp
+++ b/test/test_fitter.cpp
@@ -1,4 +1,5 @@
#include <core/fitter.hpp>
+#include <misc/bootstrap.hpp>
#include <methods/powell/powell_method.hpp>
#include <statistics/chisq.hpp>
#include <models/nbeta1d.hpp>
@@ -36,4 +37,7 @@ int main(int argc,char* argv[])
{
cout<<f.get_param_info(i).get_name()<<"="<<f.get_param_info(i).get_value()<<endl;
}
+ bootstrap<double,double,vector<double>,double,std::string> bst;
+ bst.set_fitter(f);
+ bst.sample(1000);
}
diff --git a/vmodels/beta1d.hpp b/vmodels/beta1d.hpp
new file mode 100644
index 0000000..42be3fd
--- /dev/null
+++ b/vmodels/beta1d.hpp
@@ -0,0 +1,64 @@
+#ifndef VBETA_MODEL_H_
+#define VBETA_MODEL_H_
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <misc/optvec.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename T>
+ class beta1d
+ :public model<optvec<T>,optvec<T>,optvec<T>,std::string>
+ {
+ typedef optvec<T> Tv;
+ private:
+ beta1d<T>* do_clone()const
+ {
+ return new beta1d<T>(*this);
+ }
+
+ const char* do_get_type_name()const
+ {
+ return "1d surface brightness beta model";
+ }
+ public:
+ beta1d()
+ {
+ this->push_param_info(param_info<Tv>("S0",1));
+ this->push_param_info(param_info<Tv>("rc",10));
+ this->push_param_info(param_info<Tv>("beta",2./3.));
+ this->push_param_info(param_info<Tv>("bkg",0));
+ }
+
+ public:
+ Tv do_eval(const Tv& x,const Tv& param)
+ {
+ Tv result(x.size());
+ T S0=get_element(param,0);
+ T r_c=get_element(param,1);
+ T beta=get_element(param,2);
+ T bkg=get_element(param,3);
+
+ //return x*get_element(param,0)+get_element(param,1);
+ for(size_t i=0;i!=x.size();++i)
+ {
+
+ result[i]=bkg+S0*pow(1+(x[i]*x[i])/(r_c*r_c),-3*beta+static_cast<T>(.5));
+
+ }
+ return result;
+ }
+
+ private:
+ std::string do_get_information()const
+ {
+ return "";
+ }
+ };
+}
+
+
+
+#endif
+//EOF
diff --git a/vmodels/constant.hpp b/vmodels/constant.hpp
new file mode 100644
index 0000000..c4044b0
--- /dev/null
+++ b/vmodels/constant.hpp
@@ -0,0 +1,54 @@
+#ifndef VCONSTANT_MODEL_H_
+#define VCONSTANT_MODEL_H_
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <misc/optvec.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+ template <typename T>
+ class constant
+ :public model<optvec<T>,optvec<T>,optvec<T>,std::string>
+ {
+ typedef optvec<T> Tv;
+ private:
+ constant<T>* do_clone()const
+ {
+ return new constant<T>(*this);
+ }
+ const char* do_get_type_name()const
+ {
+ return "constant";
+ }
+ public:
+ constant()
+ {
+ this->push_param_info(param_info<Tv>("c",1));
+ }
+
+ public:
+ Tv do_eval(const Tv& x,const Tv& param)
+ {
+ //return x*param[0]+param[1];
+ Tv result(x.size());
+ for(int i=0;i<result.size();++i)
+ {
+ result[i]=param[0];
+ }
+ return result;
+ }
+
+ private:
+ std::string do_get_information()const
+ {
+ return "Constant\n"
+ "y=C\n";
+ }
+ };
+}
+
+
+
+#endif
+//EOF
diff --git a/vmodels/poly1d.hpp b/vmodels/poly1d.hpp
new file mode 100644
index 0000000..adaf03b
--- /dev/null
+++ b/vmodels/poly1d.hpp
@@ -0,0 +1,78 @@
+#ifndef POLY_MODEL_H_
+#define POLY_MODEL_H_
+#define OPT_HEADER
+#include <core/fitter.hpp>
+#include <cmath>
+#include <sstream>
+#include <cassert>
+
+namespace opt_utilities
+{
+ template <typename T,int n>
+ class poly1d
+ :public model<optvec<T>,optvec<T>,optvec<T>,std::string>
+ {
+ typedef optvec<T> Tv;
+ private:
+ poly1d<T,n>* do_clone()const
+ {
+ return new poly1d<T,n>(*this);
+ }
+
+ const char* do_get_type_name()const
+ {
+ return "polynomial";
+ }
+ public:
+ poly1d()
+ {
+ assert(n>=0);
+ for(int i=0;i<=n;++i)
+ {
+ std::ostringstream ostr;
+ ostr<<"a"<<i;
+ this->push_param_info(param_info<Tv>(ostr.str().c_str(),1));
+ }
+
+ }
+
+ public:
+ Tv do_eval(const Tv& x,const Tv& param)
+ {
+ // return x*get_element(param,0)+get_element(param,1);
+ Tv result(x.size());
+ for(int m=0;m<result.size();++m)
+ {
+ result[m]=0;
+ for(int i=0;i<=n;++i)
+ {
+ typename Tv::value_type xn(1);
+ for(int j=0;j<i;++j)
+ {
+ xn*=x[m];
+ }
+ result[m]+=get_element(param,i)*xn;
+ }
+ }
+ return result;
+ }
+
+ private:
+ std::string do_get_information()const
+ {
+ std::ostringstream ostr;
+ ostr<<n<<"-order polynorminal model\n";
+ for(int i=0;i<n;++i)
+ {
+ ostr<<"a"<<i<<"*x^"<<i<<"+";
+ }
+ ostr<<"a"<<n<<"*x^"<<n;
+ return ostr.str();
+ }
+ };
+}
+
+
+
+#endif
+//EOF
Copyright © 2017-2018 Aaron LI. All Rights Reserved.