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

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