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

HAL
open science

# Confidence Intervals for Randomized Quasi-Monte Carlo Estimators

Pierre L'Ecuyer, Marvin K Nakayama, Art B Owen, Bruno Tuffin

► To cite this version:

Pierre L'Ecuyer, Marvin K Nakayama, Art B Owen, Bruno Tuffin. Confidence Intervals for Randomized Quasi-Monte Carlo Estimators. Winter Simulation Conference, Dec 2023, San Antonio, United States. hal-04088085

HAL Id: hal-04088085

https://inria.hal.science/hal-04088085

Submitted on 3 May 2023

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
# Confidence Intervals for Randomized Quasi-Monte Carlo Estimators

Pierre L'Ecuyer, Marvin K Nakayama, Art B Owen, Bruno Tuffin

► To cite this version:

Pierre L'Ecuyer, Marvin K Nakayama, Art B Owen, Bruno Tuffin. Confidence Intervals for Randomized Quasi-Monte Carlo Estimators. WSC 2023 - Winter Simulation Conference, Dec 2023, San Antonio, United States. pp.1-18. hal-04088085v3

HAL Id: hal-04088085

https://inria.hal.science/hal-04088085v3

Submitted on 19 Feb 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Confidence Intervals for Randomized Quasi-Monte
Carlo Estimators

Pierre L'Ecuyer

Université de Montréal

CANADA

Marvin K. Nakayama

New Jersey Institute of Technology

USA

Art B. Owen

Stanford University

USA

Bruno Tuffin

Inria, Univ Rennes, CNRS, IRISA

FRANCE

ABSTRACT

Randomized Quasi-Monte Carlo (RQMC) methods provide unbiased estimators whose variance often converges at a faster rate than standard Monte Carlo as a function of the sample size. However, computing valid confidence intervals is challenging because the observations from a single randomization are dependent and the central limit theorem does not ordinarily apply. The natural solution is to replicate the RQMC process independently a small number of times to estimate the variance and use a standard confidence interval based on the normal or Student distribution. In this paper we investigate bootstrap methods for getting nonparametric confidence intervals for the mean using a modest number of replicates. Our main conclusion is that intervals based on the Student *t* distribution are more reliable than even the bootstrap *t* method on the integration problems arising from RQMC. In our examples, the individual RQMC estimates often had extreme kurtosis but always had mild skewness. That situation is very favorable for the standard confidence intervals. Simulation experiments cannot cover all the cases of interest, so we add some caveats in our conclusions.
# Confidence Intervals for Randomized Quasi-Monte Carlo Estimators

Pierre L'Ecuyer

Université de Montréal

CANADA

Marvin K. Nakayama

New Jersey Institute of Technology

USA

Art B. Owen

Stanford University

USA

Bruno Tuffin

Inria, Univ Rennes, CNRS, IRISA

FRANCE

NOTE: Corrections to the published version have been made in red here, on Page 7.

## ABSTRACT

Randomized Quasi-Monte Carlo (RQMC) methods provide unbiased estimators whose variance often converges at a faster rate than standard Monte Carlo as a function of the sample size. However, computing valid confidence intervals is challenging because the observations from a single randomization are dependent and the central limit theorem does not ordinarily apply. A natural solution is to replicate the RQMC process independently a small number of times to estimate the variance and use a standard confidence interval based on a normal or Student $t$ distribution. We investigate the standard Student $t$ approach and two bootstrap methods for getting nonparametric confidence intervals for the mean using a modest number of replicates. Our main conclusion is that intervals based on the Student $t$ distribution are more reliable than even the bootstrap $t$ method on the integration problems arising from RQMC.

## 1 INTRODUCTION

Integration is a fundamental task with a ubiquitous role in scientific computation. In very special settings a definite integral can be found in closed form, but in many others we must use numerical methods. When the integration problem is of very
# 1 INTRODUCTION

Integration is a fundamental task with a ubiquitous role in scientific computation. In very special settings a definite integral can be found in closed form, but in many others we must use numerical methods. When the integration problem is of very high dimension, then stochastic simulation (Monte Carlo or MC) methods dominate practice. There has been considerable progress in replacing those MC computations by deterministic quasi-Monte Carlo (QMC) methods, which can significantly improve accuracy. Those QMC methods do not support practical error estimation strategies. Randomized QMC (RQMC) methods have been developed that provide independent and unbiased estimates with accuracy equal to or even better than QMC provides. So far those replicates have primarily been used to estimate the variance of an RQMC estimate. Much less has been written about how to use such replicates to construct approximate confidence intervals. In this paper we explore various approaches to construct confidence intervals, including bootstrap methods, in RQMC. The motivation is as follows: any practical problem requiring that an answer must be estimated with accuracy further necessitates providing some evidence of the accuracy attained. For example, [40] mandates that risk studies include an uncertainty quantification.

Stochastic simulations get their randomness from a random number generator that “simulates” realizations of independent random variables uniformly distributed over the interval (0, 1) (denoted $\mathcal{U}(0, 1)$). Typically, each simulation run consumes a vector of $\mathbf{U}$ of $s$ successive random numbers from this sequence and outputs a simulated random variable $f(\mathbf{U})$, where $f$ represents all the computations made in the simulation. The aim is to estimate the expectation

$$ \mu = \mathbb{E}[f(\mathbf{U})] = \int_{[0,1]^s} f(\mathbf{u})d\mathbf{u}. \qquad (1) $$

QMC replaces the independent replicates of the $s$-dimensional random vector $\mathbf{U}$ by a (deterministic) low-discrepancy sequence of points $\mathbf{u}_0, \mathbf{u}_1, \mathbf{u}_2, \dots$ in $[0, 1)^s$. The first $n$ of these points cover the unit hypercube $[0, 1)^s$ far more evenly than random points would. When $n$ is fixed a priori, we talk of a low-discrepancy point set. QMC points are typically constructed for special values of $n$, such as powers of 2. The value of $\mu$ is then approximated by the average $\hat{\mu}_n = \frac{1}{n} \sum_{i=0}^{n-1} f(\mathbf{u}_i)$. Here $\mathbf{u}_i = (u_{i,1}, \dots, u_{i,s}) \in [0, 1)^s$.

Under certain assumptions on the integrand $f$, there are point sets and sequences for which the approximation error $|\hat{\mu} - \mu|$ converges as $\mathcal{O}(n^{-\alpha+\varepsilon})$ for any $\varepsilon > 0$, where $\alpha > 0$ depends on the smoothness assumptions made on $f$, which are expressed in terms of the integrals over the unit hypercube of its squared mixed partial derivatives with respect to the coordinates of $\mathbf{u}$ [5, 27, 22, 14]. As a special
high dimension, then stochastic simulation (Monte Carlo or MC) methods dominate practice. There has been considerable progress in replacing those MC computations by deterministic quasi-Monte Carlo (QMC) methods, which can significantly improve accuracy. But those QMC methods do not support practical error estimation strategies. Randomized QMC (RQMC) methods have been developed that provide independent and unbiased estimates with accuracy equal to or even better than QMC provides. So far those replicates have primarily been used to estimate the variance of an RQMC estimate. Much less has been written about how to use such replicates to construct approximate confidence intervals. In this paper we explore various approaches to construct confidence intervals, including bootstrap methods, in RQMC. The motivation is as follows: any practical problem requiring that an answer must be estimated with accuracy further necessitates providing some evidence of the accuracy attained; e.g., the work [40] mandates that risk studies include an uncertainty quantification.

Stochastic simulations get their randomness from a random number generator that “simulates” realizations of independent random variables uniformly distributed over the interval (0, 1) (denoted $\mathcal{U}(0,1)$). Typically, each simulation run consumes a vector of $\mathbf{U}$ of $d$ successive random numbers from this sequence and outputs a simulated random variable $f(\mathbf{U})$, where $f$ represents all the computations made in the simulation. The aim is to estimate the expectation

$$ \mu = \mathbb{E}[f(\mathbf{U})] = \int_{[0,1]^d} f(\mathbf{u})d\mathbf{u}. \qquad (1) $$

QMC replaces the independent replicates of the $d$-dimensional random vector $\mathbf{U}$ by a (deterministic) low-discrepancy sequence of points $\mathbf{u}_0, \mathbf{u}_1, \mathbf{u}_2, \dots$ in $[0, 1)^d$. The first $n$ of these points cover the unit hypercube $[0, 1)^d$ far more evenly than random points would. When $n$ is fixed a priori, we talk of a low-discrepancy point set. QMC points are typically constructed for special values of $n$, such as powers of 2. The value of $\mu$ is then approximated by the average $\hat{\mu}_n = \frac{1}{n} \sum_{i=0}^{n-1} f(\mathbf{u}_i)$. Here $\mathbf{u}_i = (u_{i,1}, \dots, u_{i,d}) \in [0, 1)^d$.

Under certain assumptions on the integrand $f$, there are point sets and sequences for which the QMC approximation error $|\hat{\mu} - \mu|$ converges as $\mathcal{O}(n^{-\alpha+\varepsilon})$ as $n \to \infty$ for any $\varepsilon > 0$, where $\alpha > 0$ depends on the smoothness assumptions made on $f$, which are expressed in terms of the integrals over the unit hypercube of its squared mixed partial derivatives with respect to the coordinates of $\mathbf{u}$ [5, 27, 22, 14]. As a special case, if $f$ is assumed to have bounded variation in the sense of Hardy and Krause, the rate is $\mathcal{O}(n^{-1}(\ln n)^d)$. Explicit constructions of these point sets and sequences are also available, mainly in the form of lattice point sets and digital nets and sequences [5, 19, 4]. One limitation, however, is that QMC fails to provide practical (easily computable) error bounds, so it is difficult to assess its accuracy
case, if $f$ is assumed to have bounded variation in the sense of Hardy and Krause, the rate is $\mathcal{O}(n^{-1}(\ln n)^s)$. Explicit constructions of these point sets and sequences are also available, mainly in the form of lattice point sets and digital nets and sequences [5, 19, 4]. One limitation, however, is that QMC fails to provide practical (easily computable) error bounds, so it is difficult to assess its accuracy [39, 26].

RQMC methods randomize the points in a way that each randomized point has the uniform distribution over $[0, 1)^s$ while the whole point set retains its high uniformity over this unit hypercube [18, 15, 29]. The RQMC estimator

$$ \hat{\mu}_{n,\text{rqmc}} = \frac{1}{n} \sum_{i=0}^{n-1} f(U_i) \qquad (2) $$

is unbiased for $\mu$, and under similar smoothness assumptions on $f$ as for QMC, its variance converges as $\mathcal{O}(n^{-2\alpha+\epsilon})$ for any $\epsilon > 0$, where $\alpha$ depends on the smoothness of $f$ and on the randomization method [5, 31]. But equation (2) averages dependent outputs, which complicates estimation of the variance of $\hat{\mu}_{n,\text{rqmc}}$ from a single randomization. A simple solution instead constructs $R > 1$ independent and identically distributed (iid) replicates of $\hat{\mu}_{n,\text{rqmc}}$ from independent randomizations, and computes the empirical variance of the $R$ iid estimators, which is an unbiased variance estimator.

A natural next step is to compute the standard confidence interval for $\mu$ based on a central limit theorem (CLT) for the average of the $R$ replicates. Unfortunately, the justification of this step is asymptotic in $R$ and, as explained below, we would prefer not to use large $R$. More specifically, let $y_r, r = 1, 2, \dots, R$, denote the $R$ iid realizations of $\hat{\mu}_{n,\text{rqmc}}$ in (2) and let $\bar{y}_{n,R}$ be their average. Typically, $\text{Var}[y_r] \approx Kn^{-2\alpha}$ for some $\alpha > 1/2$, so the variance of the global average is $\text{Var}[\bar{y}_{n,R}] \approx KR^{-1}n^{-2\alpha}$. For a given upper bound on the total number $Rn$ of simulation runs (total computing budget), this variance is minimized by taking $R$ as small as possible ($R = 1$). On the other hand, a small $R > 1$ gives a poor estimator of the variance, and a questionable CLT approximation. A compromise must be made if we want a reasonable estimate of the variance to assess the accuracy of our estimator of the mean. In practice, $R$ is often chosen in the range from 5 to 20. The global average $\bar{y}_{n,R}$ obeys a classical CLT when $R \to \infty$ for fixed $n$, but only in a few special cases when $n \to \infty$ for fixed $R$, so the standard intervals may fail to attain their nominal coverage. This is discussed in more detail in Sections 2 and 3.

In this paper, we explore empirically the distribution of $\hat{\mu}_{n,\text{rqmc}}$ for various integrands $f$, various choices of point sets and randomizations, and different values of $n$. We then compare the standard confidence interval with two others that have a nonparametric derivation. One of those is the bootstrap $t$ [7]. It has remarkably good asymptotic coverage [8] and exceptional empirical coverage in small samples [28]. The other is the percentile bootstrap method, possibly the best known
[39, 26].

RQMC methods randomize the points in a way that each randomized point has the uniform distribution over $[0, 1)^d$ while the whole point set retains its high uniformity over this unit hypercube [18, 15, 29]. The RQMC estimator

$$ \hat{\mu}_{n,\text{rqmc}} = \frac{1}{n} \sum_{i=0}^{n-1} f(U_i) \qquad (2) $$

is unbiased for $\mu$, and under similar smoothness assumptions on $f$ as for QMC, its variance converges as $\mathcal{O}(n^{-2\alpha+\varepsilon})$ for any $\varepsilon > 0$, where $\alpha$ depends on the smoothness of $f$ and on the randomization method [5, 31]. But equation (2) averages dependent outputs, which complicates estimation of the variance of $\hat{\mu}_{n,\text{rqmc}}$ from a single randomization. A simple solution instead constructs $R > 1$ independent and identically distributed (iid) replicates of $\hat{\mu}_{n,\text{rqmc}}$ from independent randomizations, and computes the empirical variance of the R iid estimators, which is an unbiased variance estimator.

A natural next step is to compute the standard confidence interval for $\mu$ based on a central limit theorem (CLT) for the average of the R replicates. Unfortunately, the justification of this step is asymptotic in $R$ and, as explained below, we would prefer not to use large $R$. More specifically, let $y_r, r = 1, 2, \dots, R$, denote the R iid realizations of $\hat{\mu}_{n,\text{rqmc}}$ in (2) and let $\bar{y}_{n,R}$ be their average. Typically, $\text{Var}[y_r] \approx Kn^{-2\alpha}$ for some $\alpha > 1/2$, so the variance of the global average is $\text{Var}[\bar{y}_{n,R}] \approx KR^{-1}n^{-2\alpha}$. For a given upper bound on the total number $Rn$ of simulation runs (total computing budget), this variance is minimized by taking $R$ as small as possible ($R = 1$). On the other hand, a small $R > 1$ gives a poor estimator of the variance, and a questionable CLT approximation. A compromise must be made if we want a reasonable estimate of the variance to assess the accuracy of our estimator of the mean. In practice, $R$ is often chosen in the range from 5 to 20. The global average $\bar{y}_{n,R}$ obeys a classical CLT when $R \to \infty$ for fixed $n$, but only in a few special cases when $n \to \infty$ for fixed $R$, so the standard intervals may fail to attain their nominal coverage. This is discussed in more detail in Sections 2 and 3.

In this paper, we explore empirically the distribution of $\hat{\mu}_{n,\text{rqmc}}$ for various integrands $f$, various choices of point sets and randomizations, and different values of $n$. We then compare the standard confidence interval with two others that have a nonparametric derivation. One of those is the bootstrap $t$ [7]. It has remarkably good asymptotic coverage [8] and exceptional empirical coverage in small samples [28]. The other is the percentile bootstrap method, possibly the best known bootstrap. For background on the bootstrap, see [6]. In this work we see some interesting cases where the distribution of $\hat{\mu}_{n,\text{rqmc}}$ is very far from normal and has a small but non-negligible probability mass which is quite far from the mean. Such
bootstrap. For background on the bootstrap, see [6]. In this work we see some interesting cases where the distribution of $\hat{\mu}_{n,\text{rqmc}}$ is very far from normal and has a small but non-negligible probability mass which is quite far from the mean. Such cases are very far from those where nonparametric confidence intervals have been studied.

The remainder is organized as follows. In Section 2, we recall the definitions of the point sets and randomizations considered in this paper, and the variance bounds that are known for the relevant combinations. In Section 3, we discuss several ways of computing confidence intervals RQMC estimators, including via a Student $t$ approximation and by applying the bootstrap $t$. In Section 4, we illustrate numerically the behavior of these methods for some functions $f$. Our conclusions are in Section 5.

## 2 POINT SETS AND RANDOMIZATION METHODS

Here we define the point sets and the randomizations considered in this paper. For each type of RQMC points, we define a short tag (acronym) to use later when reporting results.

Our first type of point set corresponds to a *lattice rule of rank 1*. It is defined by selecting the number of points, $n$, and a generating vector $\mathbf{a} = (a_1, \dots, a_s)$ of integers $a_j \in \{1, \dots, n-1\}$. Then the $i$th point $\mathbf{u}_i$ is defined by $\mathbf{u}_i = (i\mathbf{a} \bmod n)/n$, for $i = 0, \dots, n-1$. In our experiments, the vector $\mathbf{a}$ was selected using LatNet Builder [20, 19] with the $\mathcal{P}_\alpha$ criterion with $\alpha = 2$, product weights, and a weight of $\gamma_j = 2/(2+j)$ for coordinate $j$. See [19] for the details.

The standard randomization for a lattice rule is a random shift modulo 1 [3, 38, 17]: we generate one random point $\mathbf{U}$ uniformly in $(0, 1)^s$ and add it to each point $\mathbf{u}_i$ modulo 1, coordinate-wise. This preserves the lattice structure of the point set. The tag “Lat-RS” will refer to this type of RQMC point set. It is often recommended to add a baker (or tent) transformation after the random shift. We replace each $\mathbf{u}_{i,j}$ by tent($\mathbf{u}_{i,j}$) = $1-2|\mathbf{u}_{i,j}-1/2|$. As $u$ ranges from 0 to 1, tent($u$) goes from 0 to 1 and then back to 0. With this transformation, a continuous function over $[0, 1]^s$ becomes continuous and periodic over the $s$-dimensional unit torus, and this greatly improves the convergence rate of the RQMC estimator [10, 15]. We refer to this type of point set by “Lat-RSB.” For a non-periodic integrand $f$ for which all mixed partial derivatives of order 1 with respect to any subset of coordinates are square integrable over $(0, 1)^s$, for instance, the RQMC variance converges as $O(n^{-2+\epsilon})$ with Lat-RS and $O(n^{-4+\epsilon})$ with Lat-RSB. This is a huge gain. It also changes substantially the distribution of the RQMC estimator, which can have a significant impact when computing a confidence interval.
cases are very far from those where nonparametric confidence intervals have been
previously studied.

The remainder is organized as follows. In Section 2, we recall the definitions of the point sets and randomizations considered in this paper, and the variance bounds that are known for the relevant combinations. In Section 3, we discuss several ways of computing confidence intervals RQMC estimators, including via a Student *t* approximation and by applying the bootstrap *t*. In Section 4, we illustrate numerically the behavior of these methods for some functions *f*. Our conclusions are in Section 5.

2 POINT SETS AND RANDOMIZATION METHODS

Here we define the QMC point sets and the randomizations considered in this pa-
per. For each type of RQMC points, we define a short tag (acronym) to use later
when reporting empirical results.

Our first type of QMC point set corresponds to a *lattice rule of rank 1*. It is
defined by selecting the number *n* of points and a generating vector *a* = (*a*₁, ..., *a*₉)
of integers *a*ⱼ ∈ {1, ..., *n* − 1}. Then the *i*th point *u*ᵢ is defined by *u*ᵢ = (*i*α mod *n*)/*n*, for *i* = 0, ..., *n* − 1. In our experiments, the vector *a* was selected using
LatNet Builder [20, 19] with the *P*α criterion with α = 2, product weights, and a
weight of γⱼ = 2/(2 + *j*) for coordinate *j*. See [19] for the details.

The standard randomization for a lattice rule is a random shift modulo 1 [3, 38, 17]: we generate one random point **U** uniformly in (0, 1)⁴ and add it to each point **u**ᵢ modulo 1, coordinate-wise. Adding the same uniform **U** to each of the *n* points preserves the lattice structure of the point set, with the *n* resulting points dependent. Moreover, each randomly shifted point is uniformly distributed on (0, 1)⁴, making the average over all *n* randomly shifted points an unbiased estimator of μ by the linearity of expectation. The tag “Lat-RS” will refer to this type of RQMC point set. It is often recommended to add a baker (or tent) transformation after the random shift. We replace each uᵢ,ⱼ by tent(uᵢⱼ) = 1 - 2|uᵢⱼ - 1/2|. As *u* ranges from 0 to 1, tent(*u*) goes from 0 to 1 and then back to 0. With this transformation, a continuous function over [0, 1]⁴ becomes continuous and periodic over the *d*-dimensional unit torus, and this greatly improves the convergence rate of the RQMC estimator [10, 15]. We refer to this type of point set by “Lat-RSB.” For a non-periodic integer *f* for which all mixed partial derivatives of order 1 with respect to any subset of coordinates are square integrable over (0, 1)⁴, for instance, the RQMC variance converges as O(n⁻²+ε) with Lat-RS, and as O(n⁻⁴+ε) with Lat-RSB. This is a huge gain. It also changes substantially the distribution of the RQMC estimator, which can have a significant impact when computing a confidence interval.
For a digital net in base 2, the coordinate j of point i is

$$u_{i,j} = \sum_{\ell=1}^{w} u_{i,j,\ell} 2^{-\ell}, \quad i = 0, \dots, b^k - 1 \text{ and } j = 1, \dots, s,$$

where

$$\begin{pmatrix} u_{i,j,1} \\ \vdots \\ u_{i,j,w} \end{pmatrix} = C_j \begin{pmatrix} a_{i,0} \\ \vdots \\ a_{i,k-1} \end{pmatrix} \mod 2 \quad \text{with } i = \sum_{\ell=0}^{k-1} a_{i,\ell} 2^{\ell},$$

and the $C_j$ are $w \times k$ generating matrices of rank $k$, with elements in $\{0, 1\}$. See [5, 27, 15] for more details. This gives $n = b^k$ points. In our implementation, $w = 31$. For the $C_j$, we use the popular Sobol' constructions, with the direction numbers proposed in [12].

There are many ways of randomizing a digital net while keeping the digital net structure and respecting the RQMC conditions [31, 15]. The simplest is a random digital shift in base 2: generate a single random point $U = (U_1, \dots, U_j)$ uniformly over $(0, 1)^s$, and perform a bitwise xor of each coordinate of this point with the corresponding coordinates of the points $u_i$ of the digital net. We refer to this one by “Sob-DS.” A second type of approach randomizes the generating matrices $C_j$ before generating the points [31]. The most popular way to do that is to pre-multiply each $C_j$ (mod 2) by a random lower triangular $w \times w$ matrix $L_j$ with 1's on the diagonal and independent random binary entries below the diagonal [24]. We refer to it by “Sob-LMS.” The rationale is that the Sobol' (or other given) generating matrices often have too much structure, and randomizing them improves their quality on average. Note that to obtain an unbiased RQMC estimator with LMS, it is essential to perform a random digital shift after randomizing the generating matrices. A third approach is the nested uniform scramble (NUS) [29] also described in [5] and [15], for example. We refer to it by “Sob-NUS.” For both NUS and LMS, under sufficient smoothness conditions on $f$, the RQMC variance converges as $\mathcal{O}(n^{-3}(\log n)^{s-1})$ [30, 41]. There is also a strong law of large numbers for some randomizations of digital nets [34].

For all these RQMC methods, a standard CLT with a normal limiting distribution holds when $n$ is fixed and $R \to \infty$, because the RQMC estimator is an average of $R$ iid unbiased estimators. But for $n \to \infty$ for fixed $R$, this type of CLT has been proved only for NUS [23], and only for nets with $t = 0$ and for smooth integers. For randomly-shifted lattice rules, when $n \to \infty$ for $R = 1$, it is shown that the RQMC estimator usually has a limiting distribution whose properly-scaled cumulative distribution function (CDF) is a spline function, not a normal distribution [21]. Figure 10 in [15] gives an example of an Lat-RS estimator whose distribution is far from normal. Sufficient conditions are provided in [26] for a CLT to hold in
The next QMC point set we describe is a digital net in base $b=2$. The coordinate $j$ of point $i$ is

$$u_{i,j} = \sum_{\ell=1}^{w} u_{i,j,\ell} 2^{-\ell}, \quad i=0, \dots, 2^k-1 \text{ and } j=1, \dots, d,$$

where

$$\begin{pmatrix} u_{i,j,1} \\ \vdots \\ u_{i,j,w} \end{pmatrix} = C_j \begin{pmatrix} a_{i,0} \\ \vdots \\ a_{i,k-1} \end{pmatrix} \mod 2 \quad \text{with } i = \sum_{\ell=0}^{k-1} a_{i,\ell} 2^{\ell},$$

and the $C_j$ are $w \times k$ generating matrices of rank $k$, with elements in $\{0,1\}$. See [5, 27, 15] for more details. This gives $n = 2^k$ points. Our implementation used $w=31$. For the $C_j$, we use the popular Sobol' constructions, with the direction numbers proposed in [12].

There are many ways of randomizing a digital net while keeping the digital net structure and respecting the RQMC conditions [31, 15]. The simplest is a random digital shift in base $b=2$: generate a single random point $U = (U_1, \dots, U_j)$ uniformly over $(0,1)^d$, and perform a bitwise xor of each coordinate of this point with the corresponding coordinates of the points $u_i$ of the digital net. We refer to this one by “Sob-DS.” It already provides an unbiased RQMC estimator. The variance of this estimator can often be improved by also randomizing the generating matrices $C_j$ before generating the points [31]. The most popular way to do that is to pre-multiply each $C_j$ (mod 2) by a random lower triangular $w \times w$ matrix $L_j$ with 1's on the diagonal and independent random binary entries below the diagonal [24]. We refer to it by “Sob-LMS,” where LMS stands for “left matrix scramble.” The rationale is that the Sobol' (or other given) generating matrices often have too much structure, and randomizing them improves their quality on average. Note that to obtain an unbiased RQMC estimator with LMS, it is essential to perform a random digital shift after randomizing the generating matrices. A third approach is the nested uniform scramble (NUS) [29] also described in [5] and [15], for example. We refer to it by “Sob-NUS.” For both NUS and LMS, under sufficient smoothness conditions on $f$, the RQMC variance converges as $\mathcal{O}(n^{-3}(\log n)^{d-1})$ [30, 41]. There is also a strong law of large numbers for some randomizations of digital nets [34].

For all these RQMC methods, a standard CLT with a normal limiting distribution holds when $n$ is fixed, $\hat{\mu}_{n,\text{rqmc}}$ has finite positive variance, and $R \to \infty$, because the RQMC estimator is an average of $R$ iid unbiased estimators. But for $n \to \infty$ for fixed $R$, this type of CLT has been proved only for NUS [23], and only for special constructions known as $(0,m,d)$-nets and smooth integrands; see, e.g., Section 15.7 of [33] for details. For randomly-shifted lattice rules, when $n \to \infty$ for $R=1$, the
an asymptotic regime where both $n$ and $R$ increase together ($R$ must increase fast enough).

### 3 CONFIDENCE INTERVALS FOR RQMC

As outlined above, the challenge in forming confidence intervals for RQMC comes from a harsh tradeoff. When using $R$ replicates of an RQMC method with $n$ points, at cost $nR$, we expect the greatest accuracy in our estimate for large $n$ and small $R$, while confidence intervals will become more reliable with larger $R$. The sample size for our confidence intervals is $R$. We therefore would like reliable confidence intervals for the mean based on a small sample of $R$ iid variables. It is known [1] that there do not exist nonparametric confidence intervals for the mean. The usable methods have asymptotic justifications, but our setting considers $R$ not large. Our confidence intervals take the form $[L, U]$ for statistics $L = L(y_1, \dots, y_R)$ and $U = U(y_1, \dots, y_R)$, where the nominal confidence level is $1 - \alpha$ for a given $\alpha \in (0, 1)$, e.g., $\alpha = 0.05$. Because these intervals cannot be perfect, we study their coverage error

$$ \varepsilon := \Pr(L \le \mu \le U) - (1 - \alpha). $$

Here we describe confidence intervals for the mean of IID values $y_r$. In our setting, $y_r = \hat{\mu}_r$ for $r = 1, \dots, R$. Because RQMC estimates are unbiased, the common mean $\mathbb{E}(y_r)$ is the integral estimate $\mu$ that we seek. Given $y_1, \dots, y_R$, let

$$ \bar{y} = \frac{1}{R} \sum_{r=1}^{R} y_r \quad \text{and} \quad s^2 = \frac{1}{R-1} \sum_{r=1}^{R} (y_r - \bar{y})^2. $$

The standard approximate $100(1 - \alpha)\%$ confidence interval based on the $t$-test has

$$ L = \bar{y} - t_{(R-1)}^{1-\alpha/2} s R^{-1/2} \quad \text{and} \quad U = \bar{y} + t_{(R-1)}^{1-\alpha/2} s R^{-1/2} $$

where $t_{(R-1)}^{1-\alpha/2}$ is the $1 - \alpha/2$ quantile of Student's $t$ distribution on $R - 1$ degrees of freedom. This interval is constructed to be exact for iid Gaussian data $y_r \stackrel{\text{iid}}{\sim} N(\mu, \sigma^2)$, i.e., normal with mean $\mu$ and variance $\sigma^2 > 0$. It has an asymptotic justification based on the CLT when $y_r$ has finite variance.

The second method we consider is a kind of percentile bootstrap. A bootstrap sample of $y_1, \dots, y_R$ is a list of values $y_1^*, \dots, y_R^*$ sampled independently (with replacement) from $y_1, \dots, y_R$. If there are repeated values among the $y_r$, the sampling respects that multiplicity. Operationally, we sample indices $\ell_1^*, \dots, \ell_R^* \stackrel{\text{iid}}{\sim} \mathcal{U}\{1, 2, \dots, R\}$ (i.e., uniformly from the discrete set $\{1, 2, \dots, r\}$), then take $y_r^* =$
RQMC estimator usually has a limiting distribution whose properly-scaled cumulative distribution function (CDF) is a spline function, not a normal distribution [21]. [15], Figure 10, gives an example of an Lat-RS estimator whose distribution is far from normal. [26] provide sufficient conditions for a CLT to hold in an asymptotic regime where both n and R increase together (R must increase fast enough).

### 3 CONFIDENCE INTERVALS FOR RQMC

As outlined above, the challenge in forming confidence intervals for RQMC comes from a harsh tradeoff. When using R replicates of an RQMC method with n points, at cost nR that is constrained to not exceed a given computation budget, we expect the greatest accuracy in our estimate for large n and small R, while confidence intervals will become more reliable with larger R. The sample size for our confidence intervals is R. We therefore would like reliable confidence intervals for the mean based on a small sample of R iid variables. It is known since [1] that there do not exist nonparametric confidence intervals with exact coverage for the mean. The usable methods have asymptotic justifications, but our setting considers R not large. Our confidence intervals take the form [L, U] for statistics L = L(y₁, ..., yR) and U = U(y₁, ..., yR), where the nominal confidence level is 1 − α for a given α ∈ (0, 1), e.g., α = 0.05. Because these intervals cannot be perfect, we study their coverage error

$$ \varepsilon := \Pr(L \le \mu \le U) - (1 - \alpha). $$

Here we describe confidence intervals for the mean of IID values yᵣ. In our setting, yᵣ = ūr for r = 1, ..., R. Because RQMC estimates are unbiased, the common mean E(yᵣ) is the integral value μ that we seek. Given y₁, ..., yR, let

$$ \bar{y} = \frac{1}{R} \sum_{r=1}^{R} y_{r} \quad \text{and} \quad S^{2} = \frac{1}{R-1} \sum_{r=1}^{R} (y_{r} - \bar{y})^{2}. $$

The standard approximate 100(1 − α)% confidence interval for μ based on the t-test has

$$ L = \bar{y} - t_{(R-1)}^{1-\alpha/2} SR^{-1/2} \quad \text{and} \quad U = \bar{y} + t_{(R-1)}^{1-\alpha/2} SR^{-1/2} $$

where t(1 − α/2) is the 1 − α/2 quantile of the Student t distribution on R − 1 degrees of freedom. This interval is constructed to be exact for iid Gaussian data yᵣ iid ~ N(μ, σ²), i.e., normal with mean μ and variance σ² > 0. It has an asymptotic justification based on the CLT when yᵣ has finite variance.
$y_{r}^{*}$, and compute $\bar{y}^{*} = (1/R) \sum_{r=1}^{R} y_{r}^{*}$. The idea is that in this process the sample values $y_{r}^{*}$ are to the empirical distribution of $y_{r}$ what those sample values are to the true distribution of $y_{r}$. In the percentile method we repeat this resampling $B$ times independently, getting $\bar{y}^{*b}$ for $b = 1, \dots, B$. Sorting them yields $\bar{y}^{*(1)} \le \bar{y}^{*(2)} \le \dots \le \bar{y}^{*(B)}$. The confidence interval endpoints are then the quantiles $\bar{y}^{*([B\alpha/2])}$ and $\bar{y}^{*([B(1-\alpha)/2])}$, where $\lfloor \cdot \rfloor$ and $\lceil \cdot \rceil$ are the floor and ceiling functions. Confidence intervals based on small sample sizes $R$ commonly use $B = 1000$. So, for $B = 1000$ and $\alpha = 0.05$, we would take $L = \bar{y}^{*(25)}$ and $U = \bar{y}^{*(975)}$. Some authors prefer to take $B = 999$ in this case and then use $L = \bar{y}^{*(25)}$ and $U = \bar{y}^{*(975)}$. The rationale is that the $B$ points $\bar{y}^{*b}$ partition $\mathbb{R}$ into $B+1$ intervals and the confidence interval is then the union of the central 95% of those intervals.

The third method we consider is the bootstrap $t$ method of Efron [7]. It is recommended for RQMC without much analysis in [33]. It outperforms the standard method based on the Student distribution when estimating a ratio of expectations in a regenerative simulation setting [2]. It is known to be highly accurate for bootstrapping the mean, both asymptotically [8] and for small sample sizes, where for all but one distribution, [28] finds empirically that it works reasonably well for sample sizes as low as 5. It also has a stronger rationale than the percentile method. The reasoning behind the bootstrap $t$ is that the sample distribution of the $t$ statistic $\sqrt{R}(\bar{y}-\mu)/s$ is well approximated by the sample distribution of a bootstrapped $t$ statistic $\sqrt{R}(\bar{y}^*-\bar{y})/s^*$, where $s^*$ is the standard deviation of $y_1^*, \dots, y_R^*$. We take $B$ independent bootstrap $t$ values $t_L^{*b}$ for $b = 1, \dots, B$, sort them, and then let $t_L^*$ and $t_U^*$ be the $\alpha/2$ and $1-\alpha/2$ quantiles of the $t^{*b}$ values. By taking $B$ large enough, we can make

$$ \Pr(t_L^* \le \sqrt{R}\frac{\bar{y}^* - \bar{y}}{s^*} \le t_U^*) $$

as close as we like to $1-\alpha$. Then if we reason that $t_L^* \le \sqrt{R}(\bar{y}-\mu)/s \le t_U^*$ has approximately $1-\alpha$ probability of holding, we take

$$ L = \bar{y} - st_U^* R^{-1/2} \quad \text{and} \quad U = \bar{y} - st_L^* R^{-1/2} $$

to form our approximate confidence interval for $\mu$. Note that the upper limit for $t^*$ appears in the lower limit for $\mu$ and vice versa.

Conditions under which the coverage error $\epsilon$ for these methods is $O(1/R)$ are given in [8]. It is remarkable that the confidence interval is more accurate than $\bar{y}$ itself, which has error $O_p(1/\sqrt{R})$ as $R \to \infty$. It is difficult to connect the tabular values in [8] to the names that others use for those methods. The bootstrap $t$ method is called ‘STUD’ for studentized; it is also called ‘percentile $t$’ because it used percentiles of a bootstrapped $t$ statistic. The percentile method is called “the other percentile method” in [9], who then notes that it is called ‘BACK’ in [8], who
The second method we consider is a kind of *percentile bootstrap*. A bootstrap sample of $y_1, \dots, y_R$ is a list of values $y_1^*, \dots, y_R^*$ sampled independently (with replacement) from $y_1, \dots, y_R$. If there are repeated values among the $y_r$, the sampling respects that multiplicity. Operationally, we sample indices $\ell_1^*, \dots, \ell_R^*$ iid $\mathcal{U}\{1, 2, \dots, R\}$ (i.e., uniformly from the discrete set $\{1, 2, \dots, r\}$), then take $y_r^* = y_{\ell_r^*}$, and compute $\bar{y}^* = (1/R) \sum_{r=1}^R y_r^*$. The idea is that in this process the sample values $y_r^*$ are to the empirical distribution of $y_r$ what those sample values are to the true distribution of $y_r$. In the percentile method we repeat this resampling $B$ times independently, getting $\bar{y}^{*b}$ for $b = 1, \dots, B$. Sorting them yields $\bar{y}^{*(1)} \le \bar{y}^{*(2)} \le \dots \le \bar{y}^{*(B)}$. The confidence interval endpoints are then the quantiles $\bar{y}^{*([B\alpha/2])}$ and $\bar{y}^{*([B(1-\alpha)/2])}$, where $\lfloor \cdot \rfloor$ and $\lceil \cdot \rceil$ are the floor and ceiling functions. Confidence intervals based on small sample sizes $R$ commonly use $B = 1000$. So, for $B = 1000$ and $\alpha = 0.05$, we would take $L = \bar{y}^{*(25)}$ and $U = \bar{y}^{*(975)}$. Some authors prefer to take $B = 999$ in this case and then use $L = \bar{y}^{*(25)}$ and $U = \bar{y}^{*(975)}$. The rationale is that the $B$ points $\bar{y}^{*b}$ partition $\mathbb{R}$ into $B+1$ intervals and the confidence interval is then the union of the central 95% of those intervals.

The third method we consider is the *bootstrap t method* [7]. It is recommended for RQMC without much analysis [33]. It is shown in [2] that it outperforms the standard method based on the Student $t$ distribution when estimating a ratio of expectations in a regenerative simulation setting. It is known to be highly accurate for bootstrapping the mean, both asymptotically [8] and for small sample sizes, where for all considered distributions but one, In [28], it is empirically found that it works reasonably well for sample sizes as low as 5. It also has a stronger rationale than the percentile method. The reasoning behind the bootstrap $t$ is that the sample distribution of the $t$ statistic $\sqrt{R}(\bar{y} - \mu)/S$ is well approximated by the sample distribution of a bootstrapped $t$ statistic $\sqrt{R}(\bar{y}^* - \bar{y})/S^*$, where $S^*$ is the standard deviation of $y_1^*, \dots, y_R^*$. We take $B$ independent bootstrap $t$ values $t^{*b}$ for $b = 1, \dots, B$, sort them, and then let $t_L^*$ and $t_U^*$ be the $\alpha/2$ and $1 - \alpha/2$ quantiles of the $t^{*b}$ values. By taking $B$ large enough, we can make

$$ \Pr(t_L^* \le \sqrt{R} \frac{\bar{y}^* - \bar{y}}{S^*} \le t_U^*) $$

as close as we like to $1 - \alpha$. Then if we reason that $t_L^* \le \sqrt{R}(\bar{y} - \mu)/S \le t_U^*$ has approximately $1 - \alpha$ probability of holding, we take

$$ L = \bar{y} - St_U^* R^{-1/2} \quad \text{and} \quad U = \bar{y} - St_L^* R^{-1/2} $$

to form our approximate confidence interval for $\mu$. Note that the upper limit for $t^*$ appears in the lower limit for $\mu$ and vice versa.

Conditions under which the coverage error $\epsilon$ for these methods is $\mathcal{O}(1/R)$ are given in [8]. It is remarkable that the confidence interval is more accurate than $\bar{y}$
interprets the non-reversal of quantiles to mean that it is looking at the quantile ap-
proximations backwards. Reference [8] does not include the usual t-test confidence
intervals but does include their normal theory counterparts which use quantiles of
the N(0, 1) distribution in place of the t distribution. With these assumptions, the
coverage errors from Table 1 in [8] are

$$
\text{Normal theory:} \quad (1/R)\varphi(z^{1-\alpha/2}) \left[ 0.14\kappa - 2.12\gamma^2 - 3.35 \right] + \mathcal{O}(1/R^2),
$$

$$
\text{Percentile: } (1/R)\varphi(z^{1-\alpha/2})[-0.72\kappa - 0.37\gamma^2 - 3.35] + \mathcal{O}(1/R^2), \text{ and}
$$

$$
\text{Bootstrap } t: \quad (1/R)\varphi(z^{1-\alpha/2})[-2.84\kappa + 4.25\gamma^2] + \mathcal{O}(1/R^2),
$$

where $\gamma = \mathbb{E}((\hat{\mu}_r - \mu)^3)/\sigma^3$ and $\kappa = \mathbb{E}((\hat{\mu}_r - \mu)^4)/\sigma^4 - 3$ are the skewness and kurtosis, respectively, of $\hat{\mu}_r$. These rates hold for central confidence intervals that are our focus here. For one-sided intervals, the percentile method has errors that are $\mathcal{O}(R^{-1/2})$, while the bootstrap $t$ method has one-sided errors that are $\mathcal{O}(R^{-1})$. Hall's conditions involve 8 finite moments for $\hat{\mu}_r$ and an assumption that the support of $\hat{\mu}_r$ is not an arithmetic sequence. Here $\varphi$ is the $N(0,1)$ density function and $z^{1-\alpha/2} = \Phi^{-1}(1-\alpha/2)$ is the $(1-\alpha/2)$-quantile of the $N(0,1)$ CDF $\Phi$. Hall's $\alpha$ is twice ours.

The normal theory method above has error $-3.35\varphi(z^{1-\alpha/2})/R + \mathcal{O}(1/R^2)$ for Gaussian data (because it should have used the $t$ threshold). The bootstrap $t$ has an advantage in missing the $-3.35$ component that the others have. It has error $\mathcal{O}(R^{-2})$ on Gaussian data. It has a large positive coefficient for $\gamma^2$ where the others have negative coefficients. This gives it extra coverage on highly skewed data. It is well known to produce very long intervals; in addition to the asymptotics, there is an intuitive explanation in that the denominator $s^*$ can become small. It can even be zero.

The asymptotics predict that the bootstrap *t* will undercover when *κ* is large and
*γ* = 0. In our context, *γ* and *κ* are the skewness and kurtosis of RQMC estimates
based on *n* = 2*<sup>*k*</sup>* function evaluations. In plain Monte Carlo sampling, the average
of *n* random iid variables with skewness *γ* has skewness *γ*/√*n* and the average of
*n* iid random variables with kurtosis *κ* has kurtosis *κ*/*n*. Much less is known about
the skewness and kurtosis of RQMC estimates. Indeed, the work presented here is
the most thorough exploration of that issue yet.

For continuously distributed $y_r$, one can show that $\Pr(s^* = 0) = R^{1-R}$. This is
not negligible for $R = 5$ which we consider. Such $s^*$ ordinarily provide $|t^*| = \infty$.
As long as there are fewer than $B\alpha/2$ bootstrap samples with $t^* = \infty$, the lower
confidence interval will be finite, and similarly, fewer than $B\alpha/2$ bootstrap samples
with $t^* = -\infty$ means that the upper confidence interval will be finite. We revisit this
issue in our numerical work in Section 4 and again in our conclusions in Section 5.
itself, which has error $\mathcal{O}_p(1/\sqrt{R})$ as $R \to \infty$. It is difficult to connect the tabular values in [8] to the names that others use for those methods. The bootstrap t method is called ‘STUD’ for studentized; it is also called ‘percentile t’ because it used percentiles of a bootstrapped t statistic. The percentile method is called “the other percentile method” in [9], who then notes that it is called ‘BACK’ in [8], who interprets the non-reversal of quantiles to mean that it is looking at the quantile approximations backwards. The reference [8] does not include the usual t-test confidence intervals but does include their normal theory counterparts which use quantiles of the $N(0, 1)$ distribution in place of the Student t distribution. With these assumptions, the coverage errors from Table 1 of [8] are

Normal theory: $(1/R)\phi(z^{1-\alpha/2})[0.14\kappa - 2.12\gamma^2 - 3.35] + \mathcal{O}(1/R^2)$,

Percentile: $(1/R)\phi(z^{1-\alpha/2})[-0.72\kappa - 0.37\gamma^2 - 3.35] + \mathcal{O}(1/R^2)$, and

Bootstrap t: $(1/R)\phi(z^{1-\alpha/2})[-2.84\kappa + 4.25\gamma^2] + \mathcal{O}(1/R^2)$,

where $\gamma = \mathbb{E}((\hat{\mu}_r - \mu)^3)/\sigma^3$ and $\kappa = \mathbb{E}((\hat{\mu}_r - \mu)^4)/\sigma^4 - 3$ are the skewness and kurtosis, respectively, of $\hat{\mu}_r$, $\phi$ is the $N(0, 1)$ density function, and $z^{1-\alpha/2} = \Phi^{-1}(1-\alpha/2)$ is the $(1-\alpha/2)$-quantile of the $N(0, 1)$ CDF $\Phi$. These rates hold for central confidence intervals that are our focus here. For one-sided intervals, the percentile method has errors that are $\mathcal{O}(R^{-1/2})$, while the bootstrap t method has one-sided errors that are $\mathcal{O}(R^{-1})$. Hall’s conditions involve 8 finite moments for $\hat{\mu}_r$ and an assumption that the support of $\hat{\mu}_r$ is not an arithmetic sequence. Hall’s $\alpha$ is twice ours.

The normal theory method above has error $-3.35\phi(z^{1-\alpha/2})/R + \mathcal{O}(1/R^2)$ for Gaussian data (because it should have used instead the Student t threshold). The bootstrap t has an advantage in missing the $-3.35$ component that the others have. It has error $\mathcal{O}(R^{-2})$ on Gaussian data since $\gamma = \kappa = 0$. It has a large positive coefficient for $\gamma^2$ where the others have negative coefficients. This gives it extra coverage on highly skewed data. The bootstrap t is well known to produce very long intervals; in addition to the asymptotics, there is an intuitive explanation in that the denominator $S^*$ can become small. It can even be zero.

The asymptotics predict that the bootstrap t will undercover when $\kappa$ is large and $\gamma = 0$. In our context, $\gamma$ and $\kappa$ are the skewness and kurtosis of RQMC estimates based on $n = 2^k$ function evaluations. In plain Monte Carlo sampling, the average of $n$ random iid variables with skewness $\gamma$ has skewness $\gamma/\sqrt{n}$ and the average of $n$ iid random variables with kurtosis $\kappa$ has kurtosis $\kappa/n$. Much less is known about the skewness and kurtosis of RQMC estimates. Indeed, the work presented here is the most thorough exploration of that issue yet.

For continuously distributed $y_r$, one can show that $\Pr(S^* = 0) = R^{1-R}$. This is not negligible for $R = 5$ which we consider. Such $S^*$ ordinarily provide $|t^*| = \infty$.
# 4 NUMERICAL EXPERIMENTS

For our numerical experiments, we made a strategic selection of test functions $f$ described in Section 4.1. They can all be parameterized by the dimension $s$, and the integral $\mu$ in equation (1) can be computed exactly for any $s$. Subtracting $\mu$ from $f$ gives a function with integral 0 in all dimensions that we can use to estimate the coverage of our confidence intervals. For each selected (centered) function $f$, we tried the five types of RQMC point sets described in Section 2, with $n = 2^k$ points for $k = 6, 8, 10, 12, 14$, and in $s = 4, 8, 16, 32$ dimensions.

Our computations had two steps. In the first step, for each combination of $(f, \text{method}, k, s)$, we generated a pool of $N = 10^4$ independent realizations $y_1, \dots, y_N$ of the RQMC estimator $y = \hat{\mu}_{n,\text{rqmc}}$ in (2) and stored the sorted values in files. To do this, we used the packages `hups` and `mcqmctools` from the SSJ Java library [16]. These large samples were used in the second step to visualize the distribution of $y$ and estimate its variance, skewness, and kurtosis, in each case.

In the second step, we generated small samples of $R$ replicates (without replacement) from the pool of $N$ values and used them to form all three kinds of confidence intervals from Section 3. The sample sizes we used were $R \in \{5, 10, 20, 30\}$. We did that 1000 times. This strategy lets us generate more than $N/R$ confidence intervals using our pool of RQMC estimates. For each such sample, we computed a 95% confidence interval with each of the three procedures. From these 1000 replicated confidence intervals we could record the proportion of the intervals that covered 0. We also computed the lengths of all the confidence intervals. Each of our 1000 bootstrap confidence intervals was based on $B = 1000$ bootstrap resamples. The whole process of forming 1000 confidence intervals of each type was done separately for all the integrands, methods, dimensions $s$, and RQMC sample sizes $n = 2^k$ described above.

## 4.1 Selected Functions

Here are the functions we consider. In each case, the integrand depends on a point $\mathbf{u} \in (0, 1)^s$ and the integral we seek is $\int_{(0,1)^s} f(\mathbf{u})d\mathbf{u}$. For each function, we subtract an appropriate constant to make the integral equal to zero in all cases. This constant is already included in the definitions below. Some of the test functions are constructed using $\Phi$ and $\Phi^{-1}$, the CDF and inverse CDF of the standard normal distribution. Some of the functions are easy for RQMC and some are more difficult. A function that is easy for RQMC is not necessarily one that is easy to get a confidence interval for.

1. SumUeU: $f(\mathbf{u}) = -s + \sum_{j=1}^{s} u_j \exp(u_j)$. This is smooth and additive, making it very easy for RQMC methods to integrate; e.g., see Section 17.2 of [33].
As long as there are fewer than $B\alpha/2$ bootstrap samples with $t^* = \infty$, the lower confidence interval will be finite, and similarly, fewer than $B\alpha/2$ bootstrap samples with $t^* = -\infty$ means that the upper confidence interval will be finite. We revisit this issue in our numerical work in Section 4 and again in our conclusions in Section 5.

# 4 NUMERICAL EXPERIMENTS

For our numerical experiments, we made a strategic selection of test functions $f$ described in Section 4.1. They can all be parameterized by the dimension $d$, and the integral $\mu$ in equation (1) can be computed exactly for any $d$. Subtracting $\mu$ from $f$ gives a function with integral 0 in all dimensions, which we can use to estimate the coverage of our confidence intervals. For each selected (centered) function $f$, we tried the five types of RQMC point sets described in Section 2, with $n = 2^k$ points for $k = 6, 8, 10, 12, 14$, and in $d = 4, 8, 16, 32$ dimensions.

Our computations had two steps. In the first step, for each combination of $(f, \mu)$, we generated a pool of $N = 10^4$ independent realizations $y_1, \dots, y_N$ of the RQMC estimator $\hat{y} = \hat{\mu}_{n,\text{rqmc}}$ in (2) and stored the sorted values in files. To do this, we used the packages `hups` and `mcqmctools` from the SSJ Java library [16]. These large samples were used in the second step to visualize the distribution of $\hat{y}$ and estimate its variance, skewness, and kurtosis, in each case.

In the second step, we generated small samples of $R$ replicates (without replacement) from the pool of $N$ values and used them to form all three kinds of confidence intervals from Section 3. The sample sizes we used were $R \in \{5, 10, 20, 30\}$. We did that 1000 times. This strategy lets us generate more than $N/R$ confidence intervals using our pool of RQMC estimates. For each such sample, we computed a 95% confidence interval with each of the three procedures. From these 1000 replicated confidence intervals we could record the proportion of the intervals that covered 0. We also computed the lengths of all the confidence intervals. Each of our 1000 bootstrap confidence intervals was based on $B = 1000$ bootstrap samples. The whole process of forming 1000 confidence intervals of each type was done separately for all the integrands, methods, dimensions $d$, and RQMC sample sizes $n = 2^k$ described above.

## 4.1 Selected Functions

Here are the functions we consider. In each case, the integrand depends on a point $\mathbf{u} \in (0,1)^d$ and the integral we seek is $\int_{(0,1)^d} f(\mathbf{u})d\mathbf{u}$. For each function, we subtract an appropriate constant to make the integral equal to zero in all cases. This constant is already included in the definitions below. Some of the test functions are
2. MC2: $f(\boldsymbol{u}) = -1 + (s - 1/2)^s \prod_{j=1}^{s} (x_j - 1/2)$. This function is from [25]. It is smooth and very nearly, but not exactly, additive. This makes it typical of the integrands in weighted function spaces, such as the spaces of [37].

3. PieceLinGauss: $f(\boldsymbol{u}) = \max(s^{-1/2}\sum_{j=1}^{s} \Phi^{-1}(u_j) - \tau, 0) - \varphi(\tau) + \tau\Phi(-\tau)$. This function is piecewise linear and continuous in $s$ Gaussian inputs. It has two sources of infinite variation in the sense of Hardy and Krause [27, Section 2.2]. It is unbounded and it has a cusp that brings infinite variation for $s \ge 3$. It shares these properties with some integrands encountered in computational finance, which are analytically intractable. Our construction gives this integral a known mean for all $s$. We use $\tau = 1$. Larger values for $\tau$ would give the problem a “rare events” character that is outside the scope of this article.

4. IndSumNormal: $f(\boldsymbol{u}) = -\Phi(-\tau) + \mathbb{I}\{s^{-1/2} \sum_{j=1}^{s} \Phi^{-1}(u_j) \ge \tau\}$, where $\mathbb{I}$ denotes the indicator function. This function is discontinuous with infinite variation in the sense of Hardy and Krause for $s \ge 2$. The discontinuities are not axis parallel so that they are not QMC-friendly. We use $\tau = 1$. The discreteness in the integrand produces a discrete distribution for the integral estimates. That lets us test a violation of one of Hall’s conditions.

5. SmoothGauss: $f(\boldsymbol{u}) = -\Phi(1/\sqrt{2}) + \Phi(1 + s^{-1/2} \sum_{j=1}^{s} \Phi^{-1}(u_j))$. This function is smooth and bounded and monotone in each $u_j$. It has bounded variation.

6. RidgeJohnsonSU: $f(\boldsymbol{u}) = -\eta + F^{-1}(s^{-1/2} \sum_{j=1}^{s} u_j)$ where $F$ is the CDF of the Johnson’s SU distribution [13] with parameters $\gamma = \delta = \lambda = 1$ and $\xi = 0$, and $\eta$ is the mean of that distribution. It has skewness $-5.66$ and kurtosis $96.8$ (for any $s$) making it heavy tailed.

The last four examples use “ridge functions” of the form $g(\theta^T \Phi^{-1}(\boldsymbol{u}))$ for a unit vector $\theta \in \mathbb{R}^s$ and a function $g: \mathbb{R} \to \mathbb{R}$, where a superscript T denotes transpose. In a plain MC setting they, before centering, have the distribution $g(N(0,1))$. If $g$ is Lipschitz, then the ridge function has bounded mean dimension (a favorable property for RQMC) as the nominal dimension $s \to \infty$ [11], while if $g$ is discontinuous the mean dimension can grow proportionally to $\sqrt{s}$. In the ridge functions above, every variable is equally important, which tends to make integration harder. To introduce some sparsity, let $\eta_j = 2^{-j}$ and $\theta_j = \eta_j / \sqrt{\sum_{j'=1}^{s} \eta_{j'}^2}$. Then we can have sparse versions of the last three test functions. Sparsity does not affect the value of the integral. We do not include sparse versions in this work.
constructed using $\Phi$ and $\Phi^{-1}$, the CDF and inverse CDF of the standard normal
distribution. Some of the functions are easy for RQMC, while others are more dif-
ficult. A function that is easy for RQMC is not necessarily one that is easy to get a
confidence interval for.

1. SumUeU: $f(\mathbf{u}) = -d + \sum_{j=1}^{d} u_j \exp(u_j)$. This is smooth and additive, making it very easy for RQMC methods to integrate; e.g., see Section 17.2 of [33].

2. MC2: $f(\mathbf{u}) = -1 + (d-1/2)^{-d} \prod_{j=1}^{d} (d-x_j)$. This function is from [25]. It is smooth and very nearly, but not exactly, additive. This makes it typical of the integrands in weighted function spaces, such as the spaces of [37].

3. PieceLinGauss: $f(\mathbf{u}) = \max(d^{-1/2} \sum_{j=1}^{d} \Phi^{-1}(u_j) - \tau, 0) - \varphi(\tau) + \tau\Phi(-\tau)$. This function is piecewise linear and continuous in $d$ Gaussian inputs. It has two sources of infinite variation in the sense of Hardy and Krause [27, Section 2.2]. It is unbounded and it has a cusp that brings infinite variation for $d \ge 3$. It shares these properties with some integrands encountered in computational finance, which are analytically intractable. Our construction gives this integral a known mean for all $d$. We use $\tau = 1$. Larger values for $\tau$ would give the problem a “rare events” character that is outside the scope of this article.

4. IndSumNormal: $f(\mathbf{u}) = -\Phi(-\tau) + \mathbb{I}\{d^{-1/2} \sum_{j=1}^{d} \Phi^{-1}(u_j) \ge \tau\}$, where $\mathbb{I}$ denotes the indicator function. This function is discontinuous with infinite variation in the sense of Hardy and Krause for $d \ge 2$. The discontinuities are not axis parallel so that they are not QMC-friendly. We use $\tau = 1$. The discreteness in the integrand produces a discrete distribution for the integral estimates. That lets us test a violation of one of the conditions of [8].

5. SmoothGauss: $f(\mathbf{u}) = -\Phi(1/\sqrt{2}) + d^{-1/2} \sum_{j=1}^{d} \Phi(1 + \Phi^{-1}(u_j))$. This function is smooth and bounded and monotone in each $u_j$. It has bounded variation.

6. RidgeJohnsonSU: $f(\mathbf{u}) = -\eta + F^{-1}(\Phi(d^{-1/2} \sum_{j=1}^{d} \Phi^{-1}(u_j)))$ where $F$ is the CDF of the Johnson's SU distribution [13] with parameters $\gamma = \delta = \lambda = 1$ and $\xi = 0$, and $\eta$ is the mean of that distribution. It has skewness $-5.66$ and kurtosis $96.8$ (for any $d$) making it heavy tailed.

The last four examples use “ridge functions” of the form $g(\theta^T\Phi^{-1}(\mathbf{u}))$ for a unit vector $\theta \in \mathbb{R}^d$ and a function $g: \mathbb{R} \to \mathbb{R}$, where a superscript T denotes transpose. In a plain MC setting they, before centering, have the distribution $g(N(0,1))$. If
## 4.2 Comparison of Student and Bootstrap Confidence Intervals

Our main interest is in coverage accuracy of the confidence intervals. No nonparametric method can get exactly 95% coverage in finite samples. But if we are reasonably confident that the true coverage is below 94%, then we deem the method to have failed. We judged any method that attained less than 92.7% coverage to have failed. From 1000 independent draws, the binomial distribution would give less than a 4% chance of seeing this failure if the true coverage were 94% (or higher). Our 1000 draws from the pool involve some overlap so this computation is not exact, but we think it is accurate enough.

Each CI method was tested in 2400 tasks: 6 integrands, 5 RQMC methods, 4 dimensions, 5 RQMC sample sizes and 4 values of the replication size *R*. The percentile method failed 1698 (70.75%) of those tasks. This is not very surprising. Its undercoverage was seen in [28], and it is not well suited to very small sample sizes like we would want for RQMC. Despite it being a pretty standard method for some bootstrap computations, it is not well regarded for setting confidence intervals for the mean. We do not consider it further.

The bootstrap *t* method failed 81 times, and the plain Student *t* confidence interval method failed only 3 times. The lowest 15 coverage levels for the plain Student *t* interval were all for *R* = 5 replicates. The bootstrap *t* had coverage over 95% in all those cases, but it also uses very long confidence intervals in those sets.

Of the 81 failures for bootstrap *t*, 74 were for Sob-LMS on either SumUeU (44 times) or MC2 (30 times). These are for the extremely spiky histograms mentioned earlier and included in the illustrations in Section 4.3. There were 6 failures for Sob-NUS and 1 for Sob-DS.

As mentioned above, the bootstrap *t* has the possibility of giving an interval of infinite length. When given a sample of *R* distinct values, it has probability *R*¹⁻*R* of yielding *s** = 0. For *R* = 5, this becomes 0.0016. To get an infinite confidence interval we would need that to happen at least 25 times out of 1000. A binomial calculation shows that this event has probability approximately 1.35 × 10⁻²¹. The probability of *s** = 0 is much higher if two or more of the replicate values it gets are equal to each other. There were 21 bootstrap *t* confidence intervals with infinite length. They were all for the integrand IndSumNormal with *R* = 5. The sample value for this integrand counts the number of events that happen. It has a discrete distribution that may have many fewer than 2*^k^* different values. Getting two equal values among the *R* samples makes it much more likely that we see *s** = 0 and an infinite length interval. The bootstrap *t* had many very long interval lengths for the IndSumNormal integrand. This problem could then arise in bootstrap *t* confidence intervals for the value of a digital option in finance.

A second kind of failure is having coverage too high. This means that the con-
g is Lipschitz, then the ridge function has bounded mean dimension (a favorable property for RQMC) as the nominal dimension $d \to \infty$ [11], while if g is discontinuous the mean dimension can grow proportionally to $\sqrt{d}$. In the ridge functions above, every variable is equally important, which tends to make integration harder. To introduce some sparsity, let $\eta_j = 2^{-j}$ and $\theta_j = \eta_j / \sqrt{\sum_{j'=1}^d \eta_{j'}^2}$. Then we can have sparse versions of the last three test functions. Sparsity does not affect the value of the integral. We do not include sparse versions in this work.

## 4.2 Comparison of Student and Bootstrap Confidence Intervals

Our main interest is in coverage accuracy of the confidence intervals. No non-parametric method can get exactly 95% coverage in finite samples. But if we are reasonably confident that the true coverage is below 94%, then we deem the method to have failed. In our experiments, we judged any method that attained less than 92.7% coverage to have failed. From 1000 independent draws, the binomial distribution would give less than a 4% chance of seeing this failure if the true coverage were 94% (or higher). Our 1000 draws from the pool involve some overlap so this computation is not exact, but we think it is accurate enough.

Each confidence interval method was tested in 2400 tasks: 6 integrands, 5 RQMC methods, 4 dimensions, 5 RQMC sample sizes $n$, and 4 values of the replication size $R$. The percentile method failed 1698 (70.75%) of those tasks. This is not very surprising. Its undercoverage was seen in [28], and it is not well suited to very small sample sizes like we would want for RQMC. Despite it being a pretty standard method for some bootstrap computations, it is not well regarded for setting confidence intervals for the mean. We do not consider it further.

The bootstrap $t$ method failed 81 times, and the plain Student $t$ confidence interval method failed only 3 times. The lowest 15 coverage levels for the plain Student $t$ interval were all for $R = 5$ replicates. The bootstrap $t$ had coverage over 95% in all those cases, but it also uses very long confidence intervals in those sets.

Of the 81 failures for bootstrap $t$, 74 were for Sob-LMS on either SumUeU (44 times) or MC2 (30 times). These are for extremely spiky histograms, which are included in the illustrations in Section 4.3. There were 6 failures for Sob-NUS and 1 for Sob-DS.

As noted in the last paragraph of Section 3, the bootstrap $t$ has the possibility of giving an interval of infinite length. When given a sample of $R$ distinct values, it has probability $R^{1-R}$ of yielding $S^* = 0$. For $R = 5$, this becomes 0.0016. To get an infinite confidence interval we would need that to happen at least 25 times out of $B = 1000$. A binomial calculation shows that this event has probability approximately $1.35 \times 10^{-21}$. The probability of $S^* = 0$ is much higher if two or more of
fidence intervals could have been shorter. It is usually considered a lesser failure
than having coverage too low. The standard intervals had coverage higher than 97%
81 times. This is the threshold at which we are confident that the actual coverage
is above 96%. That only happened for the smooth additive integrand SumUeU and
the smooth nearly additive integrand MC2. Of those 81 cases, 80 of them were for
Sob-LMS and one was for Sob-NUS. For Sob-MLS we know that the underlying
kurtosis of the RQMC points diverges to infinity as *n* increases from results in [35].

4.3 Some Histograms

We looked at many histograms made up of $N = 10^4$ integral estimates from the pools of values that we created. The majority of them looked nearly Gaussian. Some of that is to be expected when one averages $n = 2^k$ different random values. If those had been iid values, then the CLT would apply and the averages would have skewness $\gamma/\sqrt{n}$ and kurtosis $\kappa/n$, where $\gamma$ and $\kappa$ are the skewness and kurtosis of a single random integrand value. The *n* values in RQMC sampling are identically distributed but far from independent.

Figure 1 shows a selection of histograms, primarily the unusual ones. The upper left entry is a somewhat negatively skewed example from the RidgeJohnsonSU integrand. It is shown for the Sob-DS method, but that skewed pattern appeared for other methods too and not just for small *s* and *k*. The upper right entry shows the spike plus outliers distribution mentioned earlier. It also arises for nearly additive MC2 at somewhat higher values of *k* for each *s*. The phenomenon was pointed out in [35] in one dimension and in [36] for multiple dimensions which includes some numerical illustrations. We expect coverage length to be quite variable in settings like this because, for example, *R* = 10 replicates might not include any of the outliers. The middle left entry shows a bimodal distribution that frequently appears for lattice samples with a baker transformation. It was also illustrated in [15] and in Chapter 17.4 of [33]. This bimodality did not always happen in our simulations. In more extreme settings, the distribution is nearly *U* shaped. The middle right entry shows a phenomenon for the IndSumNormal integrand. It is roughly Gaussian plus a spike near one value. For smaller values of *k*, this integrand gives a discrete distribution with a roughly binomial appearance. The lower left entry is a roughly Gaussian histogram, similar to most of those in the data set. The lower right entry is a histogram for Sob-NUS that looks a lot like many that arise for Sob-LMS but which we have never before seen for Sob-NUS.

Some of our examples showed very extreme kurtosis. None of them showed
extreme skewness. The standard confidence interval is known to have robust cover-
age in response to kurtosis at the expense of having very long confidence intervals,
while being vulnerable to skewness. Figure 2 shows a summary of interval length
the replicate values it gets are equal to each other. There were 21 bootstrap *t* confidence intervals with infinite length. They were all for the integrand IndSumNormal with *R* = 5. The sample value for this integrand counts the number of events that happen. It has a discrete distribution that may have many fewer than 2^k^ different values. Getting two equal values among the *R* samples makes it much more likely that we see *S*\* = 0 and an infinite length interval. The bootstrap *t* had many very long interval lengths for the IndSumNormal integrand. This problem could then arise in bootstrap *t* confidence intervals for the value of a digital option in finance.

A second kind of failure is having coverage too high. This means that the confidence intervals could have been shorter. It is usually considered a lesser failure than having coverage too low. The standard intervals had coverage higher than 97% 81 times. This is the threshold at which we are confident that the actual coverage is above 96%. That only happened for the smooth additive integrand SumUeU and the smooth nearly additive integrand MC2. Of those 81 cases, 80 of them were for Sob-LMS and one was for Sob-NUS. For Sob-LMS we know that the underlying kurtosis of the RQMC points diverges to infinity as *n* increases from results in [35].

## 4.3 Some Histograms

We looked at many histograms made up of *N* = 10⁴ integral estimates from the pools of values that we created. The majority of them looked nearly Gaussian. Some of that is to be expected when one averages *n* = 2^k^ different random values. If those had been iid values, then the CLT would apply and the averages would have skewness γ/√n and kurtosis κ/n, where γ and κ are the skewness and kurtosis of a single random integrand value. The *n* values in RQMC sampling are identically distributed but far from independent.

Figure 1 shows a selection of histograms, primarily the unusual ones. The upper left entry is a somewhat negatively skewed example from the RidgeJohnsonSU integrand. It is shown for the Sob-DS method, but that skewed pattern appeared for other methods too and not just for small *d* and *k*. The upper right entry shows the “spike plus outliers” distribution mentioned earlier. It also arises for nearly additive MC2 at somewhat higher values of *k* for each *d*. The phenomenon was pointed out by in [35] in one dimension and in [36] for multiple dimensions which includes some numerical illustrations. We expect the interval length to be quite variable in settings like this because, for example, *R* = 10 replicates might not include any of the outliers. The middle left entry shows a bimodal distribution that frequently appears for lattice samples with a baker transformation. It was also illustrated in [15] and in Chapter 17.4 of [33]. This bimodality did not always happen in our simulations. In more extreme settings, the distribution is nearly *U* shaped. The middle right entry shows a phenomenon for the IndSumNormal integrand. It is roughly
<figure><img src="image_3.png" /><figcaption>Figure 1: Some histograms of $N = 10,000$ integral estimates. In each case $L$ and $U$ mark the lower and upper exteme sample values and a red vertical line marks the integrand value 0.</figcaption></figure>

and coverage probability versus kurtosis and skewness for the case $R = 10$. There we see that kurtosis brings above nominal coverage for the standard $t$ intervals. We see coverage length decreasing with extreme kurtosis. That is what happens for Sob-LMS with SumUeU and MC2 where the RQMC pool is made up of spikes plus outliers. See the upper right histogram in Figure 1. For small $R$, the method only rarely gets an outlier, so most of those confidence intervals are extremely short. It is then a pleasant outcome that it covers the true mean often enough. We are cautious about expecting that to generalize beyond our examples.
<figure><img src="image_2.png" /><figcaption>Figure 1: Some histograms of $N = 10,000$ integral estimates, with $k = \log_2(n)$. In each case $L$ and $U$ mark the lower and upper exteme sample values and a red vertical line marks the integrand value 0.</figcaption></figure>

Gaussian plus a spike near one value. For smaller values of $k$, this integrand gives a discrete distribution with a roughly binomial appearance. The lower left entry is a roughly Gaussian histogram, similar to most of those in the data set. The lower right entry is a histogram for Sob-NUS that looks a lot like many that arise for Sob-LMS but which we have never before seen for Sob-NUS.

Some of our examples showed very extreme kurtosis. None of them showed extreme skewness. The standard confidence interval is known to have robust coverage in response to kurtosis at the expense of having very long confidence intervals, while being vulnerable to skewness. Figure 2 shows a summary of interval length
<figure><img src="image_4.png" /><figcaption>Figure 2: These figures describe coverage and length of normal theory confidence intervals for *R* = 10 replicates. Vertical kurtosis lines are at 0 (for the normal distribution). Horizontal reference lines are at the nominal level as well as at levels where we conclude that coverage is outside the interval [0.94, 0.96].</figcaption></figure>

# 5 CONCLUSION

In our examples, we saw that plain normal theory confidence intervals for RQMC performed best overall. Our criterion was simply whether the two-sided interval covered well. This was somewhat surprising as the bootstrap *t* method had much better coverage in the simulations reported in [2] and [28], even for numbers of replications as low as 5. We can understand this by noticing that the RQMC estimates we considered never had extreme skewness. That would happen if a CLT applied to them. However, the CLT does not always apply and several of the histograms are strongly non-Gaussian while still nearly symmetric. A similar phenomenon of symmetric but non-Gaussian limiting distributions arises in a strategy that recycles physically generated random numbers [32]. Paper [36] show that the RQMC error in Sob-LMS is a randomly weighted sum of the average of some Walsh functions. Those Walsh function averages have a symmetric distribution, but because the random weights are dependent, it does not automatically mean that
<figure><img src="image_3.png" /><figcaption>Figure 2: These figures describe coverage and length of normal theory confidence intervals for *R* = 10 replicates. Vertical kurtosis lines are at 0 (for the normal distribution). Horizontal reference lines are at the nominal level as well as at levels where we conclude that coverage is outside the interval [0.94, 0.96].</figcaption></figure>

and coverage probability versus kurtosis and skewness for the case $R = 10$. There we see that kurtosis brings above nominal coverage for the standard *t* intervals. We see interval length decreasing with extreme kurtosis. That is what happens for Sob-LMS with SumUeU and MC2 where the RQMC pool is made up of spikes plus outliers. See the upper right histogram in Figure 1. For small *R*, the method only rarely gets an outlier, so most of those confidence intervals are extremely short. It is then a pleasant outcome that it covers the true mean often enough. We are cautious about expecting that to generalize beyond our examples.

# 5 CONCLUSION

In our examples, we saw that plain normal theory confidence intervals for RQMC performed best overall. Our criterion was simply whether the two-sided interval covered well. This was somewhat surprising as the bootstrap *t* method had much
the RQMC error has a symmetric distribution.

The case of SumNormalInd highlights an issue with the bootstrap *t* not working on discrete random variables. Discrete random variables supported on an arithmetic sequence are not included in the theory of [8]. A distribution that is a mixture of a discrete and a continuous component will also face difficulties due to the discrete component yielding more frequent *s*\* = 0 when *R* is small.

The standard normal theory intervals are known to underperform the bootstrap
*t* asymptotically in the case where one-sided intervals are desired. There the cov-
erage error is *O*(1/√*n*) for the standard method and *O*(1/*n*) for the bootstrap *t*. If
however, symmetry is a ubiquitous property of RQMC estimates, then that advan-
tage of the bootstrap *t* disappears.

We have avoided running any examples that include rare events. For those
it is necessary to employ a method such as importance sampling in order to get
reasonable accuracy from RQMC. It is possible that rare-event settings without
importance sampling will produce very skewed histograms of RQMC points.

ACKNOWLEDGMENTS

The work of P. L'Ecuyer was supported by a Discovery Grant RGPIN-2018-05795 from NSERC-Canada. The work of A. B. Owen was supported by NSF grants IIS-1837931 and DMS-2152780.

References

[1] R. R. Bahadur and L. J. Savage. The nonexistence of certain statistical procedures in nonparametric problems. *The Annals of Mathematical Statistics*, 27(4):1115–1122, 1956.

[2] D. Choquet, P. L'Ecuyer, and C. Léger. Bootstrap confidence intervals for ratios of expectations. *ACM Transactions on Modeling and Computer Simulation*, 9(4):326–348, 1999.

[3] R. Cranley and T. N. L. Patterson. Randomization of number theoretic methods for multiple integration. *SIAM Journal on Numerical Analysis*, 13(6):904–914, 1976.

[4] J. Dick, P. Kritzer, and F. Pillichshammer. *Lattice Rules: Numerical Integration, Approximation, and Discrepancy*. Springer Nature, 2022.
better coverage in the simulations reported in [2] and [28], even for numbers of
replications as low as 5. We can understand this by noticing that the RQMC es-
timates we considered never had extreme skewness. That would happen if a CLT
applied to them. However, the CLT does not always apply and several of the his-
tograms are strongly non-Gaussian while still nearly symmetric. A similar phe-
nomenon of symmetric but non-Gaussian limiting distributions arises in a strategy
that recycles physically generated random numbers [32]. It is shown in [36] that
the RQMC error in Sob-LMS is a randomly weighted sum of the average of some
Walsh functions. Those Walsh function averages have a symmetric distribution,
but because the random weights are dependent, it does not automatically mean that
the RQMC error has a symmetric distribution.

The case of SumNormalInd highlights an issue with the bootstrap *t* not working on discrete random variables. Discrete random variables supported on an arithmetic sequence are not included in the theory of [8]. A distribution that is a mixture of a discrete and a continuous component will also face difficulties due to the discrete component yielding more frequent *S*\* = 0 when *R* is small.

The standard normal theory intervals are known to underperform the bootstrap
*t* asymptotically in the case where one-sided intervals are desired. There the cov-
erage error is *O*(1/√*n*) for the standard method and *O*(1/*n*) for the bootstrap *t*. If
however, symmetry is a ubiquitous property of RQMC estimates, then that advan-
tage of the bootstrap *t* disappears.

We have avoided running any examples that include rare events. For those
it is necessary to employ a method such as importance sampling in order to get
reasonable accuracy from RQMC. It is possible that rare-event settings without
importance sampling will produce very skewed histograms of RQMC points.

ACKNOWLEDGMENTS

The work of P. L'Ecuyer was supported by a Discovery Grant RGPIN-2018-05795 from NSERC-Canada. The work of A. B. Owen was supported by NSF grants IIS-1837931 and DMS-2152780.

References

[1] R. R. Bahadur and L. J. Savage. The nonexistence of certain statistical procedures in nonpara-metric problems. *The Annals of Mathematical Statistics*, 27(4):1115–1122, 1956.

[2] D. Choquet, P. L'Ecuyer, and C. Léger. Bootstrap confidence intervals for ratios of expectations. *ACM Transactions on Modeling and Computer Simulation*, 9(4):326–348, 1999.

[3] R. Cranley and T. N. L. Patterson. Randomization of number theoretic methods for multiple integration. *SIAM Journal on Numerical Analysis*, 13(6):904–914, 1976.
[5] J. Dick and F. Pillichshammer. *Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration*. Cambridge University Press, Cambridge, U.K., 2010.

[6] B. Efron and R. J. Tibshirani. *An Introduction to the Bootstrap*. Chapman & Hall, New York, NY, 1994.

[7] Bradley Efron. Nonparametric standard errors and confidence intervals. *Canadian Journal of Statistics*, 9(2):139–158, 1981.

[8] P. Hall. Theoretical comparison of bootstrap confidence intervals. *The Annals of Statistics*, 16(3):927–953, 1988.

[9] P. Hall. *The Bootstrap and Edgeworth Expansions*. Springer, New York, 1992.

[10] F. J. Hickernell. Obtaining $O(N^{-2+\epsilon})$ convergence for lattice quadrature rules. In K.-T. Fang, F. J. Hickernell, and H. Niederreiter, editors, *Monte Carlo and Quasi-Monte Carlo Methods 2000*, pages 274–289, Berlin, 2002. Springer-Verlag.

[11] C. R. Hoyt and A. B. Owen. Mean dimension of ridge functions. *SIAM Journal on Numerical Analysis*, 58(2):1195–1216, 2020.

[12] S. Joe and F. Y. Kuo. Constructing Sobol sequences with better two-dimensional projections. *SIAM Journal on Scientific Computing*, 30(5):2635–2654, 2008.

[13] N. L. Johnson. Systems of frequency curves generated by methods of translation. *Biometrika*, 36(12):149–176, 1949.

[14] P. L'Ecuyer. Quasi-Monte Carlo methods with applications in finance. *Finance and Stochastics*, 13(3):307–349, 2009.

[15] P. L'Ecuyer. Randomized quasi-Monte Carlo: An introduction for practitioners. In P. W. Glynn and A. B. Owen, editors, *Monte Carlo and Quasi-Monte Carlo Methods: MCQMC 2016*, pages 29–52, Berlin, 2018. Springer.

[16] P. L'Ecuyer. SSJ: Stochastic simulation in Java. https://github.com/umontreal-simul/ssj, accessed 27th April 2023, 2023.

[17] P. L'Ecuyer and C. Lemieux. Variance reduction via lattice rules. *Management Science*, 46(9):1214–1235, 2000.
[4] J. Dick, P. Kritzer, and F. Pillichshammer. *Lattice Rules: Numerical Integration, Approximation, and Discrepancy*. Springer Nature, 2022.

[5] J. Dick and F. Pillichshammer. *Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration*. Cambridge University Press, Cambridge, U.K., 2010.

[6] B. Efron and R. J. Tibshirani. *An Introduction to the Bootstrap*. Chapman & Hall, New York, NY, 1994.

[7] Bradley Efron. Nonparametric standard errors and confidence intervals. *Canadian Journal of Statistics*, 9(2):139–158, 1981.

[8] P. Hall. Theoretical comparison of bootstrap confidence intervals. *The Annals of Statistics*, 16(3):927–953, 1988.

[9] P. Hall. *The Bootstrap and Edgeworth Expansions*. Springer, New York, 1992.

[10] F. J. Hickernell. Obtaining $O(N^{-2+\epsilon})$ convergence for lattice quadrature rules. In K.-T. Fang, F. J. Hickernell, and H. Niederreiter, editors, *Monte Carlo and Quasi-Monte Carlo Methods 2000*, pages 274–289, Berlin, 2002. Springer-Verlag.

[11] C. R. Hoyt and A. B. Owen. Mean dimension of ridge functions. *SIAM Journal on Numerical Analysis*, 58(2):1195–1216, 2020.

[12] S. Joe and F. Y. Kuo. Constructing Sobol sequences with better two-dimensional projections. *SIAM Journal on Scientific Computing*, 30(5):2635–2654, 2008.

[13] N. L. Johnson. Systems of frequency curves generated by methods of translation. *Biometrika*, 36(12):149–176, 1949.

[14] P. L'Ecuyer. Quasi-Monte Carlo methods with applications in finance. *Finance and Stochastics*, 13(3):307–349, 2009.

[15] P. L'Ecuyer. Randomized quasi-Monte Carlo: An introduction for practitioners. In P. W. Glynn and A. B. Owen, editors, *Monte Carlo and Quasi-Monte Carlo Methods: MCQMC 2016*, pages 29–52, Berlin, 2018. Springer.

[16] P. L'Ecuyer. SSJ: Stochastic simulation in Java. https://github.com/umontreal-simul/ssj, accessed 27th April 2023, 2023.

[17] P. L'Ecuyer and C. Lemieux. Variance reduction via lattice rules. *Management Science*, 46(9):1214–1235, 2000.

[18] P. L'Ecuyer and C. Lemieux. Recent advances in randomized quasi-Monte Carlo methods. In M. Dror, P. L'Ecuyer, and F. Szidarovszky, editors, *Modeling Uncertainty: An Examination of Stochastic Theory, Methods, and Applications*, pages 419–474. Kluwer Academic, Boston, 2002.

[19] P. L'Ecuyer, P. Marion, M. Godin, and F. Puchhammer. A tool for custom construction of QMC and RQMC point sets. In A. Keller, editor, *Monte Carlo and Quasi-Monte Carlo Methods: MC-QMC 2020*, pages 51–70, Berlin, 2022. Springer. https://arxiv.org/abs/2012.10263.

[20] P. L'Ecuyer and D. Munger. Algorithm 958: Lattice builder: A general software tool for constructing rank-1 lattice rules. *ACM Transactions on Mathematical Software*, 42(2):Article 15, 2016.

[21] P. L'Ecuyer, D. Munger, and B. Tuffin. On the distribution of integration error by randomly-shifted lattice rules. *Electronic Journal of Statistics*, 4:950–993, 2010.

[22] C. Lemieux. *Monte Carlo and Quasi-Monte Carlo Sampling*. Series in Statistics. Springer, New York, 2009.
[18] P. L'Ecuyer and C. Lemieux. Recent advances in randomized quasi-Monte Carlo methods. In M. Dror, P. L'Ecuyer, and F. Szidarovszky, editors, *Modeling Uncertainty: An Examination of Stochastic Theory, Methods, and Applications*, pages 419-474. Kluwer Academic, Boston, 2002.

[19] P. L'Ecuyer, P. Marion, M. Godin, and F. Puchhammer. A tool for custom construction of QMC and RQMC point sets. In A. Keller, editor, *Monte Carlo and Quasi-Monte Carlo Methods: MCQMC 2020*, pages 51-70, Berlin, 2022. Springer. https://arxiv.org/abs/2012.10263.

[20] P. L'Ecuyer and D. Munger. Algorithm 958: Lattice builder: A general software tool for constructing rank-1 lattice rules. *ACM Transactions on Mathematical Software*, 42(2):Article 15, 2016.

[21] P. L'Ecuyer, D. Munger, and B. Tuffin. On the distribution of integration error by randomly-shifted lattice rules. *Electronic Journal of Statistics*, 4:950-993, 2010.

[22] C. Lemieux. *Monte Carlo and Quasi-Monte Carlo Sampling*. Series in Statistics. Springer, New York, 2009.

[23] W.-L. Loh. On the asymptotic distribution of scrambled net quadrature. *Annals of Statistics*, 31(4):1282-1324, 2003.

[24] J. Matousěk. On the $L_2$-discrepancy for anchored boxes. *J. of Complexity*, 14:527-556, 1998.

[25] W. J. Morokoff and R. E. Caflisch. Quasi-Monte Carlo integration. *Journal of Computational Physics*, 122:218-230, 1995.

[26] M. K. Nakayama and B. Tuffin. Sufficient conditions for central limit theorems and confidence intervals for randomized quasi-Monte Carlo methods. techreport hal-03196085, INRIA, 2021. https://hal.inria.fr/hal-03196085.

[27] H. Niederreiter. *Random Number Generation and Quasi-Monte Carlo Methods*, volume 63. SIAM, Philadelphia, 1992.

[28] A. B. Owen. Small sample central confidence intervals for the mean. Technical Report 302, Stanford University, Department of Statistics, 1988. https://purl.stanford.edu/mz765np4744.

[29] A. B. Owen. Randomly permuted (t,m,s)-nets and (t,s)-sequences. In *Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing: Lecture Notes in Statistics*, volume 106, pages 299-317. Springer, 1995.
[23] W.-L. Loh. On the asymptotic distribution of scrambled net quadrature. *Annals of Statistics*, 31(4):1282–1324, 2003.

[24] J. Matousšk. On the $L_2$-discrepancy for anchored boxes. *J. of Complexity*, 14:527–556, 1998.

[25] W. J. Morokoff and R. E. Caflisch. Quasi-Monte Carlo integration. *Journal of Computational Physics*, 122:218–230, 1995.

[26] M. K. Nakayama and B. Tuffin. Sufficient conditions for central limit theorems and confidence intervals for randomized quasi-Monte Carlo methods. techreport hal-03196085, INRIA, 2021. https://hal.inria.fr/hal-03196085.

[27] H. Niederreiter. *Random Number Generation and Quasi-Monte Carlo Methods*, volume 63. SIAM, Philadelphia, 1992.

[28] A. B. Owen. Small sample central confidence intervals for the mean. Technical Report 302, Stanford University, Department of Statistics, 1988. https://purl.stanford.edu/mz765np4744.

[29] A. B. Owen. Randomly permuted (t,m,s)-nets and (t,s)-sequences. In *Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing: Lecture Notes in Statistics*, volume 106, pages 299–317. Springer, 1995.

[30] A. B. Owen. Monte Carlo variance of scrambled net quadrature. *SIAM Journal of Numerical Analysis*, 34:1884–1910, 1997.

[31] A. B. Owen. Variance with alternative scramblings of digital nets. *ACM Transactions on Modeling and Computer Simulation*, 13(4):363–378, 2003.

[32] A. B. Owen. Recycling physical random numbers. *Electronic Journal of Statistics*, 3:1531–1541, 2009.

[33] A. B. Owen. *Practical Quasi-Monte Carlo*. Draft available at https://artowen.su doma in gs /mc /practicalqmc.pdf, 2023.

[34] Art B Owen and Daniel Rudolf. A strong law of large numbers for scrambled net integration. *SIAM Review*, 63(2):360–372, 2021.

[35] Zexin Pan and A. B. Owen. Super-polynomial accuracy of one dimensional randomized nets using the median of means. *Mathematics of Computation*, 92(340):805–837, 2023.

[36] Zexin Pan and Art B Owen. Super-polynomial accuracy of multidimensional randomized nets using the median-of-means. Technical report, arXiv:2208.05078, 2022.

[37] I. H. Sloan and H. Woźniakowski. When are quasi-Monte Carlo algorithms efficient for high dimensional integration? *Journal of Complexity*, 14:1–33, 1998.

[38] B. Tuffin. Variance reduction order using good lattice points in Monte Carlo methods. *Computing*, 61(4):371–378, 1998.

[39] B. Tuffin. Randomization of quasi-monte carlo methods for error estimation: Survey and normal approximation. *Monte Carlo Methods and Applications*, 10(3-4):617–628, 2004.

[40] U.S. Nuclear Regulatory Commission. Guidelines for evaluating fatigue analyses incorporating the life reduction of metal components due to the effects of the light-water reactor environment for new reactors. U.S. Nuclear Regulatory Commission Regulatory Guide 1.207, U.S. Nuclear Regulatory Commission, Washington, DC, 2007.

[41] R.-X. Yue and F. J. Hickernell. The discrepancy and gain coefficients of scrambled digital nets. *Journal of Complexity*, 18(1):135–151, 2002.
[30] A. B. Owen. Monte Carlo variance of scrambled net quadrature. *SIAM Journal of Numerical Analysis*, 34:1884–1910, 1997.

[31] A. B. Owen. Variance with alternative scramblings of digital nets. *ACM Transactions on Modeling and Computer Simulation*, 13(4):363–378, 2003.

[32] A. B. Owen. Recycling physical random numbers. *Electronic Journal of Statistics*, 3:1531–1541, 2009.

[33] A. B. Owen. *Practical Quasi-Monte Carlo*. Draft available at https://artowen.su doma in s/mc/practicalqmc.pdf, 2023.

[34] Art B Owen and Daniel Rudolf. A strong law of large numbers for scrambled net integration. *SIAM Review*, 63(2):360–372, 2021.

[35] Zexin Pan and A. B. Owen. Super-polynomial accuracy of one dimensional randomized nets using the median of means. *Mathematics of Computation*, 92(340):805–837, 2023.

[36] Zexin Pan and Art B Owen. Super-polynomial accuracy of multidimensional randomized nets using the median-of-means. Technical report, arXiv:2208.05078, 2022.

[37] I. H. Sloan and H. Woźniakowski. When are quasi-Monte Carlo algorithms efficient for high dimensional integration? *Journal of Complexity*, 14:1–33, 1998.

[38] B. Tuffin. Variance reduction order using good lattice points in Monte Carlo methods. *Computing*, 61(4):371–378, 1998.

[39] B. Tuffin. Randomization of quasi-monte carlo methods for error estimation: Survey and normal approximation. *Monte Carlo Methods and Applications*, 10(3-4):617–628, 2004.

[40] U.S. Nuclear Regulatory Commission. Guidelines for evaluating fatigue analyses incorporating the life reduction of metal components due to the effects of the light-water reactor environment for new reactors. U.S. Nuclear Regulatory Commission Regulatory Guide 1.207, U.S. Nuclear Regulatory Commission, Washington, DC, 2007.

[41] R.-X. Yue and F. J. Hickernell. The discrepancy and gain coefficients of scrambled digital nets. *Journal of Complexity*, 18(1):135–151, 2002.
## AUTHOR BIOGRAPHIES

**PIERRE L'ECUYER** is a Professor in the Département d'Informatique et de Recherche Opérationnelle, Université de Montréal, Canada. He is a member of the CIRRELT and GERAD. His main research interests are random number generation, quasi-Monte Carlo methods, efficiency improvement via variance reduction, sensitivity analysis and optimization of discrete-event stochastic systems. He has published over 300 scientific articles and has developed software libraries and systems for random number generation and stochastic simulation. Web page: http://www.iro.umontreal.ca/~lecuyer. Email: lecuyer@iro.umontreal.ca.

**MARVIN K. NAKAYAMA** is a professor in the Department of Computer Science at the New Jersey Institute of Technology. He received an M.S. and Ph.D. in operations research from Stanford University and a B.A. in mathematics-computer science from U.C. San Diego. He is an associate editor for *ACM Transactions on Modeling and Computer Simulation*, and served as the simulation area editor for the *INFORMS Journal on Computing* from 2007–2016. His research interests include simulation, modeling, statistics, risk analysis, and energy. His email: marvin@njit.edu.

**ART B. OWEN** is the Max H. Stein Professor of Statistics at Stanford University. His research interests include quasi-Monte Carlo sampling, nonparametric confidence intervals, and variable importance measures. Web page: https://artowen.su.domains/. Email: owen@stanford.edu.

**BRUNO TUFFIN** received his PhD degree in applied mathematics from the University of Rennes 1, France, in 1997. Since then, he has been with Inria in Rennes. His research interests include developing Monte Carlo and quasi-Monte Carlo simulation techniques for the performance evaluation of telecommunication systems and telecommunication-related economical models. He is currently Area Editor for *INFORMS Journal on Computing* and Associate Editor for *ACM Transactions on Modeling and Computer Simulation* and *Queuing Systems*. His email: bruno.tuffin@inria.fr.