Probability Distributions and Markov Chain Monte Carlo

Probability distributions

Probability distributions are provided by the C++ standard library, but multidimensional distributions are not provided. For the time being, some experimental probability distributions are being included in O2scl .

The one-dimensional probability distributions are children of o2scl::prob_dens_func and are modeled similar to the C++ standard distributions. The following distributions are currently included

Multi-dimensional distributions are children of o2scl::prob_dens_mdim , including

Conditional probability distributions are children of o2scl::prob_cond_mdim , including

Markov chain Monte Carlo

The class o2scl::mcmc_para_base performs generic MCMC simulations and the child class o2scl::mcmc_para_table performs MCMC simulations and stores the results in a o2scl::table_units object. These classes contain support for OpenMP and MPI (or both). The class o2scl::mcmc_para_cli is a specialized class which automatically handles the command-line interface when using MCMC.

MCMC example

/* Example: ex_mcmc.cpp
-------------------------------------------------------------------
An example which demonstrates the generation of an arbitrary
distribution through Markov chain Monte Carlo.
*/
#include <o2scl/mcmc_para.h>
#include <o2scl/vec_stats.h>
#include <o2scl/test_mgr.h>
#include <o2scl/hdf_io.h>
using namespace std;
using namespace o2scl;
using namespace o2scl_hdf;
/// Convenient typedefs
typedef std::function<int(size_t,const ubvector &,double &,
std::array<double,1> &)> point_funct;
typedef std::function<int(const ubvector &,double,std::vector<double> &,
std::array<double,1> &)> fill_funct;
/// The MCMC object
class exc {
public:
/** \brief A one-dimensional bimodal distribution
Here, the variable 'log_weight' stores the natural logarithm of
the objective function based on the parameters stored in \c pars.
The object 'dat' stores any auxillary quantities which can be
computed at every point in parameter space.
*/
int bimodal(size_t nv, const ubvector &pars, double &log_weight,
std::array<double,1> &dat) {
double x=pars[0];
log_weight=log(exp(-x*x)*(sin(x-1.4)+1.0));
dat[0]=x*x;
return 0;
}
/** \brief Add auxillary quantities to 'line' so they can be
stored in the table
*/
int fill_line(const ubvector &pars, double log_weight,
std::vector<double> &line, std::array<double,1> &dat) {
line.push_back(dat[0]);
return 0;
}
exc() {
}
};
int main(int argc, char *argv[]) {
cout.setf(ios::scientific);
exc e;
// Parameter limits and initial points
ubvector low_bimodal(1), high_bimodal(1);
low_bimodal[0]=-5.0;
high_bimodal[0]=5.0;
// Function objects for the MCMC object
point_funct bimodal_func=std::bind
(std::mem_fn<int(size_t,const ubvector &,double &,
std::array<double,1> &)>(&exc::bimodal),&e,
std::placeholders::_1,std::placeholders::_2,std::placeholders::_3,
std::placeholders::_4);
fill_funct fill_func=std::bind
(std::mem_fn<int(const ubvector &,double,std::vector<double> &,
std::array<double,1> &)>(&exc::fill_line),&e,
std::placeholders::_1,std::placeholders::_2,std::placeholders::_3,
std::placeholders::_4);
vector<point_funct> bimodal_vec;
bimodal_vec.push_back(bimodal_func);
vector<fill_funct> fill_vec;
fill_vec.push_back(fill_func);
cout << "----------------------------------------------------------"
<< endl;
cout << "Plain MCMC example with a bimodal distribution:" << endl;
// Set parameter names and units
vector<string> pnames={"x","x2"};
vector<string> punits={"",""};
mct.set_names_units(pnames,punits);
// Set MCMC parameters
mct.step_fac=2.0;
mct.max_iters=20000;
mct.prefix="ex_mcmc";
// Perform MCMC
mct.mcmc(1,low_bimodal,high_bimodal,bimodal_vec,fill_vec);
// Output acceptance and rejection rate
cout << "n_accept, n_reject: " << mct.n_accept[0] << " "
<< mct.n_reject[0] << endl;
// Get results
shared_ptr<table_units<> > t=mct.get_table();
// Remove empty rows from table
t->delete_rows_func("mult<0.5");
// Compute autocorrelation length and effective sample size
std::vector<double> ac, ftom;
(*t)["x2"],(*t)["mult"],ac);
size_t ac_len=o2scl::vector_autocorr_tau(ac,ftom);
cout << "Autocorrelation length, effective sample size: "
<< ac_len << " " << t->get_nlines()/ac_len << endl;
// Create a set of fully independent samples
std::vector<double> indep;
size_t count=0;
for(size_t j=0;j<t->get_nlines();j++) {
for(size_t k=0;k<((size_t)(t->get("mult",j)+1.0e-8));k++) {
if (count==ac_len) {
indep.push_back(t->get("x2",j));
count=0;
}
count++;
}
}
// Use the independent samples to compute the final integral and
// compare to the exact result
double avg=vector_mean(indep);
double std=vector_stddev(indep);
tm.test_rel(avg,1.32513,10.0*std/sqrt(indep.size()),"ex_mcmc");
cout << avg << " " << 1.32513 << " " << fabs(avg-1.32513) << " "
<< std << " " << 10.0*std/sqrt(indep.size()) << endl;
tm.report();
return 0;
}
// End of example
boost::numeric::ublas::matrix< double >
o2scl::mcmc_para_base< func_t, std::function< int(const ubvector &, double, size_t, int, bool, data_t &)>, data_t, ubvector >::n_accept
std::vector< size_t > n_accept
The number of Metropolis steps which were accepted in each independent chain (summed over all walkers...
Definition: mcmc_para.h:269
boost::numeric::ublas::vector< double >
o2scl::mcmc_para_base< func_t, std::function< int(const ubvector &, double, size_t, int, bool, data_t &)>, data_t, ubvector >::step_fac
double step_fac
Stepsize factor (default 10.0)
Definition: mcmc_para.h:316
o2scl::mcmc_para_base< func_t, std::function< int(const ubvector &, double, size_t, int, bool, data_t &)>, data_t, ubvector >::prefix
std::string prefix
Prefix for output filenames (default "mcmc")
Definition: mcmc_para.h:310
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::mcmc_para_table::set_names_units
virtual void set_names_units(std::vector< std::string > names, std::vector< std::string > units)
Set the table names and units.
Definition: mcmc_para.h:2167
o2scl::mcmc_para_table::mcmc
virtual int mcmc(size_t n_params_local, vec_t &low, vec_t &high, std::vector< func_t > &func, std::vector< fill_t > &fill)
Perform an MCMC simulation.
Definition: mcmc_para.h:2468
o2scl::test_mgr::report
bool report() const
Provide a report of all tests so far.
o2scl::vector_mean
double vector_mean(size_t n, const vec_t &data)
Compute the mean of the first n elements of a vector.
Definition: vec_stats.h:55
o2scl::mcmc_para_base< func_t, std::function< int(const ubvector &, double, size_t, int, bool, data_t &)>, data_t, ubvector >::max_iters
size_t max_iters
If non-zero, the maximum number of MCMC iterations (default 0)
Definition: mcmc_para.h:297
o2scl::mcmc_para_table::table_prealloc
size_t table_prealloc
Number of rows to allocate for the table before the MCMC run.
Definition: mcmc_para.h:2026
o2scl::vector_stddev
double vector_stddev(size_t n, const vec_t &data)
Standard deviation with specified mean.
Definition: vec_stats.h:253
o2scl::test_mgr::test_rel
bool test_rel(data_t result, data_t expected, data_t rel_error, std::string description)
Test for .
Definition: test_mgr.h:168
o2scl::test_mgr
A class to manage testing and record success and failure.
Definition: test_mgr.h:50
o2scl_hdf
The O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl namespace ...
Definition: table.h:59
o2scl::mcmc_para_table::get_table
std::shared_ptr< o2scl::table_units<> > get_table()
Get the output table.
Definition: mcmc_para.h:2504
o2scl::vector_autocorr_vector_mult
void vector_autocorr_vector_mult(size_t n2, const vec_t &data, const vec2_t &mult, resize_vec_t &ac_vec)
Construct an autocorrelation vector using a multiplier using the first n2 elements of vectors data an...
Definition: vec_stats.h:1962
o2scl::test_mgr::set_output_level
void set_output_level(int l)
Set the output level.
Definition: test_mgr.h:119
o2scl::vector_autocorr_tau
size_t vector_autocorr_tau(const vec_t &ac_vec, resize_vec_t &five_tau_over_M)
Use the Goodman method to compute the autocorrelation length.
Definition: vec_stats.h:1857
o2scl::mcmc_para_base< func_t, std::function< int(const ubvector &, double, size_t, int, bool, data_t &)>, data_t, ubvector >::n_reject
std::vector< size_t > n_reject
The number of Metropolis steps which were rejected in each independent chain (summed over all walkers...
Definition: mcmc_para.h:276
o2scl::mcmc_para_table
A generic MCMC simulation class writing data to a o2scl::table_units object.
Definition: mcmc_para.h:1808

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