Monte Carlo Integration

Monte Carlo integration is performed by descendants of mcarlo_inte (o2scl::mcarlo_plain, o2scl::mcarlo_miser, and o2scl::mcarlo_vegas). These routines are generally superior to the direct methods for integrals over regions with large numbers of spatial dimensions.

This example computes the integral

and compares the result with the exact value, 1.3932039296.

/* Example: ex_mcarlo.cpp

-------------------------------------------------------------------

An example to demonstrate multidimensional integration

*/

#include <boost/numeric/ublas/vector.hpp>

#include <o2scl/test_mgr.h>

#include <o2scl/multi_funct.h>

#include <o2scl/inte_qng_gsl.h>

#include <o2scl/mcarlo_vegas.h>

/// For M_PI

#include <gsl/gsl_math.h>

using namespace std;

using namespace o2scl;

typedef boost::numeric::ublas::vector<double> ubvector;

double test_fun(size_t nv, const ubvector &x) {

double y=1.0/(1.0-cos(x[0])*cos(x[1])*cos(x[2]))/M_PI/M_PI/M_PI;

return y;

}

int main(void) {

test_mgr t;

t.set_output_level(1);

cout.setf(ios::scientific);

double exact=1.3932039296;

double res;

double err;

mcarlo_vegas<> gm;

ubvector a(3), b(3);

a[0]=0.0; b[0]=M_PI;

a[1]=0.0; b[1]=M_PI;

a[2]=0.0; b[2]=M_PI;

multi_funct11 tf=test_fun;

gm.n_points=100000;

gm.minteg_err(tf,3,a,b,res,err);

cout << res << " " << exact << " " << (res-exact)/err << endl;

t.test_rel(res,exact,err*10.0,"O2scl");

t.report();

return 0;

}

// End of example

The base o2scl::expval_base and its children form a set of classes useful for recording the outputs of several iterations of a numerical simulation, and summarizing the average, standard deviation, and error in the average over all iterations. After each iteration, some measurement is performed, and that measurement is added to the class with an `add()`

functions (e.g. o2scl::expval_scalar::add() ).

Autocorrelations are common in numerical simulations, and these classes allow the user to assess the impact of correlations by breaking the measurements up into blocks.

Blocks are filled one by one, moving to the next block when the requested of measurements (`n_per_block`

) in the previous block has been provided. If more measurements are given than the number of blocks times `n_per_block`

, then the number of filled blocks is halved, each block is filled with twice as much data, and `n_per_block`

is doubled. The blocks are always combined in such a way as to preserve the original ordering, i.e. the first block is combined with the second, the third with the fourth, and so on. Subsequent measurements are added to newly empty blocks.

The constructor for these classes need as input the number of blocks to record (at least one) and the number of measurements per block (also at least one). If either of these is specified to be zero then the error handler is called.

In order to properly handle autocorrelations, The number of measurements over which autocorrelation occurs should be less than the number of measurements per block.

Current averages are always reported as long as at least one measurement has been added, though not all data is always used if there are incomplete blocks or the blocks have not yet been filled evenly. Two functions `current_avg()`

and `current_avg_stats()`

are provided in children. The latter provides the number of blocks and number of measurements per block used in the currently reported statistics. If only one (full or partially full) block is available, the standard deviation and uncertainty on the average are reported as zero. If no data is available, then calls to `current_avg()`

will call the error handler. See, e.g. o2scl::expval_scalar::current_avg() and o2scl::expval_scalar::current_avg_stats().

Many of the o2scl::expval_base children also support HDF I/O.

Note that the traditional Monte Carlo integration routines (o2scl::mcarlo_plain, o2scl::mcarlo_miser, and o2scl::mcarlo_vegas) have no significant autocorrelations by construction, so long as the random number generator does not have any.

For Markov chain Monte Carlo, see the discussion in Markov chain Monte Carlo.

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