Skip to contents

Integrate a user-written function over a multidimensional box

Source: R/callback.R
quadrature_nd.Rd

Computes the definite integral of f over the hyper-rectangle [min, max] (P2 "math extras") with one of three ported stochastic multidimensional integrators: plain Monte Carlo (method = "monte_carlo", the default), Miser (recursive stratified-sampling Monte Carlo, Press et al. "Numerical Recipes" Sec. 7.9), or Vegas (Lepage's adaptive importance sampling, Sec. 7.8, with an optional Power Transform for rare tail-event sampling). min/max give both the per-dimension bounds and, via their length, the number of dimensions.

Usage

quadrature_nd(
  f,
  min,
  max,
  method = c("monte_carlo", "miser", "vegas"),
  seed = NULL,
  use_sobol = TRUE,
  max_function_evaluations = NULL,
  min_iterations = NULL,
  max_iterations = NULL,
  relative_tolerance = NULL,
  fraction = NULL,
  min_subregion_points = NULL,
  min_bisections = NULL,
  dither = NULL,
  independent_evaluations = NULL,
  function_calls = NULL,
  alpha = NULL,
  number_of_bins = NULL,
  tail_focus_parameter = NULL,
  initialize = NULL,
  check_convergence = NULL,
  target_probability = NULL
)

Arguments

f

a function taking a numeric vector and returning one number (method = "monte_carlo"/"miser"), or a function taking a numeric vector and a number (the sample weight) and returning one number (method = "vegas").

min, max

numeric vectors of the same length giving the per-dimension lower and upper bounds; their common length is the number of dimensions. Every max entry must be above the matching min entry.

method

one of "monte_carlo" (the default), "miser", or "vegas".

seed

an integer seed for the class's random number generator. NULL, the default, leaves the ported class's own clock-seeded default in force – so a NULL "miser"/"vegas" run under Sobol sampling is still reproducible (the Sobol sequence is deterministic), but a NULL "monte_carlo" run, or a "miser"/"vegas" run with use_sobol = FALSE, is not. See the Details.

use_sobol

whether to draw sample points from a Sobol low-discrepancy sequence rather than the (possibly seeded) generator. Default TRUE. Only "miser" and "vegas" read this; see the Details on why "monte_carlo" does not.

max_function_evaluations

the cap on evaluations of f. NULL, the default, leaves the ported class's own default in force. Applies to "monte_carlo" and "miser" alone – "miser"'s own recursion is bounded directly by it, where "monte_carlo"'s loop is bounded by max_iterations instead (see below); supplying it for method = "vegas" raises an error.

min_iterations, max_iterations, relative_tolerance

method = "monte_carlo" alone: the floor on iterations before the convergence check is consulted, the ceiling on iterations ("monte_carlo"'s real throttle, since max_function_evaluations is checked only after the loop ends to choose the reported status), and the relative-error convergence threshold. NULL, the default, leaves the ported class's own defaults in force. Supplying any for another method raises an error.

fraction, min_subregion_points, min_bisections, dither

method = "miser" alone: the fraction of remaining evaluations spent exploring variance at each stage, the minimum points per terminal subregion, the minimum evaluations before a subregion is bisected further, and the dither applied when the integrand's active region falls on a subdivision boundary. NULL, the default, leaves the ported class's own defaults in force. Supplying any for another method raises an error.

independent_evaluations, function_calls, alpha, number_of_bins, tail_focus_parameter, initialize, check_convergence, target_probability

method = "vegas" alone. independent_evaluations and function_calls bound the run (their product is the maximum total evaluations); alpha is the grid-refinement damping exponent; number_of_bins the stratification bin count; tail_focus_parameter the Power Transform exponent (1.0, the default, is standard uniform sampling); initialize selects a cold start (0, the default), inheriting the grid alone (1), or inheriting the grid and its answers (2); check_convergence whether to exit early on convergence. target_probability, if supplied, calls the ported configure_for_rare_events() helper – applied AFTER every other option, so it may override number_of_bins/alpha/tail_focus_parameter, exactly as the C# helper does. NULL, the default, leaves the ported class's own defaults in force. Supplying any for another method raises an error.

Value

the integral, a single number, carrying status ("Success", "MaximumFunctionEvaluationsReached", "Failure", or "None"), function_evaluations, and standard_error as attributes – and, for method = "vegas" alone, chi_squared, the Chi-Squared statistic (an approximate diagnostic for the run's own internal consistency across its independent evaluations). An upstream quirk, verified against the real C# source rather than assumed: method = "miser" always reports status "None" on success – unlike "monte_carlo" and "vegas", the ported Miser::integrate() (faithfully mirroring C#'s Miser.Integrate()) never assigns a success status, only ever writing "Failure" from its catch block.

Details

method = "vegas" takes f(x, weight) – x the sample point and weight the importance weight Vegas has already computed for it – rather than f(x), matching the upstream C# Vegas constructor's own integrand shape; a weight-ignoring wrapper (function(x, w) g(x)) reproduces an f(x)-only integrand under Vegas, exactly as the upstream unit tests wrap theirs.

The SAMPLE STREAM – which points f is called at – reproduces bit-for-bit against the same run in corehydropy, EXCEPT that seed has no effect on method = "monte_carlo" or a Sobol-sampled run of "miser"/"vegas" (the default, use_sobol = TRUE): Miser and Vegas draw their sample points from a Sobol low-discrepancy sequence rather than the Mersenne Twister seed seeds, and MonteCarloIntegration's own UseSobolSequence flag is a documented DEAD property upstream – declared but never consulted by Integrate() – so a "monte_carlo" run always draws from the generator seed seeds. use_sobol = FALSE reroutes Miser/Vegas through that same seeded generator instead. An HONEST LIMIT, measured rather than assumed: the AGGREGATED numbers this function returns are not always bit-identical between R and Python the way the sample stream is. MonteCarloIntegration/Miser/Vegas are ported CORE code, so – unlike the callback surface's own catalog tests, which the C++ side compiles with -ffp-contract=off for exactly this reason – they compile with whatever fused-multiply-add behavior each package's own build flags happen to produce (R's -O2 and corehydropy's CMake default are not guaranteed to agree), and the picture is different per method rather than uniform across the surface. method = "monte_carlo" is the one case measured to reproduce integral bit-for-bit across ALL FOUR runners – this package, corehydropy, the C++ fixture runner under both FMA settings, and the real C# library – because its own arithmetic (a running sum of hit/miss weights divided by the sample count) has no near-cancelling subtraction for a fused-multiply-add ULP to hide in; fixtures/callback/callback_cross_language.json pins it at zero tolerance for exactly that reason. "miser" and "vegas" do not reproduce that cleanly – measured directly, "miser"'s own integral misses the C# value by 1 ULP under this package's shipped build, and "vegas"'s integral, while itself exact against C#, sits beside a standard_error (both methods) and chi_squared ("vegas" only) that are not: both are built from a near-cancelling subtraction (avg2 - avg*avg-shaped for Monte Carlo/Miser, sum_chi_squared - sum_weighted_results * result for Vegas) that amplifies a fused-multiply-add ULP difference, and a live R-vs-corehydropy comparison of "vegas" at these settings showed integral itself, not just standard_error/chi_squared, moving by a few ULP language to language. This is a property of the classes' own arithmetic, not a bug in either binding, and it is why fixtures/callback/callback_cross_language.json's own quadrature_nd/quadrature_vegas digest asserts function_evaluations and status on every method, integral ADDITIONALLY on "monte_carlo" alone, and nothing else – see that file's reference note for the measurements.

Examples

# max_iterations/relative_tolerance keep this example under the CRAN 5s example budget;
# omitting them lets `"monte_carlo"`'s own (much larger) default run to full convergence.
quadrature_nd(function(x) if (x[1]^2 + x[2]^2 < 1) 1 else 0, min = c(-1, -1), max = c(1, 1),
              seed = 12345, max_iterations = 2000, relative_tolerance = 0.01)
#> [1] 3.206
#> attr(,"status")
#> [1] "Success"
#> attr(,"function_evaluations")
#> [1] 2000
#> attr(,"standard_error")
#> [1] 0.03567607
quadrature_nd(function(x, w) if (x[1]^2 + x[2]^2 < 1) 1 else 0, min = c(-1, -1), max = c(1, 1),
              method = "vegas")
#> [1] 3.141801
#> attr(,"status")
#> [1] "Success"
#> attr(,"function_evaluations")
#> [1] 19600
#> attr(,"standard_error")
#> [1] 0.002125444
#> attr(,"chi_squared")
#> [1] 0.1423939