# A novel sequential method for building upper and lower bounds of moments of distributions

Solal Martin, Emilie Chouzenoux, Víctor Elvira

<figure><img src="image_1.png" /></figure>

## To cite this version:

Solal Martin, Emilie Chouzenoux, Víctor Elvira. A novel sequential method for building upper and lower bounds of moments of distributions. 2025. hal-05392970

HAL Id: hal-05392970

https://inria.hal.science/hal-05392970v1

Preprint submitted on 2 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A novel sequential method for building upper and lower bounds of moments of distributions

Solal Martin<sup>a,*</sup>, Emilie Chouzenoux<sup>b</sup>, Víctor Elvira<sup>c</sup>

## Abstract

Approximating integrals is a fundamental task in probability theory and statistical inference, and their applied fields of signal processing, and Bayesian learning, as soon as expectations over probability distributions must be computed efficiently and accurately. When these integrals lack closed-form expressions, numerical methods must be used, from the Newton-Cotes formulas and Gaussian quadrature, to Monte Carlo and variational approximation techniques. Despite these numerous tools, few are guaranteed to preserve majoration/minoration inequalities, while this feature is fundamental in certain applications in statistics.

In this paper, we focus on the integration problem arising in the estimation of moments of scalar, unnormalized, distributions. We introduce a sequential method for constructing upper and lower bounds on the sought integral. Our approach leverages the majorization-minimization framework to iteratively refine these bounds, in an enveloped principle. The method has proven convergence, and controlled accuracy, under mild conditions. We demonstrate its effectiveness through a detailed numerical example of the estimation of a Monte-Carlo sampler variance in a Bayesian inference problem.

**Keywords:** Integral Approximation, Statistical Optimization, Majoration-Minimization methods, Importance Sampling.

<footnote><sup>a</sup>Toulouse School of Economics, France</footnote>

<footnote><sup>b</sup> OPIS, Inria Saclay, Université Paris-Saclay, France</footnote>

<footnote><sup>c</sup>School of Mathematics, University of Edinburgh, UK</footnote>

<footnote>*Corresponding author
Email address: solal.martin@tse-fr.eu (Solal Martin)</footnote>
# 1. Introduction and Problem Statement

## 1.1. Introduction

The approximation of intractable integrals has long been for long a fundamental problem in mathematical sciences, leading to extensive research and the development of various methods [1, 2, 3]. Classical approaches include Newton-Cotes formulas, such as the trapezoidal and Simpson's rules, which approximate the integrand using polynomial interpolation over equally spaced points [4]. While these methods are easy to implement, they can have low accuracy for highly oscillatory functions due to the Runge's phenomenon [5]. A more advanced approach relies on Gaussian quadrature, which optimizes deterministics integration points and weights to achieve higher accuracy with fewer function evaluations [6, 7]. MATLAB's integrator [8], for example, extends this idea using the Gauss-Kronrod quadrature rule, which enhances accuracy by adding additional nodes to an existing Gaussian quadrature formula. For integrals over infinite domains or with singularities, adaptive quadrature methods dynamically refine the integration mesh based on local error estimates, improving efficiency and precision [9].

When dealing with multi-dimensional settings, traditional quadrature methods rapidly become impractical, and Monte Carlo methods provide an alternative approach by estimating integrals through random sampling [1, 10]. Importance sampling methods, an important family of Monte Carlo algorithms, improve efficiency by favoring regions that contribute most to the integral's value [1, 11]. Another powerful approach in high-dimensional integration problems is sparse grid integration, which reduces the curse of dimensionality by combining lower-order quadrature rules adaptively [12]. Each of these methods has its advantages and limitations, making their selection highly dependent on the problem at hand. A more detailed review of these numerical integration techniques can be found in standard numerical analysis texts, such as the reference book [13], for instance.

An important class of integration problem is the computation a unidimensional moment integral of a given target distribution. Beyond that, the assessment of moments methods often require the evaluation of the variance of the computed estimator, to evaluate the method's efficiency. One widely used diagnostic in this context is the Effective Sample Size (ESS), which provides a quantitative measure of how well the weighted samples represent the target distribution. ESS is relevant in importance sampling and sequential Monte Carlo methods, where weight degeneracy can significantly affect estimator quality [14, 15]. Other techniques include variance reduction strategies
such as control variates, stratified sampling, and resampling methods, which aim to improve estimator stability and convergence [16, 17]. The variance of the estimator is often expressed as an integral as hard as the original problem. Even more the lower and upper bounds of the variance, also involve intractable integrals. Numerical approximation can be used. However, most state-of-the-art integration methods come with significant approximation errors, and, do not guarantee to preserve inequalities. This paper tackles this challenge by providing an original sequential approach to compute tight bounds for moment integrals.

## 1.2. Problem Statement

Let $f: \mathbb{R} \to \mathbb{R}$ and $\phi: \mathbb{R} \to \mathbb{R}$, differentiable. We consider in this paper the following moment integral:

$$ \mathcal{I} = \int_{-\infty}^{+\infty} f(x)\pi(x)dx, \quad \text{with} \quad \pi(x) = \exp(-\phi(x)). \qquad (1) $$

No normalization assumption is made on $\pi$, that is $Z \equiv \int_{-\infty}^{+\infty} \pi(x)dx$ might be different to one. For most choices of $f$ and $\phi$, $\mathcal{I}$ does not have a closed-form expression and should be approximated by a numerical method. The goal of this paper is to provide an efficient iterative strategy for computing tight lower and upper bounds of $\mathcal{I}$, that we define as

$$ \underline{\mathcal{I}} \le \mathcal{I} \le \bar{\mathcal{I}}. \qquad (2) $$

For the ease of presentation, we will consider throughout the paper, that function $f$ is a monomial function, i.e.,

$$ (\forall x \in \mathbb{R}) \quad f(x) = x^k, \text{ for some } k \in \mathbb{N}. \qquad (3) $$

If $k = 0$, the problem amounts to computing the normalization constant of $\pi$, since $\mathcal{I} = \mathcal{Z}$. Extension of the method to a polynomial test function $f$, will be discussed in the Remark 1, in the end of Section 4. Throughout the paper, we assume that the problem is well-posed, that is $f$ and $\pi$ are such that $\mathcal{I} \in \mathbb{R}$.

## 1.3. Paper Contributions and Outline

This paper presents a new sequential method to compute the bounds $(\underline{\mathcal{I}}, \bar{\mathcal{I}})$ in (2) with controlled precision. Our approach first stands on the construction of piecewise Gaussian upper and lower approximations for $\pi$, by relying on the powerful principle of majoration-minimization [18, 19, 20]. We then design a recursive method, to compute upper and lower envelopes, from the set of piecewise
approximations. Our approach hence yields the construction of two sequences of upper and lower bounds for the integral of interest (1), with increasing, and controlled, accuracy, and both proved to converge to the exact integral value $\mathcal{I}$. Although integral $\mathcal{I}$ could be computed, in special cases of $\pi$, using known results such as the moments of truncated Gaussian distributions, our method applies more generally to any unnormalized distribution satisfying the assumptions detailed in Section 3. We then specialize our method to the computation of bounds on the variance of an importance sampling estimator. Finally, we illustrate its practical effectiveness of our method, and compare it to the state-of-the-art envelope-based integration technique from [13, Sec.5.4], in an application of statistical inference, namely a Bayesian classification problem.

The outline of the paper is as follows. In Section 2, we introduce our construction of upper and lower Gaussian approximations for $\pi$, under reasonable assumptions on $\phi$. We deduce a preliminary formula of rough lower and upper bounds $\underline{\mathcal{I}}$ and $\overline{\mathcal{I}}$ for $\mathcal{I}$. Then, in Section 3, we propose an iterative algorithm which aims at producing $(\underline{\mathcal{I}}^{(n)})_{n \in \mathbb{N}}$ and $(\overline{\mathcal{I}}^{(n)})_{n \in \mathbb{N}}$, that are sequences of refined lower and upper bounds for the integral of interest, $\mathcal{I}$. We prove in Section 4 that both sequences converge to $\mathcal{I}$ as $n$ approaches infinity, and we provide a quantitative control of the error along iterations. In Section 5, we formulate the problem of variance estimation of a Monte-Carlo importance sampling approach and the specialization of our method to this case. Finally, in Section 6, we illustrate the performance of our method, and compare it to the state-of-the-art, in the use-case context motivated in Sec. 5.

## 2. A first construction of integral bounds

### 2.1. Introduction

Let us first define the notion of tangent majorant (resp. minorant) functions [19]. Let $t \in \mathbb{R}$. Function $x \mapsto \bar{b}(x, t)$ is said to be a tangent minorant to function $\pi$ at $t$ if

$$ (\forall x \in \mathbb{R}) \quad \pi(x) \ge \underline{b}(x; t) \quad \text{and} \quad \pi(t) = \bar{b}(t; t). \qquad (4) $$

Respectively, function $x \mapsto \bar{b}(x; t)$ is said to be a tangent majorant to function $\pi$ at $t$ if the following conditions hold:

$$ (\forall x \in \mathbb{R}) \quad \pi(x) \le \bar{b}(x; t) \quad \text{and} \quad \pi(t) = \bar{b}(t; t). \qquad (5) $$
In this section, we present how to construct tangent minorant (resp. majorant) functions for $\pi$, under structural assumptions on $\phi$. The approximations take the form of unnormalized Gaussian densities. From these, we deduce a simple strategy to derive rough lower and upper bounds for integral $\mathcal{I}$ in (1), when $f$ is the monomial given by (3).

In the following, we will make use of the following shorter notation, for the Gaussian probability density function (pdf):

$$ (\forall x \in \mathbb{R}) \quad g(x; \mu, \sigma) \equiv \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right), \qquad (6) $$

with mean $\mu \in \mathbb{R}$ and standard deviation $\sigma > \mathbb{R}^+$.

## 2.2. Gaussian tangent minorants

In order to build a tangent minorant with Gaussian shape for $\pi$, we rely on the following assumption:

**Assumption 1.** For every $t \in \mathbb{R}$, there exists $\beta(t) > 0$ such that

$$ (\forall x \in \mathbb{R}) \quad \phi(x) \le \phi(t) + \dot{\phi}(t)(x - t) + \frac{\beta(t)}{2}(x - t)^2, \qquad (7) $$

where $\dot{\phi}$ denotes the first order derivative of $\phi$.

Assumption 1 is satisfied by a wide class of differentiable functions. For instance, assume that $\phi$ is convex (i.e, $\pi$ is log-concave), twice differentiable on $\mathbb{R}$ with bounded second order derivative $\ddot{\phi}$ such that $\max_{x \in \mathbb{R}} \ddot{\phi}(x) = \beta_{\max} > 0$. Then, Assumption 1 holds with $\beta(t) = \beta_{\max}$ for every $t \in \mathbb{R}$ [21]. More generally, by the descent lemma, Assumption 1 is satisfied for any L-Lipschitz differentiable $\phi$ as soon as $\beta(t) \ge L$ for every $t \in \mathbb{R}$ [21]. Examples of functions $\phi$ fulfilling Assumption 1 are listed in Tab. 1. Other examples can be found for instance in [20, 18].

Let us emphasize that Assumption 1 is stable by summation: if $(\phi_j)_{1 \le j \le J}$ satisfy Assumption 1 for some $(\beta_j(\cdot))_{1 \le j \le J}$, then $\sum_{j=1}^J \phi_j$ satisfies Assumption 1 with $t \mapsto \beta(t) = \sum_{j=1}^J \beta_j(t)$. We will make use of this property to build our final approximation algorithm.

**Lemma 1.** Under Assumption 1, the minoration conditions (4) hold with

$$ (\forall (x, t) \in \mathbb{R}^2) \quad \underline{b}(x; t) = \underline{C}(t)g(x; \underline{u}(t), \underline{\sigma}(t)) \qquad (8) $$
<table><caption>Table 1: Examples of functions, defined on $\mathbb{R}$, satisfying Assumption 1. We provide the expression of $\phi$, of its derivative $\dot{\phi}$ and of the associated majorant curvature $\beta$. We assume that $\delta > 0$.</caption><thead><tr><th>Type</th><th>x ↦→ φ(x)</th><th>x ↦→ ˙φ(x)</th><th>t ↦→ β(t)</th></tr></thead><tbody><tr><td>Quadratic</td><td><sup>1</sup>&frasl;<sub>2</sub>x<sup>2</sup></td><td>x</td><td>1</td></tr><tr><td>Hyperbolic</td><td>(1 + x<sup>2</sup>/δ<sup>2</sup>)<sup>1/2</sup></td><td>(x/δ<sup>2</sup>)(1 + x<sup>2</sup>/δ<sup>2</sup>)<sup>-3/2</sup></td><td>δ<sup>-2</sup>(1 + t<sup>2</sup>/δ<sup>2</sup>)<sup>-3/2</sup></td></tr><tr><td>Huber</td><td><span style="font-size: 2.5em; vertical-align: middle;">{</span><table style="display: inline-table; vertical-align: middle; border: none; text-align: center; line-height: 1.2em;"><tbody><tr><td style="border: none; padding: 0 0.2em;">x<sup>2</sup>,</td><td style="border: none; padding: 0 0.2em;">if |x| &lt; δ,</td></tr><tr><td style="border: none; padding: 0 0.2em;">2δ|x| - δ<sup>2</sup>,</td><td style="border: none; padding: 0 0.2em;">otherwise</td></tr></tbody></table></td><td><span style="font-size: 2.5em; vertical-align: middle;">{</span><table style="display: inline-table; vertical-align: middle; border: none; text-align: center; line-height: 1.2em;"><tbody><tr><td style="border: none; padding: 0 0.2em;">2x,</td><td style="border: none; padding: 0 0.2em;">if |x| &lt; δ,</td></tr><tr><td style="border: none; padding: 0 0.2em;">2δsign(x),</td><td style="border: none; padding: 0 0.2em;">otherwise</td></tr></tbody></table></td><td><span style="font-size: 2.5em; vertical-align: middle;">{</span><table style="display: inline-table; vertical-align: middle; border: none; text-align: center; line-height: 1.2em;"><tbody><tr><td style="border: none; padding: 0 0.2em;">2,</td><td style="border: none; padding: 0 0.2em;">if |t| &lt; δ,</td></tr><tr><td style="border: none; padding: 0 0.2em;">2δ|t|<sup>-1</sup>,</td><td style="border: none; padding: 0 0.2em;">otherwise</td></tr></tbody></table></td></tr><tr><td>Logistic</td><td>log(1 + e<sup>x</sup>)</td><td>e<sup>x</sup>/(1 + e<sup>x</sup>)</td><td><span style="font-size: 2.5em; vertical-align: middle;">{</span><table style="display: inline-table; vertical-align: middle; border: none; text-align: center; line-height: 1.2em;"><tbody><tr><td style="border: none; padding: 0 0.2em;">1/4,</td><td style="border: none; padding: 0 0.2em;">if t = 0,</td></tr><tr><td style="border: none; padding: 0 0.2em;">1/(t(1 + e<sup>-t</sup>) - 1/2),</td><td style="border: none; padding: 0 0.2em;">elsewhere</td></tr></tbody></table></td></tr><tr><td>Cauchy</td><td>log(1 + x<sup>2</sup>/δ<sup>2</sup>)</td><td>2x(x<sup>2</sup> + δ<sup>2</sup>)<sup>-1</sup></td><td>2(t<sup>2</sup> + δ<sup>2</sup>)<sup>-1</sup></td></tr></tbody></table>

using the following variance, mean and scale parameters:

$$ (\forall t \in \mathbb{R}) \quad \begin{cases} \underline{\sigma}(t) = 1/\sqrt{\beta(t)}, \\ \underline{u}(t) = t - \underline{\sigma}(t)^2 \dot{\phi}(t), \\ \underline{C}(t) = \sqrt{2\pi\underline{\sigma}(t)^2} \exp\left(-\phi(t) + (\underline{\sigma}(t)^2/2)(\dot{\phi}(t))^2\right). \end{cases} \qquad (9) $$

105 Proof. Let $t \mapsto \beta(t)$ defined as in Assumption 1. We have:

$$ (\forall (x, t) \in \mathbb{R}^2) \quad \phi(t) + \dot{\phi}(t)(x - t) + \frac{\beta(t)}{2}(x - t)^2 \qquad (10) $$

$$ = \frac{\beta(t)}{2} x^2 + (\dot{\phi}(t) - \beta t)x + \phi(t) - t\dot{\phi}(t) + \frac{\beta(t)}{2}(t)^2 \qquad (11) $$

$$ = \frac{\beta(t)}{2} \left( x - \left(t - \frac{1}{\beta(t)} \dot{\phi}(t)\right)^2 \right) + \phi(t) - \frac{1}{2\beta(t)} (\dot{\phi}(t))^2 \qquad (12) $$

Thus, according to Assumption 1, for every $t \in \mathbb{R}$,

$$ (\forall x \in \mathbb{R}) \quad \phi(x) \le \frac{\beta(t)}{2}(x - u(t))^2 + \phi(t) - \frac{1}{2\beta(t)}(\dot{\phi}(t))^2 \qquad (13) $$

where we set

$$ (\forall t \in \mathbb{R}) \quad u(t) = t - \beta(t)^{-1} \dot{\phi}(t). \qquad (14) $$

Hence the result, using $\phi = -\log \pi$.

□
### 2.3. Gaussian tangent majorants

The majoration step requires an additional assumption on $\phi$, as follows:

**Assumption 2.** For every $t \in \mathbb{R}$, there exists $\nu(t) > 0$ such that

$$ (\forall x \in \mathbb{R}) \quad \phi(x) \ge \phi(t) + \dot{\phi}(t)(x-t) + \frac{\nu(t)}{2}(x-t)^2 \quad (15) $$

Assumption 2 is related to the notion of strong convexity (also sometimes called uniform convexity) [21]. It holds for instance for twice differentiable convex functions $\phi$ such that $\min_{x \in \mathbb{R}} \ddot{\phi}(x) = \nu > 0$, by setting $\nu(t) = \nu$ for every $t \in \mathbb{R}$. It also holds for $\phi : x \to \tilde{\phi}(x) + ax^2$, with any $a > 0$, and $\tilde{\phi}$ convex on $\mathbb{R}$. More generally, if $(\phi_j)_{1 \le j \le J}$ are convex functions, and for some $j \in \{1, \dots, J\}$, $\phi_j$ satisfies Assumption 2 with $\nu_j(\cdot)$, then $\sum_{j=1}^J \phi_j$ also does, with $t \mapsto \nu(t) = \nu_j(t)$.

**Lemma 2.** Under Assumption 2, the majoration conditions (5) are satisfied, as soon as

$$ (\forall (x,t) \in \mathbb{R}^2) \quad \bar{b}(x; t) = \bar{C}(t)g(x; \bar{u}(t), \bar{\sigma}(t)), \quad (16) $$

with variance, mean and scale parameters given by

$$ (\forall t \in \mathbb{R}) \quad \left\{ \begin{array}{l} \bar{\sigma}(t) = 1/\sqrt{\nu(t)}, \\ \bar{u}(t) = t - \bar{\sigma}(t)^2 \dot{\phi}(t), \\ \bar{C}(t) = \sqrt{2\pi\bar{\sigma}(t)^2} \exp\left(-\phi(t) + (\bar{\sigma}(t)^2/2)(\dot{\phi}(t))^2\right). \end{array} \right. \quad (17) $$

*Proof.* The proof follows similar structure than the one of Lemma 1, reverting inequalities.

### 2.4. First bounds and discussion

According to what precedes, we are able to construct, for every $t \in \mathbb{R}$, unnormalized Gaussian densities $\bar{b}(\cdot; t)$ and $\underline{b}(\cdot; t)$ that are tangent majorant and minorant functions for $\pi$ at $t$ i.e.,

$$ (\forall x \in \mathbb{R}) \quad \underline{b}(x; t) \le \pi(x) \le \bar{b}(x; t). \quad (18) $$

An example is displayed in Fig. 1, where the target function is displayed in blue line, the minorant functions are in black thin lines (left plot) and the majorant functions are in red thin lines (right plot).

We define the positive and negative parts of a function $f$ by

$$ \forall x \in \mathbb{R}, \quad f^{+}(x) = \max\{f(x), 0\}, \quad f^{-}(x) = \max\{-f(x), 0\}. \quad (19) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Example of function $\pi$ (blue thick line) and minorant/majorant functions (black/red thin lines) at various tangency points (black/red circles)</figcaption></figure>

With these definitions, we have the following identities:

$$ \forall x \in \mathbb{R}, \quad f(x) = f^{+}(x) - f^{-}(x), \quad |f(x)| = f^{+}(x) + f^{-}(x). \qquad (20) $$

Let us denote, for every $t \in \mathbb{R}$,

$$ \begin{aligned} \underline{\mathcal{I}}(t) &= \int_{\mathbb{R}} f^{+}(x) \underline{b}(x; t) dx - \int_{\mathbb{R}} f^{-}(x) \bar{b}(x; t) dx \\ &= \underline{C}(t) \int_{\mathbb{R}} f^{+}(x) g(x; \underline{u}(t), \underline{\sigma}(t)) dx - \bar{C}(t) \int_{\mathbb{R}} f^{-}(x) g(x; \bar{u}(t), \bar{\sigma}(t)) dx, \end{aligned} \qquad (21) $$

and

$$ \begin{aligned} \bar{\mathcal{I}}(t) &= \int_{\mathbb{R}} f^{+}(x) \bar{b}(x; t) dx - \int_{\mathbb{R}} f^{-}(x) \underline{b}(x; t) dx \\ &= \bar{C}(t) \int_{\mathbb{R}} f^{+}(x) g(x; \bar{u}(t), \bar{\sigma}(t)) dx - \underline{C}(t) \int_{\mathbb{R}} f^{-}(x) g(x; \underline{u}(t), \underline{\sigma}(t)) dx. \end{aligned} \qquad (22) $$

130

We deduce from (1), and (18) the following bounds,

$$ (\forall t \in \mathbb{R}) \quad \underline{\mathcal{I}}(t) \le \mathcal{I} \le \bar{\mathcal{I}}(t). \qquad (23) $$

Hence, the following inequality holds:

$$ \max_{t \in \mathbb{R}} \underline{\mathcal{I}}(t) \le \mathcal{I} \le \min_{t \in \mathbb{R}} \bar{\mathcal{I}}(t). \qquad (24) $$

135

Assuming that the test function $f$ reads as in (3), the terms $\underline{\mathcal{I}}(t)$ and $\bar{\mathcal{I}}(t)$ for every $t \in \mathbb{R}$, involve integrals that can be read as moments of (possibly truncated) Gaussian distributions. Such moments are well-defined, and can be computed easily, with high precision, based on recursive relations relying on the error function erf. These expressions are recalled in Appendix 7.
Still, computing the optimal bounds in (24) requires solving two optimization problems over the set $t \in \mathbb{R}$, whose resolution may not be straightforward. Relying on the preliminary results established in Lemmas 1 and 2, we propose in the remainder of the paper an iterative approximation strategy to refine the bounds, featuring two main ingredients: (i) a recursive and adaptive selection of a finite set of representative tangency points, and (ii) the computation of integral bounds using piecewise-Gaussian upper and lower approximations of $\pi$, delimited by the retained tangency points.

Thanks to this approach presented in Section 3, no explicit optimization is needed, and all numerical computations are straightforward, under the mild assumption that truncated Gaussian monomial moments are computable.

# 3. Sequential refinement of the bounds

## 3.1. Introduction

In Section 2, we have derived a method to construct tangent minorant and majorant functions for $\pi$, for any given value of tangency point $t \in \mathbb{R}$. As we can observe in Figure 1, the accuracy of such approximations highly depends on the distance to the tangency point. When $x = t$, the approximations are exact (because of the tangency property), and, as $x$ gets further from $t$, a degradation in terms of quality of approximation can be expected. We propose an improved approximation strategy, described in Figure 2: Given a list of $M \ge 2$ tangency points, we construct the associated tangent minorant (thin black lines) and majorant (thin red lines) functions for $\pi$ (blue line) at these points. Then, from this set of $M$ functions, improved minorant (thick blue line) and majorant (thick red line) approximations of $\pi$ are deduced. These envelopes are defined in a piecewise manner, as the upper (resp. lower) bound of the available tangent minorant (resp. majorant) functions. This idea is at the core of the proposed Algorithm 1. In order to improve again the performance of the method, we then propose an iterative selection of the tangency points, through our Algorithm 2. By construction, such approach guarantees an increasing precision of the bounds along its iterations, and, as we will show in the next section, it produces a sequence of bounds that asymptotically converge to the sought integral value.
<figure><img src="image_4.png" /><figcaption>Figure 2: Example of function $\pi$ (blue thick line), set of tangent minorant/majorant functions (black/red thin lines) at various tangency points, and associated piecewise minorant/majorant function (black/red thick line).</figcaption></figure>

## 3.2. Minorant/majorant envelopes

Let us introduce the notation that will be required below. We set $M \ge 2$ tangency points, and $\mathbf{T} = [t_1 \cdots t_m \cdots t_M] \in \mathbb{R}^M$ as the vector that contains the $M$ values of the tangency points, sorted in increasing order. The piecewise minorant function $\bar{C}(\cdot; \mathbf{T})$ and piecewise majorant function $\underline{C}(\cdot; \mathbf{T})$, associated to $\mathbf{T}$ are defined by

$$ (\forall x \in \mathbb{R}) \quad \underline{C}(x; \mathbf{T}) = \max_{m \in \{1, \dots, M\}} \underline{b}(x; t_m), \qquad (25) $$

with $\{\underline{b}(\cdot; t_m)\}_{1 \le m \le M}$ the set of $M$ tangent minorant functions with tangency points $(t_m)_{1 \le m \le M}$, and

$$ (\forall x \in \mathbb{R}) \quad \bar{C}(x; \mathbf{T}) = \min_{m \in \{1, \dots, M\}} \bar{b}(x; t_m), \qquad (26) $$

with $\{\bar{b}(\cdot; t_m)\}_{1 \le m \le M}$ the set of $M$ tangent majorant functions with tangency points $(t_m)_{1 \le m \le M}$. Therefore, by construction,

$$ (\forall x \in \mathbb{R}) \quad \underline{C}(x; \mathbf{T}) \le \pi(x) \le \bar{C}(x; \mathbf{T}). \qquad (27) $$

Note that, due to the tangency condition, and the optimality conditions (25)-(26), the function and the piecewise majorant and minorant functions take the same value at each tangent point, i.e.,

$$ (\forall m \in \{1, \dots, M\}) \quad \bar{C}(t_m; \mathbf{T}) = \underline{C}(t_m; \mathbf{T}) = \pi(t_m). \qquad (28) $$

The optimization problems (25) and (26) are known in the literature of combinatorial algorithms as the search for *upper envelope* and *lower envelope* of a collection of univariate functions.
Such problems can be solved in an efficient manner, using the recursive algorithm from [22]. We now describe this computational algorithm applied to our specific problem, i.e., with majorant/minorant functions that take the form of unnormalized Gaussian densities.

Let us focus on problem (25), which amounts to determining the upper envelope of the set of minorant functions. We introduce the short notation $\psi_{t_m}(x) \triangleq \underline{b}(x; t_m)$ the minorant built from the tangency point $t_m$, for every $m \in \{1, \dots, M\}$. The upper envelope is obtained by calling the recursive Algorithm 1 with inputs $T = \{t_1, \dots, t_M\}$ and the associated family $\{\psi_{t_m}\}_{1\le m \le M}$. Namely, the solution of (25) reads:

$$ \forall i \in \{1, \dots, N\}, \forall x \in (v_i, v_{i+1}), \quad \underline{C}(x; T) = \psi_{u_i}(x), \qquad (29) $$

where each $u_i \in T$ is a tangency point selected by the algorithm over the open set $(v_i, v_{i+1})$.
**Algorithm 1 Recursive Compound Minimization Algorithm**

**Require:** Vector **T** of size *L*, and functions {$\psi_{t_i}$}$_{1\le i\le L}$

**Ensure:** **u** of size *N*, and **v** of size *N* + 1

1: Partition **T** into **T**$^{(1)}$ and **T**$^{(2)}$ of sizes $L_1 \le \lceil \frac{L}{2} \rceil$ and $L_2 \le \lceil \frac{L}{2} \rceil$

2: **if** $L_1 > 1$ **then**

3: Call Algorithm 1 on **T**$^{(1)}$ and {$\psi_{t_i}$}$_{1\le i\le L_1}$, returning **v**$^{(1)}$ and **u**$^{(1)}$

4: **if** $L_2 > 1$ **then**

5: Call Algorithm 1 on **T**$^{(2)}$ and {$\psi_{t_i}$}$_{1\le i\le L_2}$, returning **v**$^{(2)}$ and **u**$^{(2)}$

6: Merge **v**$^{(1)}$ and **v**$^{(2)}$ into ordered list **v**' = [-$\infty$, $v'_2$, ..., $v'_{N'}$, $\infty$]

7: **for** *i* = 1 to *N'* **do**

8: Find $u_i^{*(1)} \in \mathbf{u}^{(1)}$ such that:

$$ (\forall x \in (v'_i, v'_{i+1}))(\forall j \in \mathbf{u}^{(1)}) \quad \psi_{u_i^{*(1)}}(x) \ge \psi_j(x) \qquad (30) $$

9: Find $u_i^{*(2)} \in \mathbf{u}^{(2)}$ such that:

$$ (\forall x \in (v'_i, v'_{i+1}))(\forall j \in \mathbf{u}^{(2)}) \quad \psi_{u_i^{*(2)}}(x) \ge \psi_j(x) \qquad (31) $$

10: Compute the intersection point(s) of $\psi_{u_i^{*(1)}}$ and $\psi_{u_i^{*(2)}}$ in ($v'_i, v'_{i+1}$)

11: Add the intersection point(s) to the list **v**'

12: Sort and deduplicate **v**' to get final **v** = [$v_1, v_2, \dots, v_{N+1}$]

13: **for** *i* = 1 to *N* **do**

14: Compare $\psi_{u_i^{*(1)}}$ and $\psi_{u_i^{*(2)}}$ on ($v_i, v_{i+1}$)

15: Let $u_i \in \{1, \dots, M\}$ be the index of the dominant function on this interval

16: Add $u_i$ to **u**

185

The lower envelope of the set of majorant functions, defined as the solution of Eq. (26) can be obtained similarly as in Alg. 1. More precisely, the inequalities order should be reverted in both Eqs. (30) and (31) in the step 5 of Alg. 1, while in step 7, the lowest function in each interval must be selected.
### 3.3. Integral approximation

Once the upper and lower envelopes in (27) have been obtained, using the same reasoning that we use to find the equation 23, we are now able to express new bounds for the integral $\mathcal{I}$ as

$$ \int_{\mathbb{R}} f^{+}(x) \underline{C}(x; \mathbf{T}) dx - \int_{\mathbb{R}} f^{-}(x) \overline{C}(x; \mathbf{T}) dx \leq \mathcal{I} \leq \int_{\mathbb{R}} f^{+}(x) \overline{C}(x; \mathbf{T}) dx - \int_{\mathbb{R}} f^{-}(x) \underline{C}(x; \mathbf{T}) dx. \quad (32) $$

If $f$ reads as in (3), the integrals appearing in the bounds (32) involve products of monomials with piecewise unnormalized Gaussian densities, as defined in (25) and (26). These integrals can be computed explicitly piecewise, provided that the moments of the (possibly truncated) Gaussian distribution are available (see our Appendix 7).

### 3.4. Iterative tangency point selection

In order to refine again the precision of the bounds on $\mathcal{I}$, we now define an iterative, and adaptive, strategy for the selection of the tangency points, provided in Algorithm 2. At each iteration $n \in \mathbb{N}^*$ of Alg. 2, minorant and majorant functions $\underline{C}(x; \mathbf{T}^{(n)})$ and $\overline{C}(x; \mathbf{T}^{(n)})$ are constructed, for the current set of tangency points $\mathbf{T}^{(n)}$, by Alg. 1, and the lower and upper bounds, $\underline{\mathcal{I}}^{(n)}$ and $\overline{\mathcal{I}}^{(n)}$, for $\mathcal{I}$ are computed following the framework outlined in Eq. (32):

$$ \forall n \in \mathbb{N}^*, \quad \overline{\mathcal{I}}^{(n)} = \int_{\mathbb{R}} f^{+}(x) \overline{C}(x; \mathbf{T}^{(n)}) dx - \int_{\mathbb{R}} f^{-}(x) \underline{C}(x; \mathbf{T}^{(n)}) dx, \quad (33) $$

$$ \forall n \in \mathbb{N}^*, \quad \underline{\mathcal{I}}^{(n)} = \int_{\mathbb{R}} f^{+}(x) \underline{C}(x; \mathbf{T}^{(n)}) dx - \int_{\mathbb{R}} f^{-}(x) \overline{C}(x; \mathbf{T}^{(n)}). \quad (34) $$

The maximization of the difference between the integral bounds provided by the majorant and the minorant functions in all intervals between consecutive tangent points allows to determine a new tangency point, that yields the new set $\mathbf{T}^{(n+1)}$.
**Algorithm 2 Adaptive Tangency Point Selection**

**Require:** $\mathbf{T}^{(1)} = [t_1^{(1)}, \dots, t_{M_1}^{(1)}] \in \mathbb{R}^{M_1}$ with $M_1 \ge 2$, sorted increasingly; precision value $\varepsilon > 0$

**Ensure:** $\underline{\mathcal{I}}$ and $\bar{\mathcal{I}}$, lower and upper bounds on $\mathcal{I}$

1: **for** $n = 1, 2, \dots$ **do**

2: Construct $\bar{C}(\cdot; \mathbf{T}^{(n)})$ and $\underline{C}(\cdot; \mathbf{T}^{(n)})$ using Alg. 1

3: Define subintervals $S_1^{(n)} = (-\infty, t_1^{(n)}]$, $S_2^{(n)} = [t_1^{(n)}, t_2^{(n)}]$, ..., $S_{M_n+1}^{(n)} = [t_{M_n}^{(n)}, +\infty)$

4: **for** $i = 1$ to $M_n + 1$ **do**

5: Compute:

$$ \left\{ \begin{aligned} \underline{\mathcal{I}}_i^{(n)} &= \int_{S_i^{(n)}} f^+(x) \underline{C}(x; \mathbf{T}^{(n)}) dx - \int_{S_i^{(n)}} f^-(x) \bar{C}(x; \mathbf{T}^{(n)}) dx \\ \bar{\mathcal{I}}_i^{(n)} &= \int_{S_i^{(n)}} f^+(x) \bar{C}(x; \mathbf{T}^{(n)}) dx - \int_{S_i^{(n)}} f^-(x) \underline{C}(x; \mathbf{T}^{(n)}) dx \end{aligned} \right. \qquad (35) $$

6: Compute bounds:

$$ \left\{ \begin{aligned} \underline{\mathcal{I}}^{(n)} &= \sum_{i=1}^{M_n+1} \underline{\mathcal{I}}_i^{(n)} \\ \bar{\mathcal{I}}^{(n)} &= \sum_{i=1}^{M_n+1} \bar{\mathcal{I}}_i^{(n)} \end{aligned} \right. \qquad (36) $$

7: **if** $\bar{\mathcal{I}}^{(n)} - \underline{\mathcal{I}}^{(n)} \le \varepsilon$ **then**

8: **return** $(\underline{\mathcal{I}}^{(n)}, \bar{\mathcal{I}}^{(n)})$

9: **else**

10: **for** $i = 1$ to $M_n + 1$ **do**

11: Compute:

$$ G_i^{(n)} = \bar{\mathcal{I}}_i^{(n)} - \underline{\mathcal{I}}_i^{(n)} \qquad (37) $$

12: Select $\hat{i}$:

$$ \hat{i} = \operatorname*{argmax}_{i \in \{1, \dots, M_n+1\}} G_i^{(n)} \qquad (38) $$

13: Choose new point $\hat{t} \in S_{\hat{i}}^{(n)}$

14: Update:

$$ \mathbf{T}^{(n+1)} = [t_1^{(n)}, \dots, t_{\hat{i}-1}^{(n)}, \hat{t}, t_{\hat{i}}^{(n)}, \dots, t_{M_n}^{(n)}] \qquad (39) $$
A key feature in our method is the addition of a new tangency points in the intervals where the approximation gap is maximum (step 13 of Alg. 2). As we will show in our convergence analysis in Sec. 4, the implementation of step 13, coupled with the definition of $\mathbf{T}^{(1)}$, has an impact in the theoretical guarantees of the proposed approach. We provide, in Sec. 4.3, a detailed practical solution to select new tangency points, that allows to met the convergence requirements and to reach excellent practical performance.

205

210

# 4. Convergence Analysis

We now provide our theoretical analysis of Algorithm 2. The main idea of our analysis relies on the control of the following gap between the upper and lower envelopes, defined, for every $n \in \mathbb{N}^*$ and $x \in \mathbb{R}$, as:

$$h^{(n)}(x) = \bar{C}(x; \mathbf{T}^{(n)}) - \underline{C}(x; \mathbf{T}^{(n)}), \qquad (40)$$

where $\bar{C}(x; \mathbf{T}^{(n)})$ and $\underline{C}(x; \mathbf{T}^{(n)})$ are defined in Eqs. 26 and 25, respectively.

215

We start by establishing two lemmas. Lemma 3 establishes fundamental properties of the sequence $(h^{(n)}(x))_{n \in \mathbb{N}^*}$. These properties are crucial for constructing, in Lemma 4, Gaussian bounds on the products between $f^+$ (resp. $f^-$) and $h^{(n)}$, and then deduce an explicit control of the integrals of such products, outside compact sets. Our convergence result is finally presented in Theorem 1. Under the assumption of dense set of tangency points (see Sec. 4.3 for practical strategy to meet the assumption), we show that, for $n$ large enough, the absolute difference between $\bar{I}^{(n)}$ and $\underline{I}^{(n)}$ can be made arbitrary small, and we hence deduce our convergence proof.

220

We remind the reader that the analysis is made under the assumption of a monomial $f$ (i.e., (3) holds). This assumption is discussed at the end of the section.

## 4.1. Preliminary lemmas

225

**Lemma 3.** Let $(h^{(n)}(x))_{n \in \mathbb{N}^*}$ defined in (40). The following properties hold:

(i) For every $x \in \mathbb{R}$, the sequence $(h^{(n)}(x))_{n \in \mathbb{N}^*}$ is positive-valued and monotonically decreasing.

(ii) For every $n \in \mathbb{N}^*$, $h^{(n)}$ is continuous on $\mathbb{R}$.
Proof. (i) Let $n \in \mathbb{N}^*$. We have $\mathbf{T}^{(n)} \subset \mathbf{T}^{(n+1)}$, i.e., $\{t_m^{(n)}\}_{1\le m \le T_n} \subset \{t_m^{(n+1)}\}_{1\le m \le T_{n+1}}$. Then, by definition of the piecewise bounds, we have

$$ \forall x \in \mathbb{R}, \quad \pi(x) \le \bar{C}(x; \mathbf{T}^{(n+1)}) \le \bar{C}(x; \mathbf{T}^{(n)}), \qquad (41) $$

and

$$ \forall x \in \mathbb{R}, \quad \underline{C}(x; \mathbf{T}^{(n)}) \le \underline{C}(x; \mathbf{T}^{(n+1)}) \le \pi(x). \qquad (42) $$

Therefore,

$$ 0 \le h^{(n+1)}(x) \le h^{(n)}(x). \qquad (43) $$

Hence, we have the result.

(ii) We have that for every $n \in \mathbb{N}^*$, $\bar{C}(\cdot; \mathbf{T}^{(n)})$ and $\underline{C}(\cdot; \mathbf{T}^{(n)})$ are continuous in $\mathbb{R}$. Since they are respectively the max and the min of Gaussian functions, and we know that the max and the min are operators that maintain continuity.

Hence we have the result.

<figure><img src="image_5.png" /></figure>

**Corollary 1.** Let, for every $n \in \mathbb{N}$,

$$ G^{(n)} := \bar{\mathcal{I}}^{(n)} - \underline{\mathcal{I}}^{(n)} = \int_{\mathbb{R}} |f(x)| h^{(n)}(x) dx \qquad (44) $$

where $\bar{\mathcal{I}}^{(n)}$ and $\underline{\mathcal{I}}^{(n)}$ are defined in Eqs. 33 and 34, respectively. The following properties hold:

(i) The sequence $(G^{(n)})_{n \in \mathbb{N}^*}$ is positive-valued and monotonically decreasing.

(ii) There exists a function $h: \mathbb{R} \to \mathbb{R}^+$ such that the sequence $(h^{(n)})_{n \in \mathbb{N}^*}$, defined as in Eq. 40, converges pointwise to $h$.

Proof. (i) Let $n \in \mathbb{N}^*$. By Lemma 3,

$$ 0 \le \int_{\mathbb{R}} |f(x)| \left( \bar{C}(x; \mathbf{T}^{(n+1)}) - \underline{C}(x; \mathbf{T}^{(n+1)}) \right) dx \le \int_{\mathbb{R}} |f(x)| \left( \bar{C}(x; \mathbf{T}^{(n)}) - \underline{C}(x; \mathbf{T}^{(n)}) \right) dx, \quad (45) $$

hence the result.

(ii) By Lemma 3, we know that, for every $x \in \mathbb{R}$, sequence $(h^{(n)}(x))_{n \in \mathbb{N}^*}$ is positive valued and monotonically decreasing. Consequently, there exists $h: \mathbb{R} \to \mathbb{R}^+$ such that, for every $x \in \mathbb{R}$, $\lim_{n \to \infty} h^{(n)}(x) = h(x)$, hence the result.

<figure><img src="image_6.png" /></figure>

**Lemma 4.** Set $\mathbf{T}^{(1)} = \{t_1\}$ where $t_1 \in \mathbb{R}$. We have the following statements.
(i) For every $n \in \mathbb{N}^*$:

$$ (\forall x \in \mathbb{R}) \ (\forall s \in \mathbb{R}^+) \ f^+(x)h^{(n)}(x) \le \tilde{C}(s)g(x; \tilde{\mu}(s), \tilde{\sigma}(s)), \quad (46) $$

$$ (\forall x \in \mathbb{R}) \ (\forall s \in \mathbb{R}^-) \ f^-(x)h^{(n)}(x) \le \tilde{C}(s)g(x; \hat{\mu}(s), \hat{\sigma}(s)), \quad (47) $$

with $(h^{(n)})_{n \in \mathbb{N}^*}$ defined in (40) and

$$ (\forall s > 0) \quad \left\{ \begin{aligned} \tilde{\sigma}(s) &= \bar{\sigma}(t_1), \\ \tilde{\mu}(s) &= s - \tilde{\sigma}(s)^2 \left( \frac{s-\bar{\mu}(t_1)}{\bar{\sigma}(t_1)^2} - \frac{k}{s} \right), \\ \tilde{C}(s) &= \bar{C}(t_1) \frac{\tilde{\sigma}(s)}{\bar{\sigma}(t_1)} \exp \left( -\frac{(s-\bar{\mu}(t_1))^2}{2\bar{\sigma}(t_1)^2} + k \ln(s) + \frac{\tilde{\sigma}(s)^2}{2} \left( \frac{s-\bar{\mu}(t_1)}{\bar{\sigma}(t_1)^2} - \frac{k}{s} \right)^2 \right). \end{aligned} \right. \quad (48) $$

$$ (\forall s < 0) \quad \left\{ \begin{aligned} \hat{\sigma}(s) &= \bar{\sigma}(t_1), \\ \hat{\mu}(s) &= s - \hat{\sigma}(s)^2 \left( \frac{s-\bar{\mu}(t_1)}{\bar{\sigma}(t_1)^2} - \frac{k}{s} \right), \\ \hat{C}(s) &= \bar{C}(t_1) \frac{\hat{\sigma}(s)}{\bar{\sigma}(t_1)} \exp \left( -\frac{(-s+\bar{\mu}(t_1))^2}{2\bar{\sigma}(t_1)^2} + k \ln(-s) + \frac{\hat{\sigma}(s)^2}{2} \left( \frac{-s+\bar{\mu}(t_1)}{\bar{\sigma}(t_1)^2} + \frac{k}{s} \right)^2 \right). \end{aligned} \right. \quad (49) $$

where $\bar{\mu}(t_1)$, $\bar{\sigma}(t_1)$ and $\bar{C}(t_1)$ are defined in (2).

(ii) Let $\epsilon > 0$. Let $n \in \mathbb{N}^*$.

For all $s \in \mathbb{R}^+$, there exists

$$ \tilde{K}_{\epsilon}(s) = \left[ 0, \max \left( 0, \tilde{\mu}(s) + \tilde{\sigma}(s) \sqrt{\frac{2\tilde{C}(s)}{\epsilon}} \right) \right] \quad (50) $$

such that

$$ \int_{\mathbb{R} \setminus \tilde{K}_{\epsilon}(s)} f^{+}(x)h^{(n)}(x) dx \le \frac{\epsilon}{2}. \quad (51) $$

For all $s \in \mathbb{R}^-$, there exists

$$ \hat{K}_{\epsilon}(s) = \left[ \min \left( 0, \hat{\mu}(s) + \hat{\sigma}(s) \sqrt{\frac{2\hat{C}(s)}{\epsilon}} \right), 0 \right], \quad (52) $$

such that

$$ \int_{\mathbb{R} \setminus \hat{K}_{\epsilon}(s)} f^{-}(x)h^{(n)}(x) dx \le \frac{\epsilon}{2}. \quad (53) $$
Proof. (i) Let $n \in \mathbb{N}^*$. By Eq. (3), $(h^{(n)})_{n \in \mathbb{N}^*}$ is decreasing, hence

$$ (\forall x \in \mathbb{R}) \quad f^+(x)h^{(n)}(x) \le f^+(x)h^{(1)}(x) \\ \le f^+(x)(\bar{C}(x; \mathbf{T}^{(1)}) - \underline{C}(x; \mathbf{T}^{(1)})) \quad (54) \\ \le f^+(x)\bar{C}(t_1)g(x; \bar{\mu}(t_1), \bar{\sigma}(t_1)). $$

260

Using (3),

$$ (\forall x \in \mathbb{R}) \quad f^+(x)h^{(n)}(x) \le \max(0, \bar{C}(t_1)x^k g(x; \bar{\mu}(t_1), \bar{\sigma}(t_1))), \quad (55) $$

that is, for $x > 0$,

$$ \bar{C}(t_1) x^k g(x; \bar{\mu}(t_1), \bar{\sigma}(t_1)) = \bar{C}(t_1) \frac{1}{\sqrt{2\pi} \bar{\sigma}(t_1)} \exp(-\phi_+(x)), \quad (56) $$

with

$$ \phi_+(x) = \frac{1}{2\bar{\sigma}(t_1)^2} (x - \bar{\mu}(t_1))^2 - k \ln x \quad (x > 0). \quad (57) $$

Since

$$ \phi_+(x) = \frac{1}{\bar{\sigma}(t_1)^2} + \frac{k}{x^2} \ge \frac{1}{\bar{\sigma}(t_1)^2} \quad (x > 0), \quad (58) $$

265

$\phi_+$ satisfies Assumption 2 with the (constant) modulus $\nu \equiv 1/\bar{\sigma}(t_1)^2$. By Lemma 2, for every $s > 0$,

$$ (\forall x \in \mathbb{R}) \quad f^+(x)h^{(n)}(x) \le \tilde{C}(s)g(x; \tilde{\mu}(s), \tilde{\sigma}(s)), \quad (59) $$

with the parameters given in (48), namely

$$ \tilde{\sigma}(s) = \bar{\sigma}(t_1), \quad \tilde{\mu}(s) = s - \tilde{\sigma}(s)^2 \left( \frac{s - \bar{\mu}(t_1)}{\bar{\sigma}(t_1)^2} - \frac{k}{s} \right), \quad (60) $$

$$ \tilde{C}(s) = \bar{C}(t_1) \frac{\tilde{\sigma}(s)}{\bar{\sigma}(t_1)} \exp \left( - \frac{(s - \bar{\mu}(t_1))^2}{2\bar{\sigma}(t_1)^2} + k \ln s + \frac{\tilde{\sigma}(s)^2}{2} \left( \frac{s - \bar{\mu}(t_1)}{\bar{\sigma}(t_1)^2} - \frac{k}{s} \right)^2 \right). \quad (61) $$

For the negative part, parametrize the negative support by $y = -x > 0$. For $x \le 0$, $f^-(x) = (-x)^k$, and

$$ \bar{C}(t_1) (-x)^k g(x; \bar{\mu}(t_1), \bar{\sigma}(t_1)) = \bar{C}(t_1) y^k g(-y; \bar{\mu}(t_1), \bar{\sigma}(t_1)) = \bar{C}(t_1) \frac{1}{\sqrt{2\pi} \bar{\sigma}(t_1)} \exp(-\phi_-(y)), \quad (62) $$

270

with

$$ \phi_-(y) = \frac{1}{2\bar{\sigma}(t_1)^2} (y + \bar{\mu}(t_1))^2 - k \ln y \quad (y > 0). \quad (63) $$
Since

$$ \phi''(y) = \frac{1}{\bar{\sigma}(t_1)^2} + \frac{k}{y^2} \geq \frac{1}{\bar{\sigma}(t_1)^2} \quad (y > 0), \qquad (64) $$

$\phi_-$ also satisfies Assumption 2 with $\nu \equiv 1/\bar{\sigma}(t_1)^2$. By Lemma 2, for every $u > 0$,

$$ y^k g(-y; \bar{\mu}(t_1), \bar{\sigma}(t_1)) \leq \hat{C}(u) g(y; \hat{\mu}_+(u), \hat{\sigma}(u)), \qquad (65) $$

with

$$ \hat{\sigma}(u) = \bar{\sigma}(t_1), \quad \hat{\mu}_+(u) = u - \hat{\sigma}(u)^2 \phi'_-(u) = -\bar{\mu}(t_1) + \frac{\bar{\sigma}(t_1)^2 k}{u}, \qquad (66) $$

$$ \hat{C}(u) = \bar{C}(t_1) \frac{\hat{\sigma}(u)}{\bar{\sigma}(t_1)} \exp \left( - \frac{(u + \bar{\mu}(t_1))^2}{2\bar{\sigma}(t_1)^2} + k \ln u + \frac{\hat{\sigma}(u)^2}{2} \left( \frac{u + \bar{\mu}(t_1)}{\bar{\sigma}(t_1)^2} - \frac{k}{u} \right)^2 \right). \qquad (67) $$

Returning to $x = -y$ and using $g(-x; \mu, \sigma) = g(x; -\mu, \sigma)$, and reindexing by $s := -u < 0$ (so $-s > 0$), we obtain

$$ (\forall s < 0) \ (\forall x \in \mathbb{R}) \quad f^{-}(x)h^{(n)}(x) \leq \hat{C}(s)g(x; \hat{\mu}(s), \hat{\sigma}(s)), \qquad (68) $$

with the parameters given in (49), namely

$$ \hat{\sigma}(s) = \bar{\sigma}(t_1), \quad \hat{\mu}(s) = -\hat{\mu}_+(s) = s - \hat{\sigma}(s)^2 \left( \frac{-s + \bar{\mu}(t_1)}{\bar{\sigma}(t_1)^2} + \frac{k}{s} \right), \qquad (69) $$

$$ \hat{C}(s) = \bar{C}(t_1) \frac{\hat{\sigma}(s)}{\bar{\sigma}(t_1)} \exp \left( - \frac{(-s + \bar{\mu}(t_1))^2}{2\bar{\sigma}(t_1)^2} + k \ln(-s) + \frac{\hat{\sigma}(s)^2}{2} \left( \frac{-s + \bar{\mu}(t_1)}{\bar{\sigma}(t_1)^2} + \frac{k}{s} \right)^2 \right). \qquad (70) $$

This concludes the proof of (i).

(ii) Let $\epsilon > 0$ and $n \in \mathbb{N}^*$. In (i) we have proven that (46) and (47) hold for all $n \in \mathbb{N}^*$.

From the Bienaymé-Chebyshev inequality [23], that

$$ (\forall b > 0) \quad P(|X - E(X)| \geq b) \leq \frac{\mathrm{Var}(X)}{b^2} \qquad (71) $$

with $X$ a random variable. Let us set

$$ (\forall s > 0) \quad X \sim \mathcal{N}(\tilde{\mu}(s), \tilde{\sigma}(s)^2) \qquad (72) $$

Hence (71) yields:

$$ P(|X - \tilde{\mu}(s)| \geq b) \leq \frac{\tilde{\sigma}(s)^2}{b^2} \qquad (73) $$

Let us now choose $b > 0$ such that $\frac{\epsilon}{2\tilde{C}(s)} = \frac{\bar{\sigma}(s)^2}{b^2}$, i.e.

$$ b = \sqrt{\frac{2\tilde{\sigma}(s)^2\tilde{C}(s)}{\epsilon}} \qquad (74) $$
By substituting $b$ in (73), we have:

$$P(|X - \tilde{\mu}(s)| \geq b) \leq \frac{\epsilon}{2\tilde{C}(s)}. \qquad (75)$$

Therefore,

$$\int_{\mathbb{R} \setminus [\tilde{\mu}(s)-b, \tilde{\mu}(s)+b]} g(x; \tilde{\mu}(s), \tilde{\sigma}(s)) dx \leq \frac{\epsilon}{2\tilde{C}(s)}. \qquad (76)$$

We distinguish two cases, depending on the sign of $\tilde{\mu}(s) + b$.

First, assume that, $(\forall s > 0) \tilde{\mu}(s) + b \geq 0$. If $(\forall s > 0) \tilde{\mu}(s) - b \geq 0$, then,

$$[\tilde{\mu}(s) - b, \tilde{\mu}(s) + b] \subset [0, \tilde{\mu}(s) + b].$$

By the positivity of $g(\cdot; \tilde{\mu}(s), \tilde{\sigma}(s))$,

$$\forall s > 0, \quad \int_{\mathbb{R} \setminus [0, \tilde{\mu}(s)+b]} g(x; \tilde{\mu}(s), \tilde{\sigma}(s)) dx \leq \int_{\mathbb{R} \setminus [\tilde{\mu}(s)-b, \tilde{\mu}(s)+b]} g(x; \tilde{\mu}(s), \tilde{\sigma}(s)) dx \leq \frac{\epsilon}{2}. \qquad (77)$$

If $(\forall s > 0) \tilde{\mu}(s) - b \leq 0$, by definition of $f^+$,

$$\int_{\mathbb{R}^{-}} f^{+}(x)h^{(n)}(x) dx = 0.$$

This leads to

$$\int_{\mathbb{R} \setminus [\tilde{\mu}(s)-b, \tilde{\mu}(s)+b]} f^{+}(x)h^{(n)}(x) dx \leq 0 + \int_{\tilde{\mu}(s)+b}^{\infty} f^{+}(x)h^{(n)}(x) dx \leq \tilde{C}(s) \int_{\tilde{\mu}(s)+b}^{\infty} g(x; \tilde{\mu}(s), \tilde{\sigma}(s)) dx \leq \frac{\epsilon}{2}.$$

We now jump to the second case, assuming that $(\forall s > 0) \tilde{\mu}(s) + b \leq 0$. Here, we have

$$[\tilde{\mu}(s) - b, \tilde{\mu}(s) + b] \subset [\tilde{\mu}(s) - b, 0].$$

As before, we get

$$\forall s \in \mathbb{R}^{+}, \quad \int_{\mathbb{R} \setminus [\tilde{\mu}(s)-b, 0]} g(x; \tilde{\mu}(s), \tilde{\sigma}(s)) dx \leq \int_{\mathbb{R} \setminus [\tilde{\mu}(s)-b, \tilde{\mu}(s)+b]} g(x; \tilde{\mu}(s), \tilde{\sigma}(s)) dx \leq \frac{\epsilon}{2}. \qquad (78)$$

In a nutshell, in all cases, for all $s > 0$, we have proved that there exists $\tilde{K}_{\epsilon}(s) = [0, \max(0, \tilde{\mu}(s)+\tilde{\sigma}(s)\sqrt{\frac{2\tilde{C}(s)}{\epsilon}})]$ such that

$$\int_{\mathbb{R} \setminus \tilde{K}_{\epsilon}(s)} f^{+}(x)h^{(n)}(x) dx \leq \frac{\epsilon}{2}. \qquad (79)$$

The case where $\tilde{K}_{\epsilon}(s) = \{0\}$ is avoided by picking an $s > 0$ such that $\tilde{\mu}(s) > 0$ (such $s$ exists since $\lim_{s \to 0} \tilde{\mu}(s) = +\infty$). Hence, we have proven (51).
By the same reasoning as for $f^+$ and proceeding exactly as in the proof of (i)—performing the change of variables $y = -x > 0$, applying Lemma 2 to $\phi_-(y) = \frac{(y+\bar{\mu}(t_1))^2}{2\bar{\sigma}(t_1)^2} - k \ln y$, mapping back via $x = -y$ together with $g(-x; \mu, \sigma) = g(x; -\mu, \sigma)$, and reindexing the support point $u > 0$ as $s = -u < 0$—we obtain, for every $s < 0$, the pointwise bound $f^-(x)h^{(n)}(x) \le \hat{C}(s)g(x; \hat{\mu}(s), \hat{\sigma}(s))$ with the parameters of (49); choosing $b = \hat{\sigma}(s)\sqrt{2\hat{C}(s)/\epsilon}$ and noting that $f^-(x) = 0$ for $x > 0$, we set $\hat{K}_\epsilon(s) = [\min(0, \hat{\mu}(s) + b), 0]$ and obtain $\int_{\mathbb{R}\setminus\hat{K}_\epsilon(s)} f^-(x)h^{(n)}(x)dx \le \epsilon/2$, which proves (53).

□

## 4.2. Convergence theorem

We are now ready to state our main theorem for the convergence of Algorithm 2.

**Theorem 1** (Outside-compact control and convergence). Assume that $T = \bigcup_{n \in \mathbb{N}^*} \mathbf{T}^{(n)}$ is dense in $\mathbb{R}$, and $\mathbf{T}^{(1)} = \{t_1\}$ where $t_1 \in \mathbb{R}$. Define $s \mapsto \tilde{K}_\epsilon(s)$ and $s \mapsto \hat{K}_\epsilon(s)$ as in Lemma 4. Then, for every $n \in \mathbb{N}^*$ and every $\epsilon > 0$, the following outside-compact bounds hold:

$$ \int_{\mathbb{R}\setminus\tilde{K}_\epsilon(s)} f^+(x) h^{(n)}(x) dx \le \frac{\epsilon}{2}, \qquad (80) $$

$$ \int_{\mathbb{R}\setminus\hat{K}_\epsilon(s)} f^-(x) h^{(n)}(x) dx \le \frac{\epsilon}{2}. \qquad (81) $$

Moreover, there exist $s_0 > 0$, $s_1 < 0$ and $(\tilde{n}, \tilde{n}) \in \mathbb{N}^2$ such that, for all $n > \max(\tilde{n}, \tilde{n})$,

$$ |G^{(n)}| \le \epsilon (1 + \tilde{\mu}(s_0) + \hat{\mu}(s_1)) + \sqrt{\epsilon} \left( \sqrt{2\hat{\sigma}(s_1)^2 \hat{C}(s_1)} + \sqrt{2\tilde{\sigma}(s_0)^2 \tilde{C}(s_0)} \right). \qquad (82) $$

In particular, $\lim_{n \to \infty} G^{(n)} = 0$.

**Proof.** Let $\epsilon > 0$, and $n \in \mathbb{N}^*$. According to Lemma 4, we can define $(\forall s \in \mathbb{R}^+) \tilde{K}_\epsilon(s)$ and $(\forall s \in \mathbb{R}^-) \hat{K}_\epsilon(s)$ as in (51) and (53). As stated in the proof of Lemma 4, the case where $\tilde{K}_\epsilon(s) = \{0\}$ and $\hat{K}_\epsilon(s) = \{0\}$ does not arise, simply by choosing an $s_0 > 0$ such that $\tilde{\mu}(s_0) > 0$ and an $s_1 < 0$ such that $\hat{\mu}(s_1) < 0$. Hence the proof for (80) and (81).

Let us define, for all $s \ge 0$, $\tilde{G}^{(n)}_\epsilon(s) = \mathbf{T}^{(n)} \cap \tilde{K}_\epsilon(s)$. Hence, $\tilde{G}_\epsilon(s) := T \cap \tilde{K}_\epsilon(s) = \bigcup_{n \in \mathbb{N}} \tilde{G}^{(n)}_\epsilon(s)$.

By assumption, $\tilde{G}_\epsilon(s)$ is dense in $\tilde{K}_\epsilon(s)$, therefore the complementary set, $T \cap \tilde{K}_\epsilon(s) = \mathbb{R} \cap \tilde{K}_\epsilon(s) = \tilde{K}_\epsilon(s)$.
Function $x \mapsto f^+(x)h^{(n)}(x)$ is continuous in $\mathbb{R}$ (see Eq. (3)), for every $n$ in $\mathbb{N}^*$. Consequently, $\forall x \in \mathbb{R}, x \mapsto f^+(x)h^{(n)}(x)$ is uniformly continuous in $K_\epsilon(s)$ by Heine's theorem [24].

Thus, there exists $\lambda > 0$, such that $(\forall (x, y)) \in \tilde{K}_\epsilon(s)^2$,

320

$$|x - y| < \lambda \Rightarrow |f^+(x)h^{(n)}(x) - f^+(y)h^{(n)}(y)| < \epsilon \quad (83)$$

Again, by density of the set $\tilde{G}_\epsilon(s)$,

$$\exists \tilde{n} \in \mathbb{N}^* \text{ such that } (\forall n \ge \tilde{n}) \quad \forall x \in \text{int}(\tilde{K}_\epsilon(s)), \quad |t_n - x| < \lambda, \quad (84)$$

By construction, $t_n \in \mathbf{T}^{(n)} \subset T$, and by definition, $h^{(n)}(t_n) = 0$.

Let $n \ge \tilde{n}$, by (84) and (83),

$$\forall x \in \text{int}(\tilde{K}_\epsilon(s)), \quad f^+(x)h^{(n)}(x) < \epsilon. \quad (85)$$

As a consequence,

$$\int_{\tilde{K}_{\epsilon}(s)} f^{+}(x)h^{(n)}(x) dx \le \epsilon \left( \tilde{\mu}(s_0) + \sqrt{\frac{2\tilde{\sigma}(s_0)^2 \tilde{C}(s_0)}{\epsilon}} \right) \quad (86)$$

325

Since $\tilde{K}_\epsilon(s) = [0, \tilde{\mu}(s_0) + \sqrt{\frac{2\tilde{\sigma}(s_0)^2 \tilde{C}(s_0)}{\epsilon}}]$, the length of $\tilde{K}_\epsilon(s)$ is $\tilde{\mu}(s_0) + \sqrt{\frac{2\tilde{\sigma}(s_0)^2 \tilde{C}(s_0)}{\epsilon}}$.

By (80) and (86),

$$\int_{\mathbb{R}} f^{+}(x)h^{(n)}(x) dx \le \frac{\epsilon}{2} + \epsilon \left( \tilde{\mu}(s_0) + \sqrt{\frac{2\tilde{\sigma}(s_0)^2 \tilde{C}(s_0)}{\epsilon}} \right) \quad (87)$$

One can show by applying an analogous reasoning that:

$$\exists \hat{n} \in \mathbb{N}^* \text{ such that } \forall n > \hat{n}, \quad \int_{\mathbb{R}} f^{-}(x)h^{(n)}(x) dx \le \frac{\epsilon}{2} + \epsilon \left( \hat{\mu}(s_1) + \sqrt{\frac{2\hat{\sigma}(s_1)^2 \hat{C}(s_1)}{\epsilon}} \right). \quad (88)$$

We have, $|G^{(n)}| = \int_{\mathbb{R}} f^{-}(x)h^{(n)}(x) dx + \int_{\mathbb{R}} f^{+}(x)h^{(n)}(x) dx$. Therefore, (82) holds and we have:

$$\lim_{n \to \infty} G^{(n)} = 0, \quad (89)$$

which concludes the proof. □

330

**Remark 1.** Theorem 1 establishes the convergence of the bounds when $f$ is a monomial, i.e., $f$ takes the form (3). This result could actually be extended further: since any polynomial function is
a finite linear combination of monomials, the same approximation procedure applies term-by-term, yielding convergence of the integral bounds for any polynomial integrand $f$. Interestingly, in such case, the moment computations in Algorithm 2 can still be defined in a closed form, by decomposing the polynomial into a sum of monomials and using the recursive formulas from Appendix 7.

### 4.3. From the convergence theorem to a practical algorithm

*Scope of the theorem.* Theorem 1 guarantees that the gap $G^{(n)} = \bar{I}^{(n)} - \underline{I}^{(n)}$ tends to zero provided the accumulated tangency set $T = \bigcup_{n \in \mathbb{N}^*} \mathbf{T}^{(n)}$ is dense in $\mathbb{R}$. Algorithm 2 in itself does not, by design, enforce density, since repeatedly refining the current max-gap interval does not prevent other regions from being neglected. Hence, the theorem proves convergence only in the specific (stronger) case where $T$ is dense.

*Adapting Algorithm 2 with a candidate pool (no replacement).* We implement Step 13 of Algorithm 2 to control the list of tangency points. Before inserting the suggested point $\hat{t}$, associated with the index $\hat{t}$, we project it onto the nearest available element of a candidate pool $\mathcal{U} \subset [a, b]$, and then remove it from $\mathcal{U}$ (selection *without replacement*). Conceptually, if $\mathcal{U}$ is infinite and dense on the region of interest, this selection inherits the density required by Theorem 1. In practice, we work with a finite dyadic pool; we denote the remaining candidates at an iteration $n^*$ by

$$ \mathcal{U}^{(n+1)} = \mathcal{U}^{(n)} \setminus \{\hat{t}\}, $$

where $\mathcal{U}^{(0)} = \mathcal{U}$.

*Compact used to define the pool (quantile construction).* The theoretical compacts in Lemma 4 are overly conservative in practice. We therefore define $[a, b]$ via Gaussian quantiles of the upper envelope $\bar{b}(\cdot; t_1)$ around the initial point $t_1$ (see [25, 26]). Given $\varepsilon > 0$,

$$ \int_{\mathbb{R} \setminus [a, b]} \bar{b}(x; t_1) dx \leq \bar{C}(t_1) \varepsilon, $$ 
with $a = \text{norminv}(\varepsilon/2, m(t_1), s(t_1))$, $b = \text{norminv}(1 - \varepsilon/2, m(t_1), s(t_1))$. (90)

Because $\bar{b}$ has Gaussian tails, the discrepancy between $\int_{\mathbb{R} \setminus [a, b]} \bar{b}$ and $\int_{\mathbb{R} \setminus [a, b]} |f| \bar{b}$ is of the same small order when $\varepsilon$ is small; the quantile is thus an accurate approximation of the outside-compact integral of the upper envelope.
355 Dyadic pool construction. We build the candidate set $\mathcal{U}$ from a dyadic decomposition of $[a, b]$, as described in Algorithm 3.

<caption>**Algorithm 3 Dyadic Set Construction**</caption>

**Require:** $\ell > 0$ (density parameter), $a$ (start of the interval), $b$ (end of the interval)
**Ensure:** $\mathcal{U}$ (dyadic candidate set), $M$ (number of candidates)
1: $\ell \leftarrow [b] - [a]$
2: $d_1 \leftarrow \max(1, \lfloor \ell/\max(1, \ell) \rfloor)$
3: $d_2 \leftarrow \lfloor \log_2(d_1) \rfloor$
4: $\mathcal{U} \leftarrow \emptyset$
5: $x \leftarrow [a]$
6: **while** $x \le \lfloor b \rfloor$ **do**
7: &nbsp;&nbsp;&nbsp;&nbsp;Append $x$ to $\mathcal{U}$
8: &nbsp;&nbsp;&nbsp;&nbsp;$x \leftarrow x + 2^{-d_2}$
9: $M \leftarrow |\mathcal{U}|$
10: **return** $\mathcal{U}, M$

360

**Remark 2.** The input parameter $\ell$ controls the dyadic density within $[a, b]$. The total number of candidates $M$ depends on both $\ell$ and the interval length. As $\ell$ increases, the subdivision becomes finer and $M$ increases. The construction ensures that all integers between $\lfloor a \rfloor$ and $\lfloor b \rfloor$ are included in $\mathcal{U}$, which guarantees that an integer tangency point $t_1$ (for instance, in all our experiments, $t_1 = 1$) belongs to the pool. The dyadic set is therefore aligned with the compact used for integration and offers reproducible control of its granularity.

365

**Step 13:** *pool-based insertion rule.* We replace the “choose a new point $\hat{t} \in S_t^{(n)}$” by the two-step rule described in Alg. 4 hereafter; the remaining candidates at iteration $n$ are denoted $\mathcal{U}^{(n)} \subseteq \mathcal{U}$ with $\mathcal{U}^{(0)} = \mathcal{U}$.
**Algorithm 4 Decision rule with dyadic nearest-candidate (without replacement)**

**Require:** Iteration index $n \in \mathbb{N}^*$. Current tangency set $\mathbf{T}^{(n)}$, gaps $\{G_i^{(n)}\}_{i=1}^{M_n+1}$, subintervals $\{S_i^{(n)}\}_{i=1}^{M_n+1}$, remaining dyadic candidates $\mathcal{U}^{(n)} \subseteq \mathcal{U}$, selected index $\hat{i}$ from (38)

**Ensure:** New point $\hat{t}$, updated sets $\mathbf{T}^{(n+1)}$, $\mathcal{U}^{(n+1)}$

if $\mathcal{U}^{(n)} = \emptyset$ then

return

▷ no candidate remains

if $\mathcal{U}^{(n)} \cap S_i^{(n)} = \emptyset$ then

$$ \hat{i} \leftarrow \arg\max\{ G_i^{(n)} : \mathcal{U}^{(n)} \cap S_i^{(n)} \neq \emptyset \} $$

if no such $i$ then

return

▷ no candidate remains

Compute the target location

$$ \hat{t}^* = \begin{cases} \frac{t_{\hat{i}-1}^{(n)} + t_{\hat{i}}^{(n)}}{2}, & 1 < \hat{i} < M_n + 1, \\ t_1^{(n)} - \frac{1}{M_n - 1} \sum_{m=1}^{M_n-1} (t_{m+1}^{(n)} - t_m^{(n)}), & \hat{i} = 1, \\ t_{M_n}^{(n)} + \frac{1}{M_n - 1} \sum_{m=1}^{M_n-1} (t_{m+1}^{(n)} - t_m^{(n)}), & \hat{i} = M_n + 1. \end{cases} $$

Pick the nearest available dyadic candidate in $S_i^{(n)}$ and remove it:

$$ \hat{t} \in \arg \min_{u \in \mathcal{U}^{(n)} \cap S_i^{(n)}} |u - \hat{t}^*|, \quad \mathcal{U}^{(n+1)} \leftarrow \mathcal{U}^{(n)} \setminus \{\hat{t}\}. $$

Insert $\hat{t}$ into $\mathbf{T}^{(n)}$ in sorted order to obtain $\mathbf{T}^{(n+1)}$.

*Synthesis of the proposed pipeline.*: Let us summarize below the steps of the complete proposed pipeline, satisfying the convergence Theorem 1:

1. Choose initial point $t_1 \in \mathbb{N}$.

2. Compute $[a, b]$ using (90).

370

3. Build dyadic candidate pool $\mathcal{U}$, with size $M > 0$, for a chosen density parameter $\ell > 0$, using Alg. 3.

4. Run Alg. 2, using Alg. 4 for new point selection at each call of step 13 in Alg. 2.
We terminate Algorithm 2 when the pool $\mathcal{U}$ is exhausted or as soon as

$$G^{(n)} = \bar{\mathcal{I}}^{(n)} - \underline{\mathcal{I}}^{(n)} \leq \tau,$$

for a user-chosen tolerance $\tau > 0$. The output are the integral bounds $(\underline{\mathcal{I}}^{(n)}, \bar{\mathcal{I}}^{(n)})$.

## 5. Application to importance sampler variance estimation

In this section, we provide the specific implementation of the proposed method to the problem of variance estimation of the importance sampler.

### 5.1. Importance sampling

Importance sampling (IS) is a general Monte Carlo technique for the approximation of a probability density function (pdf) of interest by a random measure composed of samples and weights [1]. Consider $\tilde{p}(x)$ a pdf in $x \in \mathbb{R}$, and $p(x)$ its unnormalized version, such $\tilde{p}(x) = \frac{p(x)}{\bar{\mathcal{I}}}$, with $\mathcal{Z}$ being the (generally) unknown normalizing constant. In an IS method, a set of $N$ samples, $\{x_n\}_{n=1}^N$, is drawn from a single proposal pdf, $q(x)$. Each sample, $x_n$, is then assigned an importance weight given by

$$w_n = \frac{p(x_n)}{q(x_n)}, \quad n = 1, \dots, N. \qquad (91)$$

The samples and the weights then form the random measure $\chi = \{x_n, w_n\}_{n=1}^N$ that approximates the measure of the target pdf as

$$\hat{p}_{\text{IS}}(x) = \frac{1}{N\hat{\mathcal{Z}}} \sum_{n=1}^{N} w_n \delta_{x_n}(x), \qquad (92)$$

where $\delta_{x_n}(x)$ is the unit delta measure concentrated at $x_n$. Given such approximation, it is direct to build approximations of moments of $p(x)$. For instance, $\hat{\mathcal{Z}} = \frac{1}{N} \sum_{n=1}^N w_n$ is an unbiased estimator of $\mathcal{Z} \triangleq \int_{\mathcal{X}} \pi(x)dx$. More generally, one can approximate $\mathcal{I} \triangleq \int_{\mathcal{X}} m(x)\tilde{p}(x)dx$ by the unbiased IS estimator

$$\hat{\mathcal{I}} = \frac{1}{N\hat{\mathcal{Z}}} \sum_{n=1}^{N} w_n m(x_n), \qquad (93)$$

assuming $\mathcal{Z}$ is known. Otherwise, if the target distribution is only known up to the normalizing constant (i.e., $\mathcal{Z}$ is unknown), one can use the self-normalized estimator

$$\tilde{\mathcal{I}} = \frac{1}{N\tilde{\mathcal{Z}}} \sum_{n=1}^{N} w_n m(x_n), \qquad (94)$$
where $\mathcal{Z}$ is approximated by its IS estimate. Under some mild assumptions, $\tilde{\mathcal{I}}$ is a consistent estimator of $\mathcal{I}$ [1, 27, 11].

395

## 5.2. Characterization of the importance sampling estimators

The variance of the estimator $\hat{\mathcal{I}}$ is directly related to the discrepancy between $\tilde{p}(x)|m(x)|$ and the proposal $q(x)$ [1, 28]. Specifically, the variance of the unbiased IS estimator of Eq. (93) is given by

400

$$ \mathrm{Var}_q(\hat{\mathcal{I}}) = \frac{1}{N\mathcal{Z}^2} \int_{\mathcal{X}} \frac{m^2(x)p^2(x)}{q(x)}dx - \frac{\mathcal{I}^2}{N}. \quad (95) $$

The computation of this variance is usually intractable, and then it should be approximated. An accurate approximation of it is crucial, to assess quantitatively the performance of the IS estimator, and to set adjustment strategies in its hyper-parameters (e.g., parameters of the proposal $q(x)$) to reduce this variance.

Using the approach introduced in Section 3, we propose to construct bounds, for each three integral terms in the variance expression, namely

405

Namely, we will apply Alg. 2, on the three above problems, that take the general form (1), with different roles for functions $\pi(x) = \exp(-\phi(x))$ and $f(x)$, identified as follows:

* Computation of the bounds in (96), using $\pi \equiv p$ and $f \equiv m$;

* Computation of the bounds in (97), using $\pi \equiv p$ and $f \equiv 1$;

* Computation of the bounds in (98), using $\pi \equiv p^2/q$ and $f \equiv m^2$.

410

Finally, combining the obtained bounds, the IS estimator variance will be approximated as:

$$ \frac{1}{N\bar{\mathcal{Z}}^2} \bar{\mathcal{J}} - \frac{\bar{\mathcal{I}}^2}{N} \le \mathrm{Var}_q(\hat{\mathcal{I}}) \le \frac{1}{N\bar{\mathcal{Z}}^2} \bar{\mathcal{J}} - \frac{\bar{\mathcal{I}}^2}{N}. \quad (99) $$

To illustrate and discuss the practical implementation, and assess the performance of this approach, we detail in the next section a practical example for the target function $\pi$. We also present a numerical comparison with a standard numerical integrator, as well as with the state-of-the-art method from [29] for integral approximation.
## 6. Numerical Results

### 6.1. Target function and proposal

We now consider an experimental study of the Bayesian logistic regression model, estimating the variance of the IS estimator as depicted in Sec. 5. Given $J \ge 1$ data, associated to $(y_j)_{1 \le j \le J} \in \{-1, 1\}$ class labels, and $(w_j)_{1 \le j \le J} \in \mathbb{R}$ features, the target function is constructed as the product of the Bayesian logistic likelihood with a Gaussian prior with zero mean and variance $s^2$ (with $s > 0$). The unnormarlized target function is

$$ (\forall x \in \mathbb{R}) \quad p(x) = A \exp(-x^2/(2s^2)) \times \prod_{j=1}^{J} \exp(-\log(1 + \exp(y_j(w_j x)))) \quad (100) $$

with $A > 0$. In all numerical experiments, we set $s = 1.2$, $J = 10$, and $A = 1$.

Following the notations in (1), we denote

$$ (\forall x \in \mathbb{R}) \quad \phi(x) = \frac{x^2}{2s^2} + \sum_{j=1}^{J} \log(1 + \exp(y_j w_j x)) - \log(A) \quad (101) $$

such that $\pi(x) = \exp(-\phi(x))$. Using [30, 31], we deduce that function $\phi$ in (101) satisfies Assump- tions 1 and 2 on $\mathbb{R}$, with the following expressions for $\beta(t)$ and $\nu(t)$:

$$ (\forall t \in \mathbb{R}) \quad \begin{cases} \beta(t) = \sum_{j=1}^{J} (y_j w_j)^2 \psi(y_j w_j t) + \frac{1}{s^2}, \\ \nu(t) = \frac{1}{s^2}, \end{cases} \quad (102) $$

where the function $\psi: \mathbb{R} \to (0, \frac{1}{4}]$ is defined as

$$ (\forall t \in \mathbb{R}) \quad \psi(t) = \begin{cases} \frac{1}{4} & \text{if } t = 0, \\ \frac{1}{t} \left( \frac{1}{1+\exp(-t)} - \frac{1}{2} \right) & \text{otherwise.} \end{cases} \quad (103) $$

We furthermore define, for the proposal of the IS sampler, a Gaussian proposal distribution $q(x)$ with mean equals 2 and variance $\theta > 0$. Figure 3 illustrates $p(x)$ and $p(x)^2/q(x)$ (using $\theta = 1.5$). The latter exhibits a sharper peak and faster decay compared to $p(x)$.
<figure><img src="image_7.png" /><figcaption>Figure 3: Densities p and p²/q as a function of x.</figcaption></figure>

We set $m(x) = x^2$. The final goal is to compute bounds on the variance estimator (95). We run Alg. 2 to construct bounds for the three intermediary quantities (96),(97), and (98), leading to the following bounds for the integral (95):

$$ \underline{V} = \frac{1}{N\bar{Z}^2} \underline{J} - \frac{\bar{I}^2}{N}, \quad \bar{V} = \frac{1}{N\bar{Z}^2} \bar{J} - \frac{\bar{I}^2}{N}. \qquad (104) $$

The bounds are compared to empirical variance estimation, obtained by running the IS simulator on $N_{runs}$ runs. Otherwise stated, we set $\theta = 1.5$, $N = 20$, and $N_{runs} = 10^6$.

## 6.2. Practical set-up

**Dyadic sets** As described in Sec. 4.3, the pool of available tangency points, for Alg. 2 is constructed by the procedure in Alg. 3, starting from an interval defined in (90). Setting $\epsilon = 10^{-6}$, we obtain $[a_p, b_p] = [-5.8655, 5.8744]$ when $\pi \equiv p$ (for (96) and (97)), and $[a_q, b_q] = [-5.9681, 4.0988]$ when $\pi \equiv p^2/q$ (for (98)). Moreover, in all experiments, we set $t_1 = 1$ and the density parameter for dyadic constructions $\ell = 10000$, which leads to $M_p = 6145$ (resp. $M_q = 5633$) candidate points in the dyadic sets obtained via Alg. 3

We structure the numerical study in three parts:

* **Validation of the obtained bounds.** We first assess in Sec. 6.2.1, the accuracy of the upper/lower envelopes produced by Algorithm 2, by comparing the midpoint

$$ \frac{\bar{X} + \underline{X}}{2} $$
against MATLAB's `integral()` for each target integral $\mathcal{X} \in \{\mathcal{Z}, \mathcal{J}, \mathcal{I}\}$. We report absolute and relative errors, together with the number of tangency points required.

* **Validation of the variance estimator.** Next, in Sec. 6.3, we compute the bounds on the variance of the Importance Sampling estimator for a range of values for the proposal variance $\theta$. This illustrates the validity of the bounds, compared to its empirical computation, and their applicative use to optimally set the proposal hyper-parameter $\theta$.

* **Comparison with state-of-the-art.** Finally, in Sec. 6.4, we compare our results with the envelope-based approximation of Evans et al. [29], in the case of the evaluation of bounds for the integrals $\mathcal{Z}$ and $\mathcal{I}$, using comparable computational budgets (same number of tangency points / evaluations), and discuss accuracy-cost trade-offs.

### 6.2.1. Numerical validation of bounds

This section evaluates the tightness and accuracy of the envelope bounds produced by Algorithm 2. First, we report, in Tables 2, 3, 4, the upper and lower bounds obtained when we stop as soon as a prescribed relative precision $\tau$ is reached. We also report the number of iterations $n_{stop}$ when reaching the stopping criterion. Note that, since $\mathbf{T}^{(1)}$ is reduced to the singleton $t_1$, we have $n_{stop} = \text{Card}(\mathbf{T}^{(n)})$, i.e., the number of evaluated tangency points. Finally, we display the difference between the estimated upper and lower bounds, as well as the absolute difference between the midpoint

$$ \frac{\bar{\mathcal{X}} + \underline{\mathcal{X}}}{2}, $$

with $\mathcal{X} \in \{\mathcal{I}, \mathcal{J}, \mathcal{Z}\}$, and the MATLAB's `integral()` result. As expected, the accuracy of the bounds increases when refining the stopping criterion, as this yields exploring more tangency points. The results are close to the ones obtained with the numerical integrator, the main difference is that our construction ensures that the bounds envelop the ground truth integral.

Second, we visualize on Figure 4 the convergence of the algorithm, until the relative-precision target $\tau = 10^{-4}$ is met.
<table><caption>Table 2: Obtained bounds for $\mathcal{I}$. The value computed with MATLAB's `integral()` is $2.641 \times 10^{-3}$.</caption><thead><tr><th colspan="6">Integral <i>I</i></th></tr><tr><th>&tau;</th><th>n<sub>stop</sub></th><th><span style="text-decoration: overline;">I</span></th><th><span style="text-decoration: overline;">I&#773;</span></th><th><span style="text-decoration: overline;">I&#773; - I&#773;</span></th><th>Diff. vs MATLAB</th></tr></thead><tbody><tr><td>1.000 &times; 10<sup>-2</sup></td><td>11</td><td>2.663 &times; 10<sup>-3</sup></td><td>2.638 &times; 10<sup>-3</sup></td><td>2.568 &times; 10<sup>-5</sup></td><td>9.271 &times; 10<sup>-6</sup></td></tr><tr><td>1.000 &times; 10<sup>-3</sup></td><td>34</td><td>2.643 &times; 10<sup>-3</sup></td><td>2.641 &times; 10<sup>-3</sup></td><td>2.605 &times; 10<sup>-6</sup></td><td>6.888 &times; 10<sup>-7</sup></td></tr><tr><td>1.000 &times; 10<sup>-4</sup></td><td>104</td><td>2.641 &times; 10<sup>-3</sup></td><td>2.641 &times; 10<sup>-3</sup></td><td>2.622 &times; 10<sup>-7</sup></td><td>7.237 &times; 10<sup>-8</sup></td></tr></tbody></table>

<table><caption>Table 3: Obtained bounds for $\mathcal{J}$. The value computed with MATLAB's `integral()` is $1.260 \times 10^{-5}$.</caption><thead><tr><th colspan="6">Integral <i>J</i></th></tr><tr><th>&tau;</th><th>n<sub>stop</sub></th><th><span style="text-decoration: overline;">J</span></th><th><span style="text-decoration: overline;">J&#773;</span></th><th><span style="text-decoration: overline;">J&#773; - J&#773;</span></th><th>Diff. vs MATLAB</th></tr></thead><tbody><tr><td>1.000 &times; 10<sup>-2</sup></td><td>12</td><td>1.271 &times; 10<sup>-5</sup></td><td>1.258 &times; 10<sup>-5</sup></td><td>1.254 &times; 10<sup>-7</sup></td><td>4.111 &times; 10<sup>-8</sup></td></tr><tr><td>1.000 &times; 10<sup>-3</sup></td><td>36</td><td>1.261 &times; 10<sup>-5</sup></td><td>1.260 &times; 10<sup>-5</sup></td><td>1.206 &times; 10<sup>-8</sup></td><td>4.042 &times; 10<sup>-9</sup></td></tr><tr><td>1.000 &times; 10<sup>-4</sup></td><td>112</td><td>1.261 &times; 10<sup>-5</sup></td><td>1.260 &times; 10<sup>-5</sup></td><td>1.256 &times; 10<sup>-9</sup></td><td>4.031 &times; 10<sup>-10</sup></td></tr></tbody></table>

<table><caption>Table 4: Obtained bounds for $\mathcal{Z}$. The value computed with MATLAB's `integral()` is $2.647 \times 10^{-3}$.</caption><thead><tr><th colspan="6">Integral <i>Z</i></th></tr><tr><th>&tau;</th><th>n<sub>stop</sub></th><th><span style="text-decoration: overline;">Z</span></th><th><span style="text-decoration: overline;">Z&#773;</span></th><th><span style="text-decoration: overline;">Z&#773; - Z&#773;</span></th><th>Diff. vs MATLAB</th></tr></thead><tbody><tr><td>1.000 &times; 10<sup>-2</sup></td><td>11</td><td>2.670 &times; 10<sup>-3</sup></td><td>2.645 &times; 10<sup>-3</sup></td><td>2.552 &times; 10<sup>-5</sup></td><td>1.027 &times; 10<sup>-5</sup></td></tr><tr><td>1.000 &times; 10<sup>-3</sup></td><td>31</td><td>2.650 &times; 10<sup>-3</sup></td><td>2.647 &times; 10<sup>-3</sup></td><td>2.602 &times; 10<sup>-6</sup></td><td>9.609 &times; 10<sup>-7</sup></td></tr><tr><td>1.000 &times; 10<sup>-4</sup></td><td>101</td><td>2.648 &times; 10<sup>-3</sup></td><td>2.647 &times; 10<sup>-3</sup></td><td>2.628 &times; 10<sup>-7</sup></td><td>8.935 &times; 10<sup>-8</sup></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 4: Evolution of the bounds along iterations $n$ of Alg. 2, until a relative precision of $\tau = 10^{-4}$ is reached, using Algorithm 2.</figcaption></figure>
Across Tables 2–4, Algorithm 4 consistently achieves the prescribed relative precisions with a modest number of tangency points (e.g., $\tau = 1.000 \times 10^{-4}$ with 104 points for $\mathcal{I}$, 112 for $\mathcal{J}$, and 101 for $\mathcal{Z}$), while keeping both the absolute gap between bounds and the deviation from MATLAB extremely small (down to the order of $10^{-8}$). The convergence plots in Figure 4 confirm the fast monotone contraction of the envelopes and their stabilization around the reference as iterations increase, illustrating that the adaptive placement of tangency points effectively targets the regions that most reduce the bound gap. We observe faster convergence for the lower bound. This is expected, as it relies on the accurate, and adaptive, curvature $\beta(t)$.

### 6.3. IS variance estimation

We present in Table 5 the obtained bounds on the variance of the IS method, computed as in (104), where the three intermediary bounds are computed with the proposed Alg. 2 with precision $\tau = 10^{-4}$. We compare these bounds with the empirical variance $V_e$ of the unbiased IS estimator. Namely, we evaluate (93), for $N_{runs}$ independent runs of IS. We use $N = 20$ samples, proposal parameter $\theta = 1.5$, and numerical intensive integration to compute the normalization constant of the target. This leads to $(\hat{I}_i)_{i \in [1, \dots, N_{runs}]}$, from which we calculate the variance $V_e$. We also calculate the theoretical value $\text{Var}_q(\hat{I})$ of the unbiased IS variance, whose expression is given in (95), computed with the numerical integrator from Matlab.

<table><caption>Table 5: Empirical variance $V_e$ of $(\hat{I}_i)_{i \in [1, \dots, N_{runs}]}$, for different $N_{runs}$, estimated bounds $(\underline{V}, \overline{V})$ by proposed approach, and theoretical value $\text{Var}_q(\hat{I})$ evaluated by intensive numerical integration in (95).</caption><thead><tr><th>V<sub>e</sub> (N<sub>runs</sub> = 10<sup>2</sup>)</th><th>V<sub>e</sub> (N<sub>runs</sub> = 10<sup>4</sup>)</th><th>V<sub>e</sub> (N<sub>runs</sub> = 10<sup>6</sup>)</th><th><u>V</u></th><th><u>V&#773;</u></th><th>Var<sub>q</sub>(I&#770;)</th></tr></thead><tbody><tr><td>4.277 &times; 10<sup>-2</sup></td><td>4.005 &times; 10<sup>-2</sup></td><td>4.018 &times; 10<sup>-2</sup></td><td>4.013 &times; 10<sup>-2</sup></td><td>4.018 &times; 10<sup>-2</sup></td><td>4.016 &times; 10<sup>-2</sup></td></tr></tbody></table>

As expected, as the number of runs $N_{runs}$ increases, $V_e$ converges towards the theoretical value (95) of the unbiased IS variance. The bounds computed with our method are valid and have good precision, while neither requiring any MC runs, nor intensive numerical integration.

Now, we run our approach to compute $(\underline{V}, \overline{V})$, for different values of the proposal variance $\theta$, dataset $J$ and prior hyper-parameter $s$. This experiment aims at illustrating a key motivation behind our work, i.e., to find the optimal proposal hyper-parameter, minimizing the variance of the estimator. The results are displayed in Figure 5(left), along with the empirical variance $V_e$.
evaluated for $N_{runs} = 10^6$ runs. We also display in Figure 5(right) the Mean Square Error (MSE) of the IS estimator to compute the second order moment. Namely, the MSE is defined as

$$ \text{MSE} = \frac{1}{N_s} \sum_{\ell=1}^{N_s} \left( \frac{\mathcal{I}}{\mathcal{Z}} - \hat{I}_{\ell} \right)^2, $$

where $\mathcal{I}$ and $\mathcal{Z}$ are computed by numerical integration. Both plots are the same, as the variance and the MSE of the estimator are identical. This illustrates that the IS estimator is unbiased. An important remark is that the range $\theta$, used in all plots, is set so as to guarantee the well definition of all involved integrals.

<figure><img src="image_9.png" /><figcaption>Figure 5: Obtained bounds (red, blue) and empirical variance $V_e$ with $N_{runs} = 10^6$ (left) and MSE on second order moment (right) for different values of $\theta$. Top: $J = 10$, $s = 1.2$, middle: $J = 10$, $s = 2$, bottom: $J = 2$, $s = 1.2$.</figcaption></figure>

We can see in Figure 5 that, depending on the values of $(J, s)$ (resp., dataset size, and prior
hyper-parameter), the plots have different aspects. In Figure 5(top), a minimum of the variance is reached, approximately at $\theta = 1.14$, which guides to tune the proposal in this case as $q = \mathcal{N}(2, 1.14)$. In Figure 5(middle), the best proposal is the smallest $\theta$ value, in the range of well-defined integrals. In contrast, in Figure 5(bottom), a large $\theta$ leads to a reduced variance. In all cases, the computed bounds are valid and accurate, and provide a direct feedback on the quality of the proposal hyper-parameter in to get a minimal variance for the estimator, without requiring expensive MC runs nor MSE computation (impossible without access to ground truth).

### 6.4. Comparison with state-of-the-art

Numerous algorithms are available to approximate integrals. The major part of them are constructed thanks to an interpolation method. This consists in approximating a function by an easier one where the two function are equal in a finite number of anchor points, for example one can use as an approximating family the well known Lagrange interpolator polynomial. However, only very few allow to compute exact bounds on integrals. We choose to compare our method to the one in [29] computing bounds relying on polynomial piece-wise envelopes. Their main result is reminded below.

**Lemma 5.** [29, Lem.2] Let $d \ge 0$. If $f^{(d)}$ is concave on $[a, b]$ then $(\forall x \in [a, b])$

$$ \sum_{k=0}^{d} \frac{f^{(k)}(a)}{(k+1)!} (b-a)^{k+1} + \frac{f^{(d)}(b) - f^{(d)}(a)}{b-a} \frac{(b-a)^{d+2}}{(d+2)!} \le \int_{a}^{b} f(x) dx \le \sum_{k=0}^{d+1} \frac{f^{(k)}(a)}{(k+1)!} (b-a)^{k+1}. \quad (105) $$

This inequality is reverse if $f^{(d)}$ is convex.

Starting from this lemma, the authors propose a piece-wise (i.e., compound) version of their approximation, by cutting the closed interval $[a, b]$ in $m$ sub-intervals, and derive a control of the precision of the obtained integral bounds, in [29, Lem.3]. In our practical implementation of this method, for a fair comparison, we use the same interval $[a, b]$ than the ones provided in Sec. 6.2 (obtained with (90)). We set a fixed number of iterations (i.e., cardinality of our pool of tangency points) for our method, namely $n = 3, 50$ or $100$. We set this same number for the compound points in Evans' method, and evaluate it with $d = 0, 1, 2$ degree values. We only provide the results for $\mathcal{I}$ and $\mathcal{Z}$, as Evans' method complexity was too large, for handling the approximation of integral $\mathcal{J}$.
<table><caption>Table 6: Estimation of Z - Absolute and relative errors for proposed method and Evans', for different compound-point budget.</caption><thead><tr><th colspan="3">Compound points = 3</th></tr><tr><th>Method</th><th>Absolute Error</th><th>Relative Error</th></tr></thead><tbody><tr><td>Ours</td><td>1.948 &times; 10<sup>-3</sup></td><td>4.306 &times; 10<sup>-1</sup></td></tr><tr><td>Evans (d=0)</td><td>1.351 &times; 10<sup>-2</sup></td><td>1.910 &times; 10<sup>0</sup></td></tr><tr><td>Evans (d=1)</td><td>1.403 &times; 10<sup>-2</sup></td><td>1.041 &times; 10<sup>0</sup></td></tr><tr><td>Evans (d=2)</td><td>2.772 &times; 10<sup>-1</sup></td><td>2.475 &times; 10<sup>1</sup></td></tr></tbody></table><table><thead><tr><th colspan="3">Compound points = 50</th></tr><tr><th>Method</th><th>Absolute Error</th><th>Relative Error</th></tr></thead><tbody><tr><td>Ours</td><td>1.180 &times; 10<sup>-6</sup></td><td>4.457 &times; 10<sup>-4</sup></td></tr><tr><td>Evans (d=0)</td><td>6.857 &times; 10<sup>-5</sup></td><td>2.557 &times; 10<sup>-2</sup></td></tr><tr><td>Evans (d=1)</td><td>1.182 &times; 10<sup>-5</sup></td><td>4.456 &times; 10<sup>-3</sup></td></tr><tr><td>Evans (d=2)</td><td>4.290 &times; 10<sup>-7</sup></td><td>1.620 &times; 10<sup>-4</sup></td></tr></tbody></table><table><thead><tr><th colspan="3">Compound points = 100</th></tr><tr><th>Method</th><th>Absolute Error</th><th>Relative Error</th></tr></thead><tbody><tr><td>Ours</td><td>2.816 &times; 10<sup>-7</sup></td><td>1.063 &times; 10<sup>-4</sup></td></tr><tr><td>Evans (d=0)</td><td>1.684 &times; 10<sup>-5</sup></td><td>6.341 &times; 10<sup>-3</sup></td></tr><tr><td>Evans (d=1)</td><td>1.433 &times; 10<sup>-6</sup></td><td>5.413 &times; 10<sup>-4</sup></td></tr><tr><td>Evans (d=2)</td><td>5.013 &times; 10<sup>-9</sup></td><td>1.894 &times; 10<sup>-6</sup></td></tr></tbody></table>
<table><caption>Table 7: Estimation of $\mathcal{I}$ - Absolute and relative errors for proposed method and Evans', for different compound-point budget.</caption><thead><tr><th colspan="3">Compound points = 3</th></tr><tr><th>Method</th><th>Absolute Error</th><th>Relative Error</th></tr></thead><tbody><tr><td>Ours</td><td>1.847 × 10<sup>-3</sup></td><td>4.353 × 10<sup>-1</sup></td></tr><tr><td>Evans (d=0)</td><td>1.042 × 10<sup>-2</sup></td><td>1.954 × 10<sup>0</sup></td></tr><tr><td>Evans (d=1)</td><td>8.886 × 10<sup>-3</sup></td><td>1.055 × 10<sup>0</sup></td></tr><tr><td>Evans (d=2)</td><td>3.671 × 10<sup>-1</sup></td><td>8.763 × 10<sup>0</sup></td></tr></tbody><thead><tr><th colspan="3">Compound points = 50</th></tr><tr><th>Method</th><th>Absolute Error</th><th>Relative Error</th></tr></thead><tbody><tr><td>Ours</td><td>1.207 × 10<sup>-6</sup></td><td>4.570 × 10<sup>-4</sup></td></tr><tr><td>Evans (d=0)</td><td>1.061 × 10<sup>-4</sup></td><td>3.936 × 10<sup>-2</sup></td></tr><tr><td>Evans (d=1)</td><td>2.291 × 10<sup>-5</sup></td><td>8.638 × 10<sup>-3</sup></td></tr><tr><td>Evans (d=2)</td><td>2.516 × 10<sup>-6</sup></td><td>9.522 × 10<sup>-4</sup></td></tr></tbody><thead><tr><th colspan="3">Compound points = 100</th></tr><tr><th>Method</th><th>Absolute Error</th><th>Relative Error</th></tr></thead><tbody><tr><td>Ours</td><td>2.916 × 10<sup>-7</sup></td><td>1.104 × 10<sup>-4</sup></td></tr><tr><td>Evans (d=0)</td><td>1.649 × 10<sup>-7</sup></td><td>6.242 × 10<sup>-5</sup></td></tr><tr><td>Evans (d=1)</td><td>2.805 × 10<sup>-6</sup></td><td>1.061 × 10<sup>-3</sup></td></tr><tr><td>Evans (d=2)</td><td>1.649 × 10<sup>-7</sup></td><td>6.242 × 10<sup>-5</sup></td></tr></tbody></table>

515

Tables 6 and 7 highlight the comparative performance of our method and the method proposed by Evans et al. Our method achieves better accuracy than Evans' method with order $d = 0$, when using a low number of compound points (3 points). This difference is particularly notable in the relative error, where Evans' method demonstrates lower accuracy due to the limited subdivision of the integration interval. However, when the number of compound points increases to 100, Evans' method closes the gap and provides results comparable to ours, though our method remains slightly more efficient in most cases.

As the order of differentiation increases to $d = 1$, Evans' method begins to show significant improvements, particularly for larger numbers of compound points. With 100 points, their method
520

achieves smaller absolute and relative errors than ours, indicating its ability to leverage higher-order
derivatives effectively to refine the bounds of the integral. This trend continues for $d = 2$, where
Evans' method significantly outperforms ours in both absolute and relative errors, particularly
when the interval is finely subdivided. These results emphasize the strength of Evans' method in
producing highly accurate bounds for higher-order approximations.

525

However, this improved accuracy comes at the cost of an increased computational complexity. Evans' method relies on higher-order derivatives and the determination of inflection points within the integration interval. Finding these inflection points involves solving root-finding problems, often using algorithms such as bisection or binary search. The complexity of Evans' method depends on two parameters: the number of compound points and the number of differentiations. In contrast, our method avoids high order derivative-based computations, offering an alternative approach to achieving accurate results. For this reason, we could not apply Evans' method to $\mathcal{J}$, as differentiating the integrand twice was overly cumbersome.

530

# 7. Conclusion

535

In this article, we introduced a new algorithm for computing accurate bounds on the moments of uni-dimensional distributions. An important application is the estimation of the variance of importance sampling estimator, allowing a facilitated tuning of proposal hyper-parameters. Our method achieves excellent precision, matching the performance of the competing approach from Evans et al. Furthermore, we provided theoretical guarantees of the iterative computation of the bounds, and a quality control on their values.

540

There are several avenues for future improvements. One potential direction is to explore other proposal families in the numerical application, such as Student's t laws. Another, more challenging, avenue is to generalize the method to approximate moments of high-dimensional distributions.

# Appendix A. Closed-form expressions for Gaussian and truncated-Gaussian monomial moments

545

In this appendix, we provide explicit formulas for the computation of monomial moments with respect to Gaussian and truncated Gaussian distributions, which will be essential tools to practically implement the proposed approach.
## A.1. Moments of the Gaussian distribution

Let $X \sim \mathcal{N}(\mu, \sigma^2)$ and $k \in \mathbb{N}$. The k-th moment of $X$ can be computed explicitly:

* If $\mu = 0$, then odd moments vanish and the even moments are given by

$$ \mathbb{E}[X^{2k}] = \sigma^{2k}(2k - 1)!! = \sigma^{2k} \cdot \frac{(2k)!}{2^k k!} $$

* For general $\mu \in \mathbb{R}$, the k-th moment is

$$ \mathbb{E}[X^k] = \sum_{j=0}^{\lfloor k/2 \rfloor} \binom{k}{2j} \mu^{k-2j} \sigma^{2j} (2j-1)!! $$

550

see e.g. [32].

## A.2. Moments of the truncated Gaussian distribution

Let $X \sim \mathcal{N}(\mu, \sigma^2)$ truncated to the interval $[a, b]$. The k-th conditional moment is defined as

$$ m_k := \mathbb{E}[X^k | a \le X \le b]. $$

These moments satisfy the recursive formula [33]:

$$ m_k = (k-1)\sigma^2 m_{k-2} + \mu m_{k-1} - \sigma \cdot \frac{b^{k-1}\phi\left(\frac{b-\mu}{\sigma}\right) - a^{k-1}\phi\left(\frac{a-\mu}{\sigma}\right)}{\Phi\left(\frac{b-\mu}{\sigma}\right) - \Phi\left(\frac{a-\mu}{\sigma}\right)}, $$

for $k \ge 1$, with base cases $m_{-1} = 0, m_0 = 1$.

Hereabove, we denoted $(\forall x \in \mathbb{R}) \phi(x) = \frac{1}{\sqrt{2\pi}}e^{-x^2/2}$, the standard normal probability distribution function, and $(\forall x \in \mathbb{R}) \Phi(x) = \int_{-\infty}^{x} \phi(t) dt$, the standard normal cumulative distribution function.

555

In particular, we have ([33, 32]):

$$ \begin{aligned} m_1 &= \mu - \sigma \cdot \frac{\phi(\beta) - \phi(\alpha)}{\Phi(\beta) - \Phi(\alpha)}, \\ m_2 &= \mu^2 + \sigma^2 - \sigma \cdot \frac{(\mu + b)\phi(\beta) - (\mu + a)\phi(\alpha)}{\Phi(\beta) - \Phi(\alpha)}, \end{aligned} $$

where $\alpha = \frac{a-\mu}{\sigma}$ and $\beta = \frac{b-\mu}{\sigma}$. When $\alpha$ and $\beta$ are both large (in absolute value), it is recommended to use numerically stabilized variants involving the scaled complementary error function $erfcx$, as detailed in [32, 34].
560 Acknowledgments

S.M. and E.C. acknowledge support from the European Research Council under Starting Grant MAJORIS ERC-2019-STG-850925.

References

[1] C. P. Robert, G. Casella, Monte Carlo Statistical Methods, Springer, 2004.

565 [2] J. S. Liu, J. S. Liu, Monte Carlo strategies in scientific computing, Vol. 10, Springer, 2001.

[3] V. I. Krylov, A. H. Stroud, Approximate calculation of integrals, Courier Corporation, 2006.

[4] P. J. Davis, P. Rabinowitz, Methods of Numerical Integration, Academic Press, 1984.

[5] J.-P. Demailly, Analyse numérique et équations différentielles, Grenoble Sciences, EDP Sciences, Les Ulis, France, 2006.

570 [6] H. Brass, K. Petras, Quadrature theory: the theory of numerical integration on a compact interval, no. 178, American Mathematical Soc., 2011.

[7] V. Elvira, L. Martino, P. Closas, Importance gaussian quadrature, IEEE Transactions on Signal Processing 69 (2020) 474–488.

575 [8] I. The MathWorks, MATLAB Integral Function Documentation, https://fr.mathworks.com/help/matlab/ref/integral.html (2025).

[9] D. Kahaner, C. Moler, S. Nash, Numerical Methods and Software, Prentice Hall, 1989.

[10] A. B. Owen, Monte Carlo theory, methods and examples, https://artowen.su.domains/mc/, 2013.

[11] V. Elvira, L. Martino, Advances in importance sampling, arXiv preprint arXiv:2102.05407 (2021).

[12] E. Novak, K. Ritter, High dimensional integration of smooth functions over cubes, Numerische Mathematik 75 (1996) 79–97.
[13] M. Evans, T. Swartz, Approximating Integrals via Monte Carlo and Deterministic Methods, Vol. 20 of Oxford Statistical Science Series, Oxford University Press, Oxford, UK, 2000.

585

[14] L. Martino, V. Elvira, F. Louzada, Effective sample size for importance sampling based on discrepancy measures, Signal Processing 131 (2017) 386-401.

[15] V. Elvira, L. Martino, C. P. Robert, Rethinking the effective sample size, International Statistical Review 90 (3) (2022) 525-550.

[16] A. Owen, Y. Zhou, Safe and effective importance sampling, Journal of the American Statistical Association 95 (449) (2000) 135-143.

590

[17] V. Elvira, L. Martino, D. Luengo, M. F. Bugallo, Generalized multiple importance sampling, Statistical Science 34 (1) (2019) 129-155.

[18] E. Chouzenoux, A. Jezierska, J.-C. Pesquet, H. Talbot, A Majorize-Minimize subspace approach for 12-10 image regularization, SIAM Journal on Imaging Sciences 6 (1) (2013) 563-591.
doi:10.1137/11085997X.
URL https://hal.science/hal-00789962

595

[19] E. Chouzenoux, J.-C. Pesquet, Optimization Methods for Signal Processing, in: C.Jutten, L. Duarte, S. Moussaoui (Eds.), Source Separation in Physical-Chemical Sensing, no. Chapter 2, IEEE Press, 2023.
URL https://inria.hal.science/hal-04250055

600

[20] Y. Sun, P. Babu, D. P. Palomar, Majorization-minimization algorithms in signal processing, communications, and machine learning, IEEE Transactions on Signal Processing 65 (3) (2017) 794-816. doi:10.1109/TSP.2016.2601299.

[21] H. Bauschke, P. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, 2nd Edition, Springer, New York, 2017.

605

[22] J.-D. Boissonnat, M. Yvinec, Algorithmic Geometry, Cambridge University Press, 1998.

[23] D. P. Bertsekas, J. N. Tsitsiklis, Introduction to Probability, Athena Scientific, 2008.

[24] W. Rudin, Principles of Mathematical Analysis, McGraw-Hill, 1976.
610

[25] M. Abramowitz, I. A. Stegun, Handbook of mathematical functions: with formulas, graphs, and mathematical tables, Dover Pub., number 55., 1972.

615

[26] The MathWorks, Inc., norminv: Inverse of the normal cumulative distribution function, mATLAB R2025a, Retrieved from https://fr.mathworks.com/help/stats/norminv.html#References (2025).

[27] A. Owen, Monte Carlo Theory, Methods and Examples, http://statweb.stanford.edu/~owen/mc/, 2013.

[28] H. Kahn, A. W. Marshall, Methods of reducing sample size in Monte Carlo computations, Journal of the Operations Research Society of America 1 (5) (1953) 263-278.

[29] M. Evans, T. Swartz, An algorithm for the approximation of integrals with exact error bounds (11 1997).

620

[30] G. Bouchard, Efficient bounds for the softmax function and applications to approximate inference in hybrid models, in: Proceedings of the Neural Information Processing Systems (NIPS 2008), Vol. 31, Vancouver, B.C., Canada, Dec. 2008.

625

[31] E. Chouzenoux, J.-B. Fest, SABRINA: A Stochastic Subspace Majorization-Minimization Algorithm, Journal of Optimization Theory and Applications 195 (2022) 919-952.
URL https://hal.science/hal-03793623

[32] J. F. de Cossio-Diaz, Moments of the univariate truncated normal distribution, https://arxiv.org/abs/1803.09541, march 26 (2018).

[33] E. Orjebin, A recursive formula for the moments of a truncated univariate normal distribution (September 26 2014).

630

[34] S. Chevillard, The functions erf and erfc computed with arbitrary precision and explicit error bounds, Information and Computation 216 (2012) 72-95. doi:10.1016/j.ic.2011.09.001.