Monte Carlo Integration

O2scl

Monte Carlo contents

Monte Carlo introduction

Monte Carlo integration is performed by descendants of mcarlo (mcarlo_plain, mcarlo_miser, and 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

\[\int_0^{\pi} \int_0^{\pi} \int_0^{\pi} \frac{1}{\pi^3} \left(1-\cos x \cos y \cos z\right)^{-1}\]

and compares the result with the exact value, 1.3932039296.

/* Example: ex_mcarlo.cpp
   -------------------------------------------------------------------
   An example to demonstrate multidimensional integration. See "License 
   Information" section of the documentation for license information.
*/

#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_funct 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

Note

Jan 2021: These classes are not well-developed and may be superceded by the new GSL running stats functions.

The base 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. 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 expval_base children also support HDF I/O.

Note that the traditional Monte Carlo integration routines (mcarlo_plain, mcarlo_miser, and 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 Probability Distributions and Markov Chain Monte Carlo.