Integration with an adaptive harmonic mean algorithm

Numerically estimating the integral of functions in high dimensional spaces is a non-trivial task. A oft-encountered example is the calculation of the marginal likelihood in Bayesian inference, in a context where a sampling algorithm such as a Markov Chain Monte Carlo provides samples of the function. We present an Adaptive Harmonic Mean Integration (AHMI) algorithm. Given samples drawn according to a probability distribution proportional to the function, the algorithm will estimate the integral of the function and the uncertainty of the estimate by applying a harmonic mean estimator to adaptively chosen regions of the parameter space. We describe the algorithm and its mathematical properties, and report the results using it on multiple test cases.


Introduction
Sampling algorithms, such as Markov Chain Monte Carlo (MCMC), 1 are often used to generate samples distributed according to nontrivial densities in high dimensional spaces.Many algorithms have been developed that allow generating samples {Λ} from an unnormalized target density f (λ): In many applications, it is desirable or even necessary to be able to normalize the target density; i.e. to calculate where Ω is the support of f .This integral can be computationally very costly or impossible to compute using standard techniques; e.g. if the volume of the support where the target f is non-negligible occupies a very small region of the full support.
An important area where such integration is necessary is Bayesian inference. 2,3Bayes' formula, for a given model M , is P (λ|Data, M ) = P (Data|λ, M )P 0 (λ|M ) where here λ denotes the parameters of the model and the data are used to update probabilities for possible values of λ from prior probabilities P 0 (λ|M ) to posterior probabilities P (λ|Data, M ).The denominator is usually expanded using the Law of Total Probability and written in the form Z = P (Data|M ) = P (Data|λ, M )P 0 (λ|M )dλ .
Z goes by the name "evidence," or "marginal likelihood," and is an example of the type of integral that we want to be able to calculate (here the data are fixed and f (λ) = P (Data|λ, M )P 0 (λ|M )).An example use of Z is the calculation of a Bayes factor to compare two models: We are specifically interested in providing an algorithm applicable in a setting where samples are available from the target density f (λ) but with possibly no further recourse to generating more samples.For this purpose, we investigate the use of a modified harmonic mean estimator (HME). 4We introduce the use of a reduced integration region to improve the HME performance.After a description of the technique, we report on numerical investigations using Metropolis-Hastings MCMC samples.Our work has in common with 5 the use of the ratio of samples found in the limited support region to the number found in the full support, but we use a substantially different integral estimation technique.

Reduced Volume HME
We are interested in estimating the integral I from Eq. (1).We start by defining the integral 2050142-2 Int.J. Mod.Phys.A 2020.35.Downloaded from www.worldscientific.comby TECHNICAL UNIVERSITY OF MUNICH on 01/11/21.Re-use and distribution is strictly not permitted, except for Open Access articles.
Integration with an adaptive harmonic mean algorithm with ∆ ⊂ Ω a finite integration region, and the ratio Given our assumption that the sampling algorithm has successfully sampled from f (λ), we use the following as an estimator to our ratio which is the fraction of the total number of samples that fall within ∆ ⊂ Ω. Defining the normalized density over ∆ f∆ (λ) = f (λ) allows us to perform a harmonic mean calculation as follows: where V ∆ is the volume of the region defined by ∆.An estimator for this expectation value is the harmonic mean The HME for the reduced volume integral then follows as This calculation is performed directly from the values of the target density f (λ i ) given by the sampling algorithm, and does not require any extra sampling.An estimator for the integral over the full space Ω can then be written down as The task of estimating our integral, therefore reduces to choosing one or several subspaces ∆ -typically small regions around local modes of f (λ).The full space Ω over which the integration ought to be performed can be large or even infinite, while this does not affect the outcome of our integral estimate.We discuss the bias and uncertainty of this estimator in the following subsection.
In general, MCMC samples come with weights (e.g.repeated samples, with the weight being the number of repetitions).We therefore rewrite Eq. (11) as variation in the integral results is largest for small windows due to the small number of samples used, and for large windows due to the divergence of the harmonic mean estimator.We discuss the biases in the next section.

Bias and Uncertainty of the Estimator
We can estimate the bias and uncertainty on Î given in Eq. 11 by separately analyzing the behavior of r and X.As described below, we choose regions ∆ for which the range of target density values is moderate.For i.i.d.sampling, this would imply, via the Central Limit Theorem, that X follows a Normal distribution.Assuming we can approximate the distribution of X with a Normal distribution, we have where P ( X) is the probability distribution for X with mean µ and variance σ 2 estimated from the observed samples.Since X appears in the denominator in Eq. 11, this produces a bias in our integral of size σ2 X /μ 2 X .The fractional uncertainty in our integral estimator is σX /μ X .
The estimator r will also typically follow approximately a Normal distribution with parameters that can be estimated from i.i.d.sampling and Binomial statistics with w i the weights assigned to the samples at parameter values λ i and W Ω = i w i the sum of all weights.The use of weights also allows this technique to be applied to samples obtained from, for example, importance sampling.

Illustration of the technique
To illustrate our technique for applying harmonic mean integration, we consider the unit normal distribution p(x) = N (x|µ = 0, σ = 1).A fixed number of samples (3 • 10 3 ) was generated from directly sampling the unit normal distribution, and Eq. ( 11) was used to calculate the integral for different subregions ∆.These regions are defined as windows of x centered on 0 and varied in width from 0.02 up to 9.This was repeated 1.5 • 10 4 times, and the mean and standard deviation were evaluated.Figure 1 shows the results of the integration as a function of window size.As is seen in the figure, harmonic mean integration applied to a finite region gives an accurate value for the integral over a wide range of sampling windows.The variation in the integral results is largest for small windows due to the small number of samples used, and for large windows due to the divergence of the harmonic mean estimator.We discuss the biases in the next section.

Bias and uncertainty of the estimator
We can estimate the bias and uncertainty on Î given in Eq. ( 11) by separately analyzing the behavior of r and X.As described below, we choose regions ∆ for which the range of target density values is moderate.For i.i.d.sampling, this would imply, via the Central Limit Theorem, that X follows a normal distribution.Assuming we can approximate the distribution of X with a normal distribution, we have , where P ( X) is the probability distribution for X with mean µ and variance σ 2 estimated from the observed samples.Since X appears in the denominator in Eq. (11), this produces a bias in our integral of size σ2 X /μ 2 X .The fractional uncertainty in our integral estimator is σX /μ X .
The estimator r will also typically follow approximately a normal distribution with parameters that can be estimated from i.i.d.sampling and binomial statistics as Since r also appears in the denominator, it will also produce a bias in our integral of size σ2 r /μ 2 r .The fractional uncertainty in our integral estimator from r is, in the approximation of i.i.d.sampling and a normal distribution, σr /μ r .
We can therefore write down an explicit correction factor b that we apply to the integral estimate Î by multiplication.The correction is illustrated in Fig. 2   as from r dominates, as the binomial uncertainty is largest for small numbers of samples.With the bias correction applied, already for as few as ≈ 20 samples inside the integration volume (corresponding to roughly |x| < 0.01), accurate results are produced, while without the correction factor applied much large windows, starting at around |x| < 0.5, are necessary.
Toward larger windows the bias produced from X become dominant, as the range of values of f of the contained samples grows.The bias correction successfully mitigates this effect as illustrated.
Once the window size exceeds the space where samples are present, the integral starts diverging.For our example using 3 • 10 3 samples, we expect those to cover a region up to only |x| ≈ 3.58, which explains well the observed trend.
Many samplers such as MCMC algorithms generate strong correlations among samples and using the binomial uncertainty discussed here can be inaccurate.We therefore also numerically evaluate the uncertainty as described in detail below.The integration regions are chosen such that the bias correction can be neglected.

Relation to other techniques for evidence calculation
A variety of techniques to calculate the marginal likelihood in Bayesian calculations have been successfully developed.A summary can be found in Ref. 8, where a number of MCMC related techniques are reviewed, including Laplace's method, 9 harmonic mean estimation, 4 Chib's method, 10 annealed importance sampling techniques, 11,12 nested sampling 13 and thermodynamic integration methods. 14,15Only the HME and Laplace techniques allow the direct estimation of the evidence from available samples, and the Laplace technique makes the unwanted assumption that the target density is a single multivariate Gaussian.The Laplace approximation fails badly for multimodal distributions or distributions with significant probability mass in the tails of the distribution.
In the Bayesian literature, 4 the HME for the evidence, Z, is formulated as follows: 1 The difference with our formulation is that the prior P 0 (λ|M ) is not included in the expectation value, and the integration volume does not appear.In order to perform this HME calculation, it is necessary to know separately the value of the likelihood and the prior at the sampling points.The HME method has been strongly criticized (even called "worst Monte Carlo method ever" 16 ), since the evaluation of the denominator in Eq. ( 14) can have very large variance.Reducing the variance by limiting the integration region is not as straightforward as in our formulation since it requires an integral of the prior function over the reduced support whereas this is not needed in our formulation.
Integration with an adaptive harmonic mean algorithm

An Adaptive Harmonic Mean Integration Algorithm
Adaptive Harmonic Mean Integration (AHMI) uses the HME on multiple subregions ∆ i to estimate the integral of f (λ) over its full support Ω.In this section we present our example algorithm in detail and will show benchmark tests on several distributions in Sec. 4. As discussed previously, defining a set of suitable regions is crucial in obtaining a robust and unbiased estimate of the integral of f (λ).In particular, to avoid biasing the result, it is essential not to use the same elements of the sample set {Λ} for both the definition of ∆ i and estimates of the integral Îi .
(2) Split into two mutually exclusive, uncorrelated and equally sized subsets A and B.
(3) Generate seed points λ seed i (via space partitioning tree).The general flow of the AHMI algorithm, including the procedure of defining the ∆ i , is summarized in Fig. 3.The various involved steps are discussed in more technical detail in the following subsections.
(   The general flow of the AHMI algorithm, including the procedure of defining the ∆ i , is summarized in Fig. 3.The various involved steps are discussed in more technical detail in the following subsections.

Samples, preprocessing and splitting
We start with a given set of samples {Λ} that we assume are drawn according to the probability distribution proportional to our function f , obtained for example from MCMC sampling.
In order to de-correlate the sample space we apply a whitening transformation.In general, a whitening transformation maps a set of random variables with a known nonsingular covariance matrix to a new set of variables with a covariance matrix equal to I. A Cholesky Decomposition is used to whiten the samples, and the AHMI estimator for the integral becomes where det R is the determinant of the whitening matrix and the primed symbols represent the quantities in the transformed space.In the following we drop this explicit addition of prime symbols and work in the whitened space (unless otherwise stated).The full set of samples is then divided into two equally sized and mutually exclusive subsets A and B.

Hyper-rectangle generation
We illustrate the hyper-rectangle generation steps in more detail using a twodimensional Gaussian shell example with distribution In our examples, we use the following settings: radius r = 5, width ω = 2 and c = 0.The integration region extends from [−25, 25] in each dimension.Samples from this distribution are shown in Fig. 4(a).
The algorithm starts by creating seed points around which to construct the integration regions ∆ i .These points should lie in areas of high density and should result in broadly distributed starting points.In order to limit computation time, a simple space-partitioning tree is used to divide the whitened space into subsets of nonoverlapping regions.A tree is constructed by performing cuts in every dimension in such a way to have an equal number of samples on the left and right leaves.The number of cuts in each parameter axes is determined by the total number of samples and the number of dimensions and is chosen in a way, that the number of samples in each leaf does not exceed 200.An example of such a space-partitioning tree is   the following condition: Although this threshold is user-selectable, by default it is equal to 500 (an example of integration with different threshold values is shown in Sec.4.1).
To create M regions for integration, we select the seed point λ seed  the following condition: Although this threshold is user-selectable, by default it is equal to 500 (an example of integration with different threshold values is shown in Sec.4.1).
To create M regions for integration, we select the seed point λ seed i with the overall largest value f (λ seed ) and follow the steps below.Then we recursively repeat the  the following condition: Although this threshold is user-selectable, by default it is equal to 500 (an example of integration with different threshold values is shown in Sec.4.1).
To create M regions for integration, we select the seed point λ seed i with the overall largest value f (λ seed i ) and follow the steps below.Then we recursively repeat the  the following condition: Although this threshold is user-selectable, by default it is equal to 500 (an example of integration with different threshold values is shown in Sec.4.1).
To create M regions for integration, we select the seed point λ seed i with the overall largest value f (λ seed i ) and follow the steps below.Then we recursively repeat the  shown in Fig. 4(b).For each partition i, the sample with the largest function value contained in that partition is defined as seed point λ seed i .
The following steps produce hyper-rectangle-shaped regions suitable for AHMI.In order to limit this variance and ensure numerical stability, the ratio between the highest and the lowest probability of samples inside a hyper-rectangle is bound by the following condition: Although this threshold is user-selectable, by default it is equal to 500 (an example of integration with different threshold values is shown in Subsec.4.1).
To create M regions for integration, we select the seed point λ seed i with the overall largest value f (λ seed i ) and follow the steps below.Then we recursively repeat the same procedure M − 1 times using the remaining seed points.
(1) The algorithm starts by building a small hyper-cube ∆i around the selected seed point λ seed i , see Fig. 4(c).
(2) This hyper-cube is then incrementally either increased or decreased in size, until the probability ratio of contained samples matches the threshold t within some tolerance, or until it contains more than 1% of the total samples.(3) The faces of the hyper-cube are then iteratively adjusted (expand or contract), to adapt to the density of the contained samples, while enforcing the condition f max /f min ≤ t.This step turns the D-dimensional hyper-cube into a Ddimensional hyper-rectangle.This hyper-rectangle adaptation algorithm continues as long as changes to the hyper-rectangle's faces are accepted.The stopping criterion is based on the fraction of samples accepted or rejected comto expectation from the volume change.However, the hyper-rectangle adaption algorithm always ensures that no modification to the hyper-rectangle's faces are made if such a modification would result in f max /f min > t.An in-depth description of the algorithm can be found in Ref. 19 and an opensource implementation of it is available as part of the BAT.jl software package. 20

Integral estimates
Once M integration regions are defined, we can compute the integral estimates ÎA i for each ∆ B i , according to Eq. ( 12).The procedure is the same for ÎB i , so we shall drop the superscripts A and B in the following two subsections.The two resulting, separate estimates ÎA and ÎB will then be combined in Subsec.3.5 to obtain the final estimate.
From the distribution of all estimates Îi we select only the 68% central percentile to reject outliers -a procedure that was empirically found to work well.This is indicated in Fig. 4(d) labeled as "accepted" and "rejected" rectangles.We proceed to combine the remaining estimates Îi into a single estimate Î using a robust and unbiased estimator for the combination of correlated measurements as suggested in Ref. 17 where the weights w i are defined as: The variances σ2 i ≡ σii and covariances σij of the mean assigned to integration regions ∆ i and ∆ j are estimated in the next section.

Covariance estimate
The following procedure is used to estimate the covariance between individual integral estimates Îi : (1) We partition {Λ} into a number S of subsets {Λ 1 }, {Λ 2 } • • • {Λ S }, chosen in a way that reduces their correlation.The default value for S is 10.(2) Separate estimates Îi,k (k enumerates the S partitions) of the integral are then performed for all sample subsets {Λ k } resulting in S integral estimates for each subspace ∆ i : (3) The covariance of the separate integral values then is Îi,k .
(4) The estimate of the covariance of Īi ≈ Îi and Īj ≈ Îj is then σ2

Alternative covariance estimation method
As S cannot be large without requiring a very large number of samples, the entries of σij tend to be noisy.In practice, we have found this to be tolerable.However, we have developed an optional alternative procedure to estimate σij , for use cases where the total number of samples does not allow for a sufficient number of subsets.

The diagonal elements σ2
ii can be estimated based on the theoretical uncertainties derived in Subsec.2.2, using the effective number of samples (estimated from the autocorrelation of the samples).
The off-diagonal elements σi =j are estimated heuristically, based on the overlap of the hyper-rectangles (quantified as the total weight of the samples they share): where with W ∩ (i, j) and W ∪ (i, j) being the number of weighted samples in ∆ i ∩ ∆ j and ∆ i ∪ ∆ j , respectively.Some details of this alternative procedure are still under development, so we do not present an in-depth description and formal arguments for it is validity here.We have observed close agreement with the empirical covariance (step 3) in many cases and hope to further refine this method in the future.

Final AHMI integral estimate
As we have now obtained two values of Î and two variances σ 2 ( Î) from the two sets {Λ A } and {Λ B }, we combine those into the final result like with variance estimate

Benchmark Examples
To validate our algorithm, we apply it to estimate the integral of several test functions in varying dimensionality (up to d = 25) for which an analytic (or accurate numerical) solution of the integral value is available.The test functions were chosen to pose different challenges to the algorithm.For our first test problem, we start with the canonical example of a multivariate normal distribution.For the second test case we look at Gaussian shells, for which the mode of the distribution does not lie on a single point but has infinite modes on (d − 1)-dimensional surface.Next we study the heavy-tailed Cauchy distribution with multiple modes, and in the end explore the asymmetric "Funnel" function.Additional information of the number of integration volumes used and the computation time are only provided for the first example, as these are very similar for the other three examples.
The samples {Λ} on which our integration is based are obtained from Metropolis-Hastings MCMC, 20 and in the case of the multivariate normal from i.i.d.sampling.The sample size is fixed to 2 • 10 6 for Gaussian shells example and it is equal to 10 6 for other examples.

Multivariate normal distribution
The first test case is a unit normal distribution (centered at zero, width one) in two up to 25 dimensions.The samples input to AHMI are obtained from i.i.d.sampling.The resulting integral estimates are shown in Fig. 5  standard error for t = 500 (220 trials) is 0.52 ± 0.03.For the threshold value t = 100 we get an unbiased integral estimate up to 14 dimensions, and after that hyperrectangles can no longer be created.The default value of t = 500 was used for the remaining examples below.
In Fig. 6 we show the execution time of the algorithm, and the number of   For both threshold values, t = [500, 1000], we get unbiased and consistent results up to around 21 dimensions, after which results start to become positively biased for t = 1000.This bias seems to be intrinsic to the method itself or our implementation of the algorithm and will be the subject of future studies.The coverage for one standard error for t = 500 (220 trials) is 0.52 ± 0.03.For the threshold value t = 100 we get an unbiased integral estimate up to 14 dimensions, and after that 13 t = 1000.This bias seems to be intrinsic to the method itself or our implementation of the algorithm and will be the subject of future studies.The coverage for one standard error for t = 500 (220 trials) is 0.52 ± 0.03.For the threshold value t = 100 we get an unbiased integral estimate up to 14 dimensions, and after that hyperrectangles can no longer be created.The default value of t = 500 was used for the remaining examples below.
In Fig. 6 we show the execution time of the algorithm, and the number of   hyper-rectangles can no longer be created.The default value of t = 500 was used for the remaining examples below.
In Fig. 6 we show the execution time of the algorithm, and the number of hyper-rectangles used for integration for the threshold value t = 500.The execution time rises with the number of dimensions almost linearly.The change in slope at low dimensionality is likely due to CPU caching behavior.The number of hyperrectangles starts to decay after 18 dimensions indicating that there exist fewer hyper-rectangles that satisfy Eq. ( 18).

Gaussian shell distribution
The functional form was given in Eq. ( 17) and an example distribution in the first two dimensions is shown in Fig. 7.The AHMI algorithm results (Fig. 8) show a ur as for the multivariate normal distribution for this more complicate n.However, integration was possible up to 17 dimensions.Up to that gral estimates, including errors, are well behaved, with a coverage f d error of 0.48 ± 0.04 (160 trials).
ultimodal Cauchy Distribution uchy distribution, with its heavy tails, is a notoriously difficult proble re to point out possible weaknesses of our algorithm.We further in xity for the hyper-volume creation process by using four separate, s distributions creating multiple modes.The functional form can be w    similar behavior as for the multivariate normal for this more complicated test function.However, integration was possible up to 17 dimensions.Up to that point the integral estimates, including errors, are well behaved, with a coverage for one standard error of 0.48 ± 0.04 (160 trials).

Multimodal Cauchy distribution
The Cauchy distribution, with its heavy tails, is a notoriously difficult problem and used here to point out possible weaknesses of our algorithm.We further increase complexity for the hyper-volume creation process by using four separate, shifted Cauchy distributions creating multiple modes.The functional form can be written as where µ = 1, σ = 0.2 and n is a dimensionality of λ.An example of this target distributions is provided in Fig. 9.The integration region extends from [−8, 8] in each dimension.Our results are collected in Fig. 10 and indicate that integration is possible up to seven dimensions, given the fixed sample size of 10

# Dimensions
Fig. 8. AHMI value for the Gaussian shell distribution as ratio to the true integral shown as a function of dimensionality.The solid line gives the mean result over ten independent trials and the shaded band show the standard deviation of these trials.The dashed lines show the average errors reported by AHMI.The samples are obtained from Metropolis-Hastings MCMC and the sample size is 2 • 10 6 .as

Conclusion
We have developed an Adaptive Harmonic Mean Integration (AHMI) algorithm that can be used to integrate a non-normalized density function using the samples

Funnel distribution
The final problem we study is the so-called "Funnel" distribution, that is described in Ref. 21.The functional form of this distribution can be written as where a = 1, b = 0.5 and n is a dimensionality of λ.An example of the distribution in its first three dimensions on the parameter range [−50, 50] is provided in Fig. 11.
The results (Fig. 12) show a similar performance as for the previous example with reliable estimates up to seven dimensions, with a coverage of 0.41 ± 0.04 (120 trials, one standard error).

Conclusion
We have developed an Adaptive Harmonic Mean Integration (AHMI) algorithm that can be used to integrate a non-normalized density function using the samples drawn according to the probability distribution proportional to the function.The fundamental assumption is that the sampling algorithm has faithfully produced samples from this distribution.Given this, the AHMI algorithm can be used to produce both an estimate of the integral of the function over its full support as well as an estimate of the uncertainty of the integral.In this first implementation of the AHMI algorithm, finite hyper-rectangles are generated in the whitened space of the samples covering the full support.The adaptive algorithm ensures that the range of function values enclosed by the hyper-rectangles is limited such that the variance of the integral results are moderate.This allows for reliable results both for the integral values as well as for reliable uncertainty estimates.
The algorithm has been tested on a number of examples and found to produce reliable and unbiased integral estimates up to around 20 dimensions for the first two test problems, and up to seven for the latter two, more difficult test cases.The reported errors provide a useful measure of uncertainty, while slightly under covering at around 40-50% (expecting 68%).The use of hyper-rectangles however limits the applicability to a not-too-large number of dimensions (≈ 20 in the case of the multivariate normal distribution) because a large fraction of the volume is in the corners of the hyper-rectangles.In cases where we have a unimodal distribution, we anticipate that regions defined as shells between hyper-spheres of different radii would be more appropriate.We leave this for a future development of the algorithm.

) 2050142- 3 Fig. 1 .
Fig. 1.Demonstration of the reduced harmonic mean technique with the unit normal distribution.The left panel shows in blue the mean (solid line) and the region covering 68 % of results (from 1.5 • 10 4 repeated trials) of the integral estimates as a function of the window extent.The true integral value is 1.0, indicated in red.The gray shaded distribution shows the unit normal pdf (right y-axis scale).For the three window extents indicated by the black, vertical lines (Window 1: |x| < 0.24, Window 2: |x| < 1.61, Window 3: |x| < 4.20) the full distribution of integral estimates from the 1.5 • 10 4 trials are plotted on the right side as histograms.

Fig. 1 .
Fig. 1. (Color online) Demonstration of the reduced harmonic mean technique with the unit normal distribution.The left panel shows in blue the mean (solid line) and the region covering 68% of results (from 1.5 • 10 4 repeated trials) of the integral estimates as a function of the window extent.The true integral value is 1.0, indicated in red.The gray shaded distribution shows the unit normal pdf (right y-axis scale).For the three window extents indicated by the black, vertical lines (Window 1: |x| < 0.24, Window 2: |x| < 1.61, Window 3: |x| < 4.20) the full distribution of integral estimates from the 1.5 • 10 4 trials are plotted on the right side as histograms.

Fig. 2 .
Fig.2.Left: The average uncorrected (Eq.11) and corrected (b • Î) integrals as a function of the extent of the window used to accept samples.The results are from averaging 1.5 • 10 4 integration results, where each integration test was performed from 3 • 10 3 samples in the full function range.Right: Individual contributions to the bias correction from the Binomial (r) and the 1/f terms.

Fig. 2 .
Fig. 2. Left: The average uncorrected (Eq.(11)) and corrected (b• Î) integrals as a function of the extent of the window used to accept samples.The results are from averaging 1.5 • 10 4 integration results, where each integration test was performed from 3 • 10 3 samples in the full function range.Right: Individual contributions to the bias correction from the binomial (r) and the 1/f terms.2050142-5 Int.J. Mod.Phys.A 2020.35.Downloaded from www.worldscientific.comby TECHNICAL UNIVERSITY OF MUNICH on 01/11/21.Re-use and distribution is strictly not permitted, except for Open Access articles.

( 4 ).( 5 )( 6 )
Create a small hyper-cube ∆i around each seed point λ seed i Adjust the faces of the hyper-cubes, resulting in hyper-rectangles.Use the resulting hyper-rectangles ∆ B i created with sample {Λ B } for the harmonic mean estimates ÎA i on {Λ A } and vice versa.

( 7 )( 8 )
Use individual estimates Îi to compute combined estimate and variance for A and B. Compute final estimate by combining the independent estimates ÎA and ÎB weighted by their variance.August 3, 2020 10:32 WSPC/INSTRUCTION FILE main 7

Fig. 3 .
Fig. 3. Overview of the different steps in the AHMI algorithm, including the procedure of finding subvolumes ∆ i and computing integral estimates Îi .2050142-7 Int.J. Mod.Phys.A 2020.35.Downloaded from www.worldscientific.comby TECHNICAL UNIVERSITY OF MUNICH on 01/11/21.Re-use and distribution is strictly not permitted, except for Open Access articles.
(a) MCMC samples after whitening transformation.The size of each point is proportional to its weight.(b) Two dimensional space-partitioning tree.All regions contain an equal number of unique samples.(c) Initial hypercubes ∆ around seed points.(d) Hyper-rectangles after adjustments.Only the red hyper-rectangles are used in the final integral estimate.

Fig. 4 .
Fig. 4. Process of finding integration regions in a two dimensional example for the Gaussian shells test function.

i 9 (
with the overall seed (a) MCMC samples after whitening transformation.The size of each point is proportional to its weight.a) MCMC samples after whitening transformation.The size of each point is proportional to its weight.(b) Two dimensional space-partitioning tree.All regions contain an equal number of unique samples.(c) Initial hypercubes ∆ around seed points.(d) Hyper-rectangles after adjustments.Only the red hyper-rectangles are used in the final integral estimate.

Fig. 4 .
Fig. 4. Process of finding integration regions in a two dimensional example for the Gaussian shells test function.
(b) Two-dimensional space-partitioning tree.All regions contain an equal number of unique samples.9 (a) MCMC samples after whitening transformation.The size of each point is proportional to its weight.(b) Two dimensional space-partitioning tree.All regions contain an equal number of unique samples.(c) Initial hypercubes ∆ around seed points.(d) Hyper-rectangles after adjustments.Only the red hyper-rectangles are used in the final integral estimate.

Fig. 4 .
Fig. 4. Process of finding integration regions in a two dimensional example for the Gaussian shells test function.
(c) Initial hypercubes ∆ around seed points.9 (a) MCMC samples after whitening transformation.The size of each point is proportional to its weight.(b) Two dimensional space-partitioning tree.All regions contain an equal number of unique samples.(c) Initial hypercubes ∆ around seed points.(d) Hyper-rectangles after adjustments.Only the red hyper-rectangles are used in the final integral estimate.

Fig. 4 .
Fig. 4. Process of finding integration regions in a two dimensional example for the Gaussian shells test function.
(d) Hyper-rectangles after adjustments.Only the red hyper-rectangles are used in the final integral estimate.

Fig. 4 .
Fig. 4. Process of finding integration regions in a two-dimensional example for the Gaussian shells test function.

Figure 4 (
Figure 4(d) shows the resulting set of M hyper-rectangles for our example.An in-depth description of the algorithm can be found in Ref.19 and an opensource implementation of it is available as part of the BAT.jl software package.20
as a function of the dimensionality for three different threshold values t = [100, 500, 1000].2050142-12 Int.J. Mod.Phys.A 2020.35.Downloaded from www.worldscientific.comby TECHNICAL UNIVERSITY OF MUNICH on 01/11/21.Re-use and distribution is strictly not permitted, except for Open Access articles.

Fig. 5 .
Fig. 5. AHMI value for the multivariate normal distribution as ratio to the true integral shown as a function of dimensionality for three threshold values t = [100, 500, 1000].The solid black lines give the mean result over ten independent trials and the shaded bands show the standard deviation of these trials.The dashed lines show the average errors reported by AHMI.Samples are obtained from i.i.d.sampling, 10 6 for each run.

Fig. 6 .
Fig.6.Left: Average AHMI execution time for the multivariate normal distribution (t = 500) in total CPU seconds (run on a system with a 2.3 GHz Intel Xeon 6140 processor, SPECrate2017-FP equivalent ca.180) as a function of dimensionality.Right: Average number of hyper-rectangles used by AHMI for computing the integral estimate for the multivariate normal distribution as a function of dimensionality.

Fig. 5 .
Fig. 5. AHMI value for the multivariate normal distribution as ratio to the true integral shown as a function of dimensionality for three threshold values t = [100, 500, 1000].The solid black lines give the mean result over ten independent trials and the shaded bands show the standard deviation of these trials.The dashed lines show the average errors reported by AHMI.Samples are obtained from i.i.d.sampling, 10 6 for each run.

Fig. 5 .
Fig. 5. AHMI value for the multivariate normal distribution as ratio to the true integral shown as a function of dimensionality for three threshold values t = [100, 500, 1000].The solid black lines give the mean result over ten independent trials and the shaded bands show the standard deviation of these trials.The dashed lines show the average errors reported by AHMI.Samples are obtained from i.i.d.sampling, 10 6 for each run.
Fig.6.Left: Average AHMI execution time for the multivariate normal distribution (t = 500) in total CPU seconds (run on a system with a 2.3 GHz Intel Xeon 6140 processor, SPECrate2017-FP equivalent ca.180) as a function of dimensionality.Right: Average number of hyper-rectangles used by AHMI for computing the integral estimate for the multivariate normal distribution as a function of dimensionality.

2
One and two dimensional distributions of samples along the first two dimension shell target function.

Fig. 7 .
Fig. 7. One-and two-dimensional distributions samples along the first two dimensions of the Gaussian shell target function.2050142-14 Int.J. Mod.Phys.A 2020.35.Downloaded from www.worldscientific.comby TECHNICAL UNIVERSITY OF MUNICH on 01/11/21.Re-use and distribution is strictly not permitted, except for Open Access articles.

Fig. 8 . 4 Fig. 9 .
Fig. 8. AHMI value for the Gaussian shell distribution as ratio to the true integral shown as function of dimensionality.The solid line gives the mean result over ten independent trials and the shaded band show the standard deviation of these trials.The dashed lines show the average errors reported by AHMI.The samples are obtained from Metropolis-Hastings MCMC and the sample size is 2 • 10 6 .

Fig. 8 .
Fig. 8. AHMI value for the Gaussian shell distribution as ratio to the true integral shown as a function of dimensionality.The solid line gives the mean result over ten independent trials and the shaded band show the standard deviation of these trials.The dashed lines show the average reported by AHMI.The samples are obtained from Metropolis-Hastings MCMC and the sample size is 2 • 6 .

Fig. 9 .
Fig. 9.One and two dimensional distributions of samples along the first four dimensions of the multimodal Cauchy target function.

Fig. 9 .Fig. 10 .
Fig. 9. One-and two-dimensional distributions of samples along the first four dimensions of the multimodal Cauchy target function.

Fig. 10 .
Fig. 10.AHMI value for the Cauchy distribution as ratio to the true integral shown as a function of dimensionality.The solid line gives the mean result over twenty independent trials and the shaded band show the standard deviation of these trials.The dashed lines show the average errors reported by AHMI.2050142-16 Int.J. Mod.Phys.A 2020.35.Downloaded from www.worldscientific.comby TECHNICAL UNIVERSITY OF MUNICH on 01/11/21.Re-use and distribution is strictly not permitted, except for Open Access articles.

Fig. 11 .Fig. 12 .
Fig. 11.One and two dimensional distributions of samples along the first three dimensions of the Funnel target function.

Fig. 11 .Fig. 11 .
Fig. 11.One-and two-dimensional distributions of samples along the first three dimensions of the Funnel target function.

Fig. 12 .
Fig.12.AHMI value for the "Funnel" distribution as ratio to the true integral shown as a function of dimensionality.The solid line gives the mean result over twenty independent trials, with the shaded band showing the standard deviation of these trials.The dashed lines show the average errors reported by AHMI.

Fig. 12 .
Fig. 12. AHMI value for the "Funnel" distribution as ratio to the true integral shown as a function of dimensionality.The solid line gives the mean result over twenty independent trials, with the shaded band showing the standard deviation of these trials.The dashed lines show the average errors reported by AHMI.2050142-17 Int.J. Mod.Phys.A 2020.35.Downloaded from www.worldscientific.comby TECHNICAL UNIVERSITY OF MUNICH on 01/11/21.Re-use and distribution is strictly not permitted, except for Open Access articles.
6. Left: Average AHMI execution time for the multivariate normal distribution (t = 500) in total CPU seconds (run on a system with a 2.3 GHz Intel Xeon 6140 processor, SPECrate2017-FP equivalent ca.180) as a function of dimensionality.Right: Average number of hyper-rectangles used by AHMI for computing the integral estimate for the multivariate normal distribution as a function of dimensionality.Int.J. Mod.Phys.A 2020.35.Downloaded from www.worldscientific.comby TECHNICAL UNIVERSITY OF MUNICH on 01/11/21.Re-use and distribution is strictly not permitted, except for Open Access articles.
Int. J. Mod.Phys.A 2020.35.Downloaded from www.worldscientific.comby TECHNICAL UNIVERSITY OF MUNICH on 01/11/21.Re-use and distribution is strictly not permitted, except for Open Access articles.
6.For this range, the AHMI results are very reliable.The coverage (120 trials, one standard error) for this function is 0.42 ± 0.05.2050142-15