Markov-chain based Vegas algorithm

function montecarlo(config::Configuration{N,V,P,O,T}, integrand::Function, neval,
    verbose=0, debug=false;
    measurefreq::Int=1, measure::Union{Nothing,Function}=nothing) where {N,V,P,O,T}

This function applies a Markov-chain Monte Carlo (MCMC) technique combined with the Vegas algorithm to compute integrals. In addition to calculating the original integrals, it also introduces a normalization integral with an integrand ~ 1, which enhances the efficiency and robustness of high-dimensional integration tasks.

Overview

Given multiple integrands involving multiple variables, the algorithm finds the best distribution ansatz that fits all integrands together. For instance, consider we want to calculate two integrals: $f_1(x)$ and $f_2(x, y)$, where $x$ and $y$ are two different types of variables. The algorithm learns distributions $\rho_x(x)$ and $\rho_y(y)$ such that $f_1(x)/\rho_x(x)$ and $f_2(x, y)/\rho_x(x)/\rho_y(y)$ are as flat as possible.

Then, it samples variables $x$ and $y$ using the Metropolis-Hastings algorithm with a joint distribution p(x, y),

\[p(x, y) = r_0 \cdot \rho_x(x) \cdot \rho_y(y) + r_1 \cdot |f_1(x)| \cdot \rho_y(y) + r_2 \cdot |f_2(x, y)|\]

where $r_i$ are certain reweighting factor to make sure all terms contribute same weights.

One can then estimate the integrals by averaging the observables $f_1(x)\rho_y(y)/p(x, y)$ and $f_2(x, y)/p(x, y)$.

The algorithm defaults to the standard Vegas algorithm if $r_0 = 1$ and $r_{i>0} = 0$.

Arguments

• integrand: A user-defined function evaluating the integrand. The function should be either integrand(var, config) or integrand(var, weights, config) depending on whether inplace is false or true respectively. Here, var are the random variables and weights is an output array to store the calculated weights. The last parameter passes the MC Configuration struct to the integrand, so that user has access to userdata, etc.

• measure: An optional user-defined function to accumulate the integrand weights into the observable. The function signature should be measure(var, obs, relative_weights, config). Here, obs is a vector of observable values for each component of the integrand and relative_weights are the weights calculated from the integrand multiplied by the probability of the corresponding variables.

The following are the snippets of the integrand and measure functions:

function integrand(var, config)
    # calculate your integrand values
    # return integrand1, integrand2, ...
end

function measure(var, obs, weights, config)
    # accumulates the weight into the observable
    # For example,
    # obs[1] = weights[1] # integral 1
    # obs[2] = weights[2] # integral 2
    # ...
end

Examples

The following command calls the VegasMC solver,

julia> integrate((x, c)->(x[1]^2+x[2]^2); var = Continuous(0.0, 1.0), dof = [[2,],], verbose=-1, solver=:vegasmc)
Integral 1 = 0.6640840471808533 ± 0.000916060916265263   (reduced chi2 = 0.945)