26 #ifndef O2SCL_FIT_BAYES_H 27 #define O2SCL_FIT_BAYES_H 31 #include <boost/numeric/ublas/vector.hpp> 33 #include <o2scl/hist.h> 34 #include <o2scl/rng_gsl.h> 35 #include <o2scl/interp.h> 36 #include <o2scl/fit_base.h> 37 #include <o2scl/expval.h> 38 #include <o2scl/mcarlo.h> 39 #include <o2scl/mcarlo_vegas.h> 40 #include <o2scl/vector.h> 42 #ifndef DOXYGEN_NO_O2NS 49 template<
class vec_t=boost::numeric::ublas::vector<
double> >
57 virtual double operator()(
size_t npar,
const vec_t &par) {
76 template<
class fit_func_t=fit_funct,
class multi_func_t=uniform_prior<>,
77 class vec_t=boost::numeric::ublas::vector<
double> >
class fit_bayes {
89 virtual double integrand(
size_t npar,
const vec_t &par) {
90 return weight_fun(lndat,*lxdat,*lydat,*lyerr,npar,par)*(*pri)(npar,par);
140 virtual void evidence(
size_t ndat, vec_t &xdat, vec_t &ydat,
141 vec_t &yerr,
size_t npar, vec_t &plo2,
142 vec_t &phi2, multi_func_t &prior_fun,
143 double &evi,
double &err) {
150 std::bind(std::mem_fn<
double(
size_t,
const ubvector &)>
152 this,std::placeholders::_1,std::placeholders::_2);
153 def_inte.
minteg_err(mfm,npar,plo2,phi2,evi,err);
161 const vec_t &ydat,
const vec_t &yerr,
162 size_t npar,
const vec_t &par) {
165 for(
size_t i=0;i<ndat;i++) {
166 weight*=exp(-pow((*ff)(npar,par,xdat[i])-ydat[i],2.0)/
167 2.0/yerr[i]/yerr[i]);
178 virtual int fit(
size_t ndat, vec_t &xdat, vec_t &ydat, vec_t &yerr,
179 size_t npar, vec_t &plo2, vec_t &pmax, vec_t &phi2,
180 vec_t &plo_err, vec_t &pmax_err, vec_t &phi_err,
181 fit_func_t &fitfun, multi_func_t &prior_fun) {
186 double weight, weight_old;
190 par_old.resize(npar);
194 for(
size_t i=0;i<npar;i++) {
195 par_old[i]=(plo2[i]+phi2[i])/2.0;
197 weight_old=weight_fun(ndat,xdat,ydat,yerr,npar,par_old)*
198 (*pri)(npar,par_old);
199 if (weight_old==0.0) {
205 std::vector<hist> harr(npar);
206 for(
size_t i=0;i<npar;i++) {
208 harr[i].clear_wgts();
212 size_t meas_size=((size_t)(((
double)n_iter)/((double)nmeas)));
213 size_t imeas=meas_size-1;
217 ubvector_int nonzero_bins(npar);
220 for(
size_t it=0;it<n_iter+n_warm_up;it++) {
223 for (
size_t i=0;i<npar;i++) {
224 par[i]=gr.
random()*(phi2[i]-plo2[i])+plo2[i];
226 weight=weight_fun(ndat,xdat,ydat,yerr,npar,par)*(*pri)(npar,par);
230 if (weight>0.0 && (weight>weight_old || r<weight/weight_old)) {
231 for (
size_t i=0;i<npar;i++) {
240 for(
size_t i=0;i<npar;i++) {
242 if (par_old[i]>plo2[i] && par_old[i]<phi2[i]) {
243 harr[i].update(par_old[i]);
247 if (it-n_warm_up==imeas) {
251 std::vector<double> xlo, xhi, xmax;
252 for(
size_t i=0;i<npar;i++) {
255 std::vector<double> edge_x, edge_y;
256 edge_x.push_back(2.0*harr[i].get_rep_i(0)-
257 harr[i].get_rep_i(1));
258 edge_y.push_back(0.0);
259 for(
size_t j=0;j<hsize;j++) {
260 if (harr[i].get_wgt_i(j)>0.0) {
263 edge_x.push_back(harr[i].get_rep_i(j));
264 edge_y.push_back(harr[i].get_wgt_i(j));
266 edge_x.push_back(2.0*harr[i].get_rep_i(hsize-1)-
267 harr[i].get_rep_i(hsize-2));
268 edge_y.push_back(0.0);
272 edge_y[edge_y.size()-2]/=2.0;
278 std::vector<double> locs;
282 double lmin=*min_element(locs.begin(),locs.end());
283 double lmax=*max_element(locs.begin(),locs.end());
287 size_t imax=vector_max_index<std::vector<double>,
double>
290 (edge_x[imax-1],edge_x[imax],edge_x[imax+1],
291 edge_y[imax-1],edge_y[imax],edge_y[imax+1]));
300 for(
size_t i=0;i<npar;i++) {
301 harr[i].clear_wgts();
313 ubvector sd1(npar), sd2(npar), sd3(npar);
318 for(
size_t i=0;i<npar;i++) {
319 if (((
double)nonzero_bins[i])<((
double)nmeas)*2.5) {
320 O2SCL_ERR2(
"Not enough data for accurate parameter ",
335 virtual int fit_hist(
size_t ndat, vec_t &xdat, vec_t &ydat, vec_t &yerr,
336 size_t npar, vec_t &plo2, vec_t &phi2,
337 std::vector<hist> &par_hist, fit_func_t &fitfun,
338 multi_func_t &prior_fun) {
343 double weight, weight_old;
347 par_old.resize(npar);
351 for(
size_t i=0;i<npar;i++) {
352 par_old[i]=(plo2[i]+phi2[i])/2.0;
354 weight_old=weight_fun(ndat,xdat,ydat,yerr,npar,par_old)*
355 (*pri)(npar,par_old);
356 if (weight_old==0.0) {
362 par_hist.resize(npar);
363 std::vector<hist> harr(npar);
364 for(
size_t i=0;i<npar;i++) {
366 harr[i].clear_wgts();
368 par_hist[i].clear_wgts();
372 size_t meas_size=((size_t)(((
double)n_iter)/((double)nmeas)));
373 size_t imeas=meas_size-1;
375 ubmatrix mat_tmp(npar,hsize);
378 for(
size_t it=0;it<n_iter+n_warm_up;it++) {
381 for (
size_t i=0;i<npar;i++) {
382 par[i]=gr.
random()*(phi2[i]-plo2[i])+plo2[i];
384 weight=weight_fun(ndat,xdat,ydat,yerr,npar,par)*(*pri)(npar,par);
388 if (weight>0.0 && (weight>weight_old || r<weight/weight_old)) {
389 for (
size_t i=0;i<npar;i++) {
398 for(
size_t i=0;i<npar;i++) {
400 if (par_old[i]>plo2[i] && par_old[i]<phi2[i]) {
401 harr[i].update(par_old[i]);
405 if (it-n_warm_up==imeas) {
407 for(
size_t i=0;i<npar;i++) {
408 for(
size_t j=0;j<hsize;j++) {
409 mat_tmp(i,j)=harr[i].get_wgt_i(j);
412 hist_ev.
add(mat_tmp);
415 for(
size_t i=0;i<npar;i++) {
416 harr[i].clear_wgts();
427 ubmatrix avg(npar,hsize), std_dev(npar,hsize), avg_err(npar,hsize);
429 for(
size_t i=0;i<npar;i++) {
430 for(
size_t j=0;j<hsize;j++) {
431 par_hist[i].set_wgt_i(j,avg(i,j));
440 #ifndef DOXYGEN_NO_O2NS std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct
Multi-dimensional function typedef.
multi_func_t * pri
Prior distribution.
Matrix expectation value.
size_t hsize
Histogram size (default 20)
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
void current_avg(vec_t &avg, vec2_t &std_dev, vec3_t &avg_err)
Report current average, standard deviation, and the error in the average.
int vector_invert_enclosed_sum(double sum, size_t n, vec_t &x, vec2_t &y, double &lev, int boundaries=0, int verbose=0, bool err_on_fail=true)
Compute the endpoints which enclose the regions whose integral is equal to sum.
virtual int minteg_err(func_t &func, size_t ndim, const vec_t &a, const vec_t &b, double &res, double &err)
Integrate function func from x=a to x=b.
invalid argument supplied by user
virtual int fit(size_t ndat, vec_t &xdat, vec_t &ydat, vec_t &yerr, size_t npar, vec_t &plo2, vec_t &pmax, vec_t &phi2, vec_t &plo_err, vec_t &pmax_err, vec_t &phi_err, fit_func_t &fitfun, multi_func_t &prior_fun)
Fit ndat data points in xdat and ydat with errors yerr to function fitfun with npar parameters...
double random()
Return a random number in .
virtual double integrand(size_t npar, const vec_t &par)
The integrand for the evidence.
size_t lndat
Number of data points.
rng_gsl gr
Random number generator.
fit_func_t * ff
User-specified function.
virtual double weight_fun(size_t ndat, const vec_t &xdat, const vec_t &ydat, const vec_t &yerr, size_t npar, const vec_t &par)
The weight function (based on a distribution)
mcarlo_vegas def_inte
Default Monte Carlo integrator.
data_t quadratic_extremum_x(const data_t x1, const data_t x2, const data_t x3, const data_t y1, const data_t y2, const data_t y3)
Return the x value of the extremum of a quadratic defined by three pairs.
void current_avg(mat_t &avg, mat2_t &std_dev, mat3_t &avg_err)
Report current average, standard deviation, and the error in the average.
double vector_integ_interp(size_t n, ovec_t &v, size_t interp_type)
Integral of a vector from interpolation object.
size_t nmeas
Number of measurements (default 20)
virtual int fit_hist(size_t ndat, vec_t &xdat, vec_t &ydat, vec_t &yerr, size_t npar, vec_t &plo2, vec_t &phi2, std::vector< hist > &par_hist, fit_func_t &fitfun, multi_func_t &prior_fun)
Desc.
Vector expectation value.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
size_t n_warm_up
Number of warmup iterations (default 100)
size_t n_iter
Number of total iterations (default 1000)
Multidimensional integration using Vegas Monte Carlo (GSL)
Random number generator (GSL)
Fit a function to data using Bayesian methods.
void add(vec_t &val)
Add measurement of value val.
void add(mat_t &val)
Add measurement of value val.
void vector_find_level(double level, size_t n, vec_t &x, vec2_t &y, std::vector< double > &locs)
Perform inverse linear interpolation.
virtual void evidence(size_t ndat, vec_t &xdat, vec_t &ydat, vec_t &yerr, size_t npar, vec_t &plo2, vec_t &phi2, multi_func_t &prior_fun, double &evi, double &err)
Compute the evidence.