fit_bayes.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2011-2017, 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_funct11, 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_funct11 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_linear(hsize+2,edge_x,edge_y);
275 
276  double lev;
277  std::vector<double> locs;
278  vector_invert_enclosed_sum(0.68*total,hsize+2,edge_x,edge_y,lev);
279  vector_find_level(lev,hsize+2,edge_x,edge_y,locs);
280 
281  double lmin=*min_element(locs.begin(),locs.end());
282  double lmax=*max_element(locs.begin(),locs.end());
283  xlo.push_back(lmin);
284  xhi.push_back(lmax);
285 
286  size_t imax=vector_max_index<std::vector<double>,double>
287  (hsize+2,edge_y);
288  xmax.push_back(quadratic_extremum_x
289  (edge_x[imax-1],edge_x[imax],edge_x[imax+1],
290  edge_y[imax-1],edge_y[imax],edge_y[imax+1]));
291  }
292 
293  // Add the limits to the expval_vector object
294  low.add(xlo);
295  high.add(xhi);
296  max.add(xmax);
297 
298  // Clear the histogram data but leave the grid
299  for(size_t i=0;i<npar;i++) {
300  harr[i].clear_wgts();
301  }
302 
303  // Setup for the next measurement
304  imeas+=meas_size;
305  }
306 
307  }
308 
309  }
310 
311  // Now get parameter values and errors from the expval_vector objects
312  ubvector sd1(npar), sd2(npar), sd3(npar);
313  low.current_avg(plo2,sd1,plo_err);
314  high.current_avg(phi2,sd2,phi_err);
315  max.current_avg(pmax,sd3,pmax_err);
316 
317  for(size_t i=0;i<npar;i++) {
318  if (((double)nonzero_bins[i])<((double)nmeas)*2.5) {
319  O2SCL_ERR2("Not enough data for accurate parameter ",
320  "estimation in fit_bayes::fit().",exc_efailed);
321  }
322  }
323 
324  return 0;
325  }
326 
327  /** \brief Desc
328 
329  For each measurement block, this function collects the data for
330  all the parameters into 1d histogram objects. Then, at the end
331  of the block, the histogram information is added to a hist
332  object for each parameter.
333  */
334  virtual int fit_hist(size_t ndat, vec_t &xdat, vec_t &ydat, vec_t &yerr,
335  size_t npar, vec_t &plo2, vec_t &phi2,
336  std::vector<hist> &par_hist, fit_func_t &fitfun,
337  multi_func_t &prior_fun) {
338 
339  ff=&fitfun;
340  pri=&prior_fun;
341 
342  double weight, weight_old;
343 
344  // Create par_old and set up initial weight in weight_old
345  vec_t par_old, par;
346  par_old.resize(npar);
347  par.resize(npar);
348 
349  // Choose the first point as the midpoint of the parameter space
350  for(size_t i=0;i<npar;i++) {
351  par_old[i]=(plo2[i]+phi2[i])/2.0;
352  }
353  weight_old=weight_fun(ndat,xdat,ydat,yerr,npar,par_old)*
354  (*pri)(npar,par_old);
355  if (weight_old==0.0) {
356  O2SCL_ERR2("Initial weight zero in ",
357  "fit_bayes::fit().",exc_einval);
358  }
359 
360  // Set up the histograms for marginal estimation
361  par_hist.resize(npar);
362  std::vector<hist> harr(npar);
363  for(size_t i=0;i<npar;i++) {
364  harr[i].set_bin_edges(uniform_grid_end<>(plo2[i],phi2[i],hsize));
365  harr[i].clear_wgts();
366  par_hist[i].set_bin_edges(uniform_grid_end<>(plo2[i],phi2[i],hsize));
367  par_hist[i].clear_wgts();
368  }
369 
370  // Create ev objects for lower and upper limits of each parameter
371  size_t meas_size=((size_t)(((double)n_iter)/((double)nmeas)));
372  size_t imeas=meas_size-1;
373  expval_matrix hist_ev(npar,hsize,nmeas,1);
374  ubmatrix mat_tmp(npar,hsize);
375 
376  // Main iteration
377  for(size_t it=0;it<n_iter+n_warm_up;it++) {
378 
379  // Select new point in parameter space and evaluate it's weight
380  for (size_t i=0;i<npar;i++) {
381  par[i]=gr.random()*(phi2[i]-plo2[i])+plo2[i];
382  }
383  weight=weight_fun(ndat,xdat,ydat,yerr,npar,par)*(*pri)(npar,par);
384 
385  // Metropolis sampling
386  double r=gr.random();
387  if (weight>0.0 && (weight>weight_old || r<weight/weight_old)) {
388  for (size_t i=0;i<npar;i++) {
389  par_old[i]=par[i];
390  }
391  weight_old=weight;
392  }
393 
394  // If we're not in the warm up phase
395  if (it>=n_warm_up) {
396 
397  for(size_t i=0;i<npar;i++) {
398  // Only update if its in range
399  if (par_old[i]>plo2[i] && par_old[i]<phi2[i]) {
400  harr[i].update(par_old[i]);
401  }
402  }
403 
404  if (it-n_warm_up==imeas) {
405 
406  for(size_t i=0;i<npar;i++) {
407  for(size_t j=0;j<hsize;j++) {
408  mat_tmp(i,j)=harr[i].get_wgt_i(j);
409  }
410  }
411  hist_ev.add(mat_tmp);
412 
413  // Clear the histogram data but leave the grid
414  for(size_t i=0;i<npar;i++) {
415  harr[i].clear_wgts();
416  }
417 
418  // Setup for the next measurement
419  imeas+=meas_size;
420  }
421 
422  }
423 
424  }
425 
426  ubmatrix avg(npar,hsize), std_dev(npar,hsize), avg_err(npar,hsize);
427  hist_ev.current_avg(avg,std_dev,avg_err);
428  for(size_t i=0;i<npar;i++) {
429  for(size_t j=0;j<hsize;j++) {
430  par_hist[i].set_wgt_i(j,avg(i,j));
431  }
432  }
433 
434  return 0;
435  }
436 
437  };
438 
439 #ifndef DOXYGEN_NO_O2NS
440 }
441 #endif
442 
443 #endif
multi_func_t * pri
Prior distribution.
Definition: fit_bayes.h:130
Matrix expectation value.
Definition: expval.h:579
size_t hsize
Histogram size (default 20)
Definition: fit_bayes.h:124
vec_t * lydat
Y-values.
Definition: fit_bayes.h:103
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
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:479
vec_t * lyerr
Y-errors.
Definition: fit_bayes.h:106
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
invalid argument supplied by user
Definition: err_hnd.h:59
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
double random()
Return a random number in .
Definition: rng_gsl.h:82
virtual double integrand(size_t npar, const vec_t &par)
The integrand for the evidence.
Definition: fit_bayes.h:89
An unnormalized uniform prior distribution for several variables.
Definition: fit_bayes.h:50
size_t lndat
Number of data points.
Definition: fit_bayes.h:97
rng_gsl gr
Random number generator.
Definition: fit_bayes.h:133
fit_func_t * ff
User-specified function.
Definition: fit_bayes.h:94
generic failure
Definition: err_hnd.h:61
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
mcarlo_vegas def_inte
Default Monte Carlo integrator.
Definition: fit_bayes.h:136
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:323
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:789
size_t nmeas
Number of measurements (default 20)
Definition: fit_bayes.h:127
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:334
Vector expectation value.
Definition: expval.h:304
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
size_t n_warm_up
Number of warmup iterations (default 100)
Definition: fit_bayes.h:118
double vector_integ_linear(size_t n, vec_t &x, vec2_t &y)
Compute the integral over y(x) using linear interpolation.
Definition: interp.h:2011
vec_t * lxdat
X-values.
Definition: fit_bayes.h:100
size_t n_iter
Number of total iterations (default 1000)
Definition: fit_bayes.h:121
void vector_invert_enclosed_sum(double sum, size_t n, vec_t &x, vec2_t &y, double &lev, int verbose=0)
Compute the endpoints which enclose the regions whose integral is equal to sum.
Definition: interp.h:2058
Multidimensional integration using Vegas Monte Carlo (GSL)
Definition: mcarlo_vegas.h:115
Random number generator (GSL)
Definition: rng_gsl.h:55
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
Fit a function to data using Bayesian methods.
Definition: fit_bayes.h:77
void add(vec_t &val)
Add measurement of value val.
Definition: expval.h:354
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct11
Multi-dimensional function typedef.
Definition: multi_funct.h:45
void add(mat_t &val)
Add measurement of value val.
Definition: expval.h:633
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:1965
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
Linear grid with fixed number of bins and fixed endpoint.
Definition: uniform_grid.h:315

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