fit_bayes.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2011-2020, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file fit_bayes.h
24  \brief File defining \ref o2scl::fit_bayes and related classes
25 */
26 #ifndef O2SCL_FIT_BAYES_H
27 #define O2SCL_FIT_BAYES_H
28 
29 #include <vector>
30 
31 #include <boost/numeric/ublas/vector.hpp>
32 
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>
41 
42 #ifndef DOXYGEN_NO_O2NS
43 namespace o2scl {
44 #endif
45 
46  /** \brief An unnormalized uniform prior distribution for several
47  variables
48  */
49  template<class vec_t=boost::numeric::ublas::vector<double> >
50  class uniform_prior {
51 
52  public:
53 
54  /** \brief Return the value of the prior at the specified
55  point in parameter space
56  */
57  virtual double operator()(size_t npar, const vec_t &par) {
58  return 1.0;
59  }
60 
61  };
62 
63  /** \brief Fit a function to data using Bayesian methods
64 
65  This class is experimental.
66 
67  This class uses Markov Chain Monte Carlo (MCMC) and
68  marginal estimation to give a probability distribution
69  for parameters in a fit.
70 
71  \future Also make weight_fun() an object of type multi_func_t?
72  \future Offer two ways to do the evidence: direct MC or the
73  interpolation method from SLB13
74  \future Build upon gen_fit_funct instead of fit_funct?
75  */
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 {
78 
79  public:
80 
84 
85  protected:
86 
87  /** \brief The integrand for the evidence
88  */
89  virtual double integrand(size_t npar, const vec_t &par) {
90  return weight_fun(lndat,*lxdat,*lydat,*lyerr,npar,par)*(*pri)(npar,par);
91  }
92 
93  /// User-specified function
94  fit_func_t *ff;
95 
96  /// Number of data points
97  size_t lndat;
98 
99  /// X-values
100  vec_t *lxdat;
101 
102  /// Y-values
103  vec_t *lydat;
104 
105  /// Y-errors
106  vec_t *lyerr;
107 
108  public:
109 
110  fit_bayes() {
111  n_warm_up=100;
112  n_iter=10000;
113  hsize=20;
114  nmeas=20;
115  }
116 
117  /// Number of warmup iterations (default 100)
118  size_t n_warm_up;
119 
120  /// Number of total iterations (default 1000)
121  size_t n_iter;
122 
123  /// Histogram size (default 20)
124  size_t hsize;
125 
126  /// Number of measurements (default 20)
127  size_t nmeas;
128 
129  /// Prior distribution
130  multi_func_t *pri;
131 
132  /// Random number generator
134 
135  /// Default Monte Carlo integrator
137 
138  /** \brief Compute the evidence
139  */
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) {
144  lndat=ndat;
145  lxdat=&xdat;
146  lydat=&ydat;
147  lyerr=&yerr;
148  pri=&prior_fun;
149  multi_funct mfm=
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);
154  return;
155  }
156 
157  /** \brief The weight function (based on a \f$ \chi^2 \f$
158  distribution)
159  */
160  virtual double weight_fun(size_t ndat, const vec_t &xdat,
161  const vec_t &ydat, const vec_t &yerr,
162  size_t npar, const vec_t &par) {
163 
164  double weight=1.0;
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]);
168  }
169  return weight;
170  }
171 
172  /** \brief Fit \c ndat data points in \c xdat and \c ydat with
173  errors \c yerr to function \c fitfun with \c npar parameters.
174 
175  The initial values of the parameters should be specified in
176  \c par.
177  */
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) {
182 
183  ff=&fitfun;
184  pri=&prior_fun;
185 
186  double weight, weight_old;
187 
188  // Create par_old and set up initial weight in weight_old
189  vec_t par_old, par;
190  par_old.resize(npar);
191  par.resize(npar);
192 
193  // Choose the first point as the midpoint of the parameter space
194  for(size_t i=0;i<npar;i++) {
195  par_old[i]=(plo2[i]+phi2[i])/2.0;
196  }
197  weight_old=weight_fun(ndat,xdat,ydat,yerr,npar,par_old)*
198  (*pri)(npar,par_old);
199  if (weight_old==0.0) {
200  O2SCL_ERR2("Initial weight zero in ",
201  "fit_bayes::fit().",exc_einval);
202  }
203 
204  // Set up the histograms for marginal estimation
205  std::vector<hist> harr(npar);
206  for(size_t i=0;i<npar;i++) {
207  harr[i].set_bin_edges(uniform_grid_end<>(plo2[i],phi2[i],hsize));
208  harr[i].clear_wgts();
209  }
210 
211  // Create ev objects for lower and upper limits of each parameter
212  size_t meas_size=((size_t)(((double)n_iter)/((double)nmeas)));
213  size_t imeas=meas_size-1;
214  expval_vector low(npar,nmeas,1);
215  expval_vector high(npar,nmeas,1);
216  expval_vector max(npar,nmeas,1);
217  ubvector_int nonzero_bins(npar);
218 
219  // Main iteration
220  for(size_t it=0;it<n_iter+n_warm_up;it++) {
221 
222  // Select new point in parameter space and evaluate it's weight
223  for (size_t i=0;i<npar;i++) {
224  par[i]=gr.random()*(phi2[i]-plo2[i])+plo2[i];
225  }
226  weight=weight_fun(ndat,xdat,ydat,yerr,npar,par)*(*pri)(npar,par);
227 
228  // Metropolis sampling
229  double r=gr.random();
230  if (weight>0.0 && (weight>weight_old || r<weight/weight_old)) {
231  for (size_t i=0;i<npar;i++) {
232  par_old[i]=par[i];
233  }
234  weight_old=weight;
235  }
236 
237  // If we're not in the warm up phase
238  if (it>=n_warm_up) {
239 
240  for(size_t i=0;i<npar;i++) {
241  // Only update if its in range
242  if (par_old[i]>plo2[i] && par_old[i]<phi2[i]) {
243  harr[i].update(par_old[i]);
244  }
245  }
246 
247  if (it-n_warm_up==imeas) {
248 
249  // From histogram, determine parameter value and error by finding
250  // 68% confidence region
251  std::vector<double> xlo, xhi, xmax;
252  for(size_t i=0;i<npar;i++) {
253 
254  // Correct for non-zero values at the edges
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) {
261  nonzero_bins[i]++;
262  }
263  edge_x.push_back(harr[i].get_rep_i(j));
264  edge_y.push_back(harr[i].get_wgt_i(j));
265  }
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);
269 
270  // Ensure the integral is unchanged
271  edge_y[1]/=2.0;
272  edge_y[edge_y.size()-2]/=2.0;
273 
274  double total=vector_integ_interp(hsize+2,edge_x,edge_y,
275  itp_linear);
276 
277  double lev;
278  std::vector<double> locs;
279  vector_invert_enclosed_sum(0.68*total,hsize+2,edge_x,edge_y,lev);
280  vector_find_level(lev,hsize+2,edge_x,edge_y,locs);
281 
282  double lmin=*min_element(locs.begin(),locs.end());
283  double lmax=*max_element(locs.begin(),locs.end());
284  xlo.push_back(lmin);
285  xhi.push_back(lmax);
286 
287  size_t imax=vector_max_index<std::vector<double>,double>
288  (hsize+2,edge_y);
289  xmax.push_back(quadratic_extremum_x
290  (edge_x[imax-1],edge_x[imax],edge_x[imax+1],
291  edge_y[imax-1],edge_y[imax],edge_y[imax+1]));
292  }
293 
294  // Add the limits to the expval_vector object
295  low.add(xlo);
296  high.add(xhi);
297  max.add(xmax);
298 
299  // Clear the histogram data but leave the grid
300  for(size_t i=0;i<npar;i++) {
301  harr[i].clear_wgts();
302  }
303 
304  // Setup for the next measurement
305  imeas+=meas_size;
306  }
307 
308  }
309 
310  }
311 
312  // Now get parameter values and errors from the expval_vector objects
313  ubvector sd1(npar), sd2(npar), sd3(npar);
314  low.current_avg(plo2,sd1,plo_err);
315  high.current_avg(phi2,sd2,phi_err);
316  max.current_avg(pmax,sd3,pmax_err);
317 
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 ",
321  "estimation in fit_bayes::fit().",exc_efailed);
322  }
323  }
324 
325  return 0;
326  }
327 
328  /** \brief Desc
329 
330  For each measurement block, this function collects the data for
331  all the parameters into 1d histogram objects. Then, at the end
332  of the block, the histogram information is added to a hist
333  object for each parameter.
334  */
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) {
339 
340  ff=&fitfun;
341  pri=&prior_fun;
342 
343  double weight, weight_old;
344 
345  // Create par_old and set up initial weight in weight_old
346  vec_t par_old, par;
347  par_old.resize(npar);
348  par.resize(npar);
349 
350  // Choose the first point as the midpoint of the parameter space
351  for(size_t i=0;i<npar;i++) {
352  par_old[i]=(plo2[i]+phi2[i])/2.0;
353  }
354  weight_old=weight_fun(ndat,xdat,ydat,yerr,npar,par_old)*
355  (*pri)(npar,par_old);
356  if (weight_old==0.0) {
357  O2SCL_ERR2("Initial weight zero in ",
358  "fit_bayes::fit().",exc_einval);
359  }
360 
361  // Set up the histograms for marginal estimation
362  par_hist.resize(npar);
363  std::vector<hist> harr(npar);
364  for(size_t i=0;i<npar;i++) {
365  harr[i].set_bin_edges(uniform_grid_end<>(plo2[i],phi2[i],hsize));
366  harr[i].clear_wgts();
367  par_hist[i].set_bin_edges(uniform_grid_end<>(plo2[i],phi2[i],hsize));
368  par_hist[i].clear_wgts();
369  }
370 
371  // Create ev objects for lower and upper limits of each parameter
372  size_t meas_size=((size_t)(((double)n_iter)/((double)nmeas)));
373  size_t imeas=meas_size-1;
374  expval_matrix hist_ev(npar,hsize,nmeas,1);
375  ubmatrix mat_tmp(npar,hsize);
376 
377  // Main iteration
378  for(size_t it=0;it<n_iter+n_warm_up;it++) {
379 
380  // Select new point in parameter space and evaluate it's weight
381  for (size_t i=0;i<npar;i++) {
382  par[i]=gr.random()*(phi2[i]-plo2[i])+plo2[i];
383  }
384  weight=weight_fun(ndat,xdat,ydat,yerr,npar,par)*(*pri)(npar,par);
385 
386  // Metropolis sampling
387  double r=gr.random();
388  if (weight>0.0 && (weight>weight_old || r<weight/weight_old)) {
389  for (size_t i=0;i<npar;i++) {
390  par_old[i]=par[i];
391  }
392  weight_old=weight;
393  }
394 
395  // If we're not in the warm up phase
396  if (it>=n_warm_up) {
397 
398  for(size_t i=0;i<npar;i++) {
399  // Only update if its in range
400  if (par_old[i]>plo2[i] && par_old[i]<phi2[i]) {
401  harr[i].update(par_old[i]);
402  }
403  }
404 
405  if (it-n_warm_up==imeas) {
406 
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);
410  }
411  }
412  hist_ev.add(mat_tmp);
413 
414  // Clear the histogram data but leave the grid
415  for(size_t i=0;i<npar;i++) {
416  harr[i].clear_wgts();
417  }
418 
419  // Setup for the next measurement
420  imeas+=meas_size;
421  }
422 
423  }
424 
425  }
426 
427  ubmatrix avg(npar,hsize), std_dev(npar,hsize), avg_err(npar,hsize);
428  hist_ev.current_avg(avg,std_dev,avg_err);
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));
432  }
433  }
434 
435  return 0;
436  }
437 
438  };
439 
440 #ifndef DOXYGEN_NO_O2NS
441 }
442 #endif
443 
444 #endif
o2scl::fit_bayes::lydat
vec_t * lydat
Y-values.
Definition: fit_bayes.h:103
boost::numeric::ublas::matrix< double >
o2scl::rng_gsl
Random number generator (GSL)
Definition: rng_gsl.h:55
o2scl::fit_bayes::n_iter
size_t n_iter
Number of total iterations (default 1000)
Definition: fit_bayes.h:121
boost::numeric::ublas::vector< double >
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::fit_bayes::lyerr
vec_t * lyerr
Y-errors.
Definition: fit_bayes.h:106
o2scl::fit_bayes::lxdat
vec_t * lxdat
X-values.
Definition: fit_bayes.h:100
o2scl::fit_bayes::evidence
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.
Definition: fit_bayes.h:140
o2scl::fit_bayes::n_warm_up
size_t n_warm_up
Number of warmup iterations (default 100)
Definition: fit_bayes.h:118
o2scl::fit_bayes
Fit a function to data using Bayesian methods.
Definition: fit_bayes.h:77
o2scl::itp_linear
@ itp_linear
Linear.
Definition: interp.h:71
o2scl::fit_bayes::hsize
size_t hsize
Histogram size (default 20)
Definition: fit_bayes.h:124
o2scl::fit_bayes::pri
multi_func_t * pri
Prior distribution.
Definition: fit_bayes.h:130
o2scl::rng_gsl::random
double random()
Return a random number in .
Definition: rng_gsl.h:82
o2scl::expval_matrix::current_avg
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.
Definition: expval.h:799
o2scl::fit_bayes::integrand
virtual double integrand(size_t npar, const vec_t &par)
The integrand for the evidence.
Definition: fit_bayes.h:89
o2scl::quadratic_extremum_x
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.
Definition: misc.h:401
o2scl::fit_bayes::weight_fun
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)
Definition: fit_bayes.h:160
o2scl::mcarlo_vegas
Multidimensional integration using Vegas Monte Carlo (GSL)
Definition: mcarlo_vegas.h:115
o2scl::multi_funct
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct
Multi-dimensional function typedef in src/base/multi_funct.h.
Definition: multi_funct.h:45
o2scl::fit_bayes::fit
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.
Definition: fit_bayes.h:178
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::expval_matrix::add
void add(mat_t &val)
Add measurement of value val.
Definition: expval.h:643
o2scl::uniform_grid_end
Linear grid with fixed number of bins and fixed endpoint.
Definition: uniform_grid.h:333
o2scl::fit_bayes::fit_hist
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.
Definition: fit_bayes.h:335
o2scl::mcarlo_vegas::minteg_err
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.
Definition: mcarlo_vegas.h:845
o2scl::expval_vector
Vector expectation value.
Definition: expval.h:304
o2scl::expval_matrix
Matrix expectation value.
Definition: expval.h:584
o2scl::fit_bayes::gr
rng_gsl gr
Random number generator.
Definition: fit_bayes.h:133
o2scl::vector_find_level
void vector_find_level(double level, size_t n, vec_t &x, vec2_t &y, std::vector< double > &locs)
Perform inverse linear interpolation.
Definition: interp.h:2225
o2scl::vector_integ_interp
double vector_integ_interp(size_t n, ovec_t &v, size_t interp_type)
Integral of a vector from interpolation object.
Definition: interp.h:2120
o2scl::fit_bayes::def_inte
mcarlo_vegas def_inte
Default Monte Carlo integrator.
Definition: fit_bayes.h:136
o2scl::expval_vector::add
void add(vec_t &val)
Add measurement of value val.
Definition: expval.h:359
o2scl::fit_bayes::lndat
size_t lndat
Number of data points.
Definition: fit_bayes.h:97
o2scl::fit_bayes::ff
fit_func_t * ff
User-specified function.
Definition: fit_bayes.h:94
o2scl::uniform_prior
An unnormalized uniform prior distribution for several variables.
Definition: fit_bayes.h:50
o2scl::vector_invert_enclosed_sum
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.
Definition: interp.h:2295
o2scl::fit_bayes::nmeas
size_t nmeas
Number of measurements (default 20)
Definition: fit_bayes.h:127
o2scl::expval_vector::current_avg
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.
Definition: expval.h:484
o2scl::uniform_prior::operator()
virtual double operator()(size_t npar, const vec_t &par)
Return the value of the prior at the specified point in parameter space.
Definition: fit_bayes.h:57

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).