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 .
#include <o2scl/mcmc_para.h>
#include <o2scl/vec_stats.h>
#include <o2scl/test_mgr.h>
#include <o2scl/hdf_io.h>
using namespace std;
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;
class exc {
public:
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;
}
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;
ubvector low_bimodal(1), high_bimodal(1);
low_bimodal[0]=-5.0;
high_bimodal[0]=5.0;
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;
vector<string> pnames={"x","x2"};
vector<string> punits={"",""};
mct.
mcmc(1,low_bimodal,high_bimodal,bimodal_vec,fill_vec);
cout <<
"n_accept, n_reject: " << mct.
n_accept[0] <<
" "
shared_ptr<table_units<> > t=mct.
get_table();
t->delete_rows_func("mult<0.5");
std::vector<double> ac, ftom;
(*t)["x2"],(*t)["mult"],ac);
cout << "Autocorrelation length, effective sample size: "
<< ac_len << " " << t->get_nlines()/ac_len << endl;
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++;
}
}
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;
return 0;
}