# Monte Carlo Integration¶

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

and compares the result with the exact value, 1.3932039296.

```
/*
-------------------------------------------------------------------
Copyright (C) 2006-2022, Andrew W. Steiner
This file is part of O2scl.
O2scl is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
O2scl is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with O2scl. If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------
*/
/* 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_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.