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) {
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);
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;
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]);
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;
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]);
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();
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