Struct MISER 

pub struct MISER<F, V, T> { /* private fields */ }

The MISER Algorithm

The Monte Carlo integration algorithm MISER by Press and Farrar (See [1]) uses stratified sampling based on variance of the integrand’s value on the sampled points. More specifically, the original subregion (which must be a Cartesian product of intervals [a, b] for real a, b) is bisected along some dimension such that the integrand’s value in those two subregions gives the smallest total variance amongst all other bisections. The total samples points are then distributed to those subregions according to the variance. This procedure is continued recursively until either the recursion depth is reached or the number of sample points allocated to that subregion fall under an (adjustable) threshold. The procedure is implemented in (private function) miser_recurse.

There are several adjustable parameters in this implementation of MISER. The parameters are inspired by those used in [1]. The parameters in MISER are:

• Sample Count: (Default 1000) A lower bound on the number of points to sample in total. The reason this is a lower bound is because if the minimum threshold is fallen under then the minimum will be samples. For example, if a region is allocated 5 points and the (total) threshold value is 32, then 32 values will be sampled. Thus, more than then sample count may be allocated in total. For more details, read Minimum Dimensional Sample Count below or read the docstrings of MISER::with_sample_count and MISER::with_minimum_dimensional_sample_count.
• Maximum Depth: (Default 5) The maximum recursion depth, i.e. the recursive calls will be no deeper than the number given by this parameter.
• Minimum Dimensional Sample Count: (Default 32) The minimum number of points that will be sampled for a subregion, per dimension. For example if this number is 32 and the input is 4 dimensional (in R^4) then, this number of 4x32 = 128 and if a subregion is ever allocated less than 128 points, then the normal Monte Carlo integration algorithm will be used on that subregion with 128 samples points. Please also read the docstring of MISER::with_minimum_dimensional_sample_count.
• Alpha (Default 2.): When a bisection is chosen the number of allocated points for each of the two subregions is chosen based on the variance. The exact formula is for the first region is N_1 = N * (Var_1^beta / (Var_1^beta + Var_2^beta)) where N is the total sample points, Var_1 and Var_2 are the variances of the integrand on the first and second subregions, respectively, and beta = 1 / (1 + alpha). Simiarly, for the second subregion the formula is: N_2 = N * (Var_2^beta / (Var_1^beta + Var_2^beta)). In (1), the authors recommended alpha being 2. and hence this is the default used here. See also the docstring of MISER::with_alpha.
• Dither (Default 0.05): Instead of bisecting exactly in on the midpoint of the boundries of the chosen dimension, a small dither is added to break the symmetry of certain functions. Thus, instead of having the regions [a, 0.5*(a + b)] and [0.5*(a + b), b], a random value, call it r, in the range [-dither, dither] is chosen and the two subregions become: [a, (0.5 + r)*(a + b)] and [(0.5 + r)*(a + b), b]. NOTE, therefore that the dither should be in the interval (0, 0.5). See also the docstring of MISER::with_dither.

Examples

MISER implements the OrthotopeRandomIntegrator and hence we can use it thusly:

use eqsolver::integrators::{MISER, OrthotopeRandomIntegrator};
use nalgebra::{Vector2, vector};
use rand_chacha::{ChaCha8Rng, rand_core::SeedableRng}; // or any rng you like
let mut rng = ChaCha8Rng::seed_from_u64(1729);
let f = |v: Vector2<f64>| (v[1].cos() + v[0]).sin();
let lower = vector![0., 0.];
let upper = vector![1., 1.];
let result = MISER::new(f).integrate_with_rng(lower, upper, &mut rng);
assert!((result.unwrap().mean - 0.924846).abs() <= 0.001);

References:

(1) Press, W.H. and Farrar, G.R., 1990. Recursive stratified sampling for multidimensional Monte Carlo integration. Computers in Physics, 4(2), pp.190-195. (2) https://www.gnu.org/software/gsl/doc/html/montecarlo.html

Implementations

impl<F, V, T> MISER<F, V, T>

where T: Float, F: Fn(V) -> T,

pub fn new(f: F) -> Self

Creates a new instance of the algorithm

Instantiates the MISER integrator using a function f: R -> R

pub fn with_sample_count(&mut self, sample_count: usize) -> &mut Self

Set the total number of points that will be samples in the algorithm.

NOTE! This number if a lower bound on the number of sampled points. This is because when a bisection leads to a subregion having less than the minimum sample count allocated, then the minimum sample count is sampled instead. See MISER::with_minimum_dimensional_sample_count.

Default sample count: 1000

Examples

let result = MISER::new(f)
                .with_sample_count(2000)
                .integrate_with_rng(lower, upper, &mut rng);

pub fn with_minimum_dimensional_sample_count( &mut self, minimum_dimensional_sample_count: usize, ) -> &mut Self

Set the minimum number of points per dimension that will be sampled in a region.

In other words, if a subregion has less than this number times the number of dimensions then the recursion stops and the plain Monte Carlo algorithm is ran on that subregion using the number of samples points specified by this function times the number of dimensions.

Default minimum sample count per dimension: 32

Examples

let result = MISER::new(f)
                .with_minimum_dimensional_sample_count(50)
                .integrate_with_rng(lower, upper, &mut rng);

pub fn with_maximum_depth(&mut self, maximum_depth: usize) -> &mut Self

Set the maximum recursion depth of the algorithm.

In other words, the algorithm’s recursion tree is no deeper than the number specified using this function. The recursion in this algorithm happens when a region (orthotope) is bisected along an axis. When the maximum recursion depth is reached, the plain Monte Carlo algorithm is used on that subregion (which was one or more bisections of a larger region).

Default maximum depth: 5

let result = MISER::new(f)
                .with_maximum_depth(8)
                .integrate_with_rng(lower, upper, &mut rng);

pub fn with_alpha(&mut self, alpha: T) -> &mut Self

Set the alpha value of the algorithm.

The alpha value changes the proportions of the allocated samples points to each subregion of the bisection that happens in every recursion call of MISER. More precisely, the allocated points for each subregions follow the following formulas:

• N_1 = N * (Var_1^beta / (Var_1^beta + Var_2^beta))
• N_2 = N * (Var_2^beta / (Var_1^beta + Var_2^beta))

where N is the total sample points, Var_1 and Var_2 are the variances of the integrand on the first and second subregions, respectively, and beta = 1 / (1 + alpha).

For more details read the docstring of MISER.

Default alpha: 2.

let result = MISER::new(f)
                .with_alpha(1.5)
                .integrate_with_rng(lower, upper, &mut rng);

pub fn with_dither(&mut self, dither: T) -> &mut Self

Set the dither value of the algorithm.

The dither value controls value that makes the bisection to not be cut exactly at the midpoint. This value helps break certain symmetries of functions. For example, if a region is about to be cut along an axis whose bounds are a < b then, instead of having the regions [a, 0.5*(a + b)] and [0.5*(a + b), b], a random value, call it r, in the range [-dither, dither] is chosen and the two subregions become:

• [a, (0.5 + r)*(a + b)], and
• [(0.5 + r)*(a + b), b].

NOTE! The dither must be in the range (0, 0.5) (this is a precondition).

Default dither: 0.05

let result = MISER::new(f)
                .with_dither(0.01)
                .integrate_with_rng(lower, upper, &mut rng);

Trait Implementations

impl<F, T, D> OrthotopeRandomIntegrator<Matrix<T, D, Const<1>, <DefaultAllocator as Allocator<D>>::Buffer<T>>, T> for MISER<F, Vector<T, D, <DefaultAllocator as Allocator<D>>::Buffer<T>>, T>

where T: Float + Scalar + ComplexField<RealField = T> + SampleUniform, D: Dim, F: Fn(Vector<T, D, <DefaultAllocator as Allocator<D>>::Buffer<T>>) -> T, DefaultAllocator: Allocator<D>,

fn integrate_with_rng( &self, from: Vector<T, D, <DefaultAllocator as Allocator<D>>::Buffer<T>>, to: Vector<T, D, <DefaultAllocator as Allocator<D>>::Buffer<T>>, rng: &mut impl Rng, ) -> SolverResult<MeanVariance<T>>

Run the integrator over an interval or a Cartesian product of intervals. Represented as either a Float or a nalgebra::Vector. Since the integrator is random, a random number generator is also an input to this function of type rng::Rng.

fn integrate_to_mean_variance( &self, from: V, to: V, ) -> SolverResult<MeanVariance<T>>

Run the integrator over an interval or a Cartesian product of intervals using the rand::rng() as the underlying random number generator. The output of this function is the mean and variance of the output of the randomised integration algorithm represented as a MeanVariance.

fn integrate(&self, from: V, to: V) -> SolverResult<T>

Run the integrator over an interval or a Cartesian product of intervals using rand::rng() as the underlying random number generator and return only the mean.