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.

# Monte Carlo integration example

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;
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) {
cout.setf(ios::scientific);
double exact=1.3932039296;
double res;
double err;
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

# Analysis of results from numerical simulations

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).