Asymptotic Normality of Extensible Grid Sampling

TL;DR

This paper proves the asymptotic normality of Hilbert space filling curve-based estimates in numerical integration, extending previous results to discontinuous functions and providing variance bounds.

Contribution

It establishes the asymptotic normality of HSFC-based estimates for a broader class of functions, including discontinuous ones, and derives variance lower bounds.

Findings

Asymptotic normality holds for functions in C^1 and certain discontinuous functions.

Variance of the estimate has a lower bound of order n^{-1-2/d}.

Results extend previous normality findings from smooth functions to more general cases.

Abstract

Recently, He and Owen (2016) proposed the use of Hilbert's space filling curve (HSFC) in numerical integration as a way of reducing the dimension from $d>1$ to $d=1$. This paper studies the asymptotic normality of the HSFC-based estimate when using scrambled van der Corput sequence as input. We show that the estimate has an asymptotic normal distribution for functions in $C^1([0,1]^d)$, excluding the trivial case of constant functions. The asymptotic normality also holds for discontinuous functions under mild conditions. It was previously known only that scrambled $(0,m,d)$-net quadratures enjoy the asymptotic normality for smooth enough functions, whose mixed partial gradients satisfy a H\"older condition. As a by-product, we find lower bounds for the variance of the HSFC-based estimate. Particularly, for nontrivial functions in $C^1([0,1]^d)$, the low bound is of order $n^{-1-2/d}$, which matches the rate of the upper bound established in He and Owen (2016).

Equations

$$\mu ={\int }_{[0,1{]}^d}f(X)\mathrm{d}X.$$

$${\hat{\mu }}_n=\frac{1}{n}\sum _{i=1}^{n}f(H(x_i)),$$

$$\frac{{\hat{\mu }}_n-\mu }{\sqrt{\mathrm{Var}({\hat{\mu }}_n)}}\to N(0,1)$$

$${\tilde{\mu }}_n=\frac{1}{n}\sum _{i=1}^{n}f(U_i),$$

$$\underset{n\to \infty }{\lim }{n}^{1+\frac{2}{d}}\mathrm{Var}({\tilde{\mu }}_n)=\frac{1}{12}{\int }_{[0,1{]}^d}\Vert \nabla f(X){\Vert }^2\mathrm{d}X,$$

$$f(U_i)=L_i+R_i,$$

$$\mathrm{Var}(L_i)=\frac{1}{12m^2}\sum _{i=1}^d{\left(\frac{\partial }{\partial X_i}f(c_i)\right)}^2=\frac{1}{12m^2}\Vert \nabla f(c_i){\Vert }^2,$$

$$\mathrm{Var}({\tilde{\mu }}_n)$$

$$=\frac{1}{n^2}\sum _{i=1}^n\mathrm{Var}(L_i)+o\left(\frac{1}{n{m}^2}\right)$$

$$=\frac{1}{12n^2{m}^2}\sum _{i=1}^n\Vert \nabla f(c_i){\Vert }^2+o\left(\frac{1}{n{m}^2}\right).$$

$$\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{i=1}^n\Vert \nabla f(c_i){\Vert }^2={\int }_{[0,1{]}^d}\Vert \nabla f(X){\Vert }^2\mathrm{d}X,$$

$$\frac{{\tilde{\mu }}_n-\mu }{\sigma n^{-1/2-1/d}}\to N(0,1),$$

$$\mathbb{E}[|f(U_i)-\mathbb{E}[f(U_i)]{|}^{2+\delta }]\le C_{d,\delta }{n}^{-\frac{2+\delta }{d}},$$

$$\underset{n\to \infty }{\lim }\frac{s_n^2}{n^{1-\frac{2}{d}}}=\frac{1}{12}{\int }_{[0,1{]}^d}\Vert \nabla f(X){\Vert }^2\mathrm{d}X>0.$$

$$\underset{n\to \infty }{\lim }\frac{1}{s_n^{2+\delta }}\sum _{i=1}^n\mathbb{E}[|f(U_i)-\mathbb{E}[f(U_i)]{|}^{2+\delta }]$$

$$=\underset{n\to \infty }{\mathrm{limsup}}{C}_{d,\delta }{n}^{-\frac{\delta }{2}}{\left(\frac{n^{1-\frac{2}{d}}}{s_n^2}\right)}^{\frac{2+\delta }{2}}=0$$

$$a_i=\sum _{j=1}^{\infty }{a}_{ij}{b}^{-j}.$$

$$x_{ij}={\pi }_{a_{i1}\dots a_{ij-1}}(a_{ij}).$$

$${\hat{\mu }}_n=\frac{1}{n}\sum _{i=1}^{n}f(X^{(i)}),$$

$$\mathrm{Var}({\hat{\mu }}_n)\ge \frac{{\sigma }^2}{96}{2}^{-2/d-d}{n}^{-1-2/d}.$$

$$\frac{{\hat{\mu }}_n-\mu }{\sqrt{\mathrm{Var}({\hat{\mu }}_n)}}\to N(0,1),$$

$${\sigma }_i^2$$

$$\ge \frac{1}{{\lambda }_d(E_i)}{\int }_{{\tilde{E}}_i}[f(X)-{\mu }_i{]}^2\mathrm{d}X$$

$$=\frac{1}{{\lambda }_d(E_i)}{\int }_{{\tilde{E}}_i}\left([f(X)-{\mu }_i^{\prime }{]}^2+2(f(X)-{\mu }_i^{\prime })({\mu }_i^{\prime }-{\mu }_i)+({\mu }_i^{\prime }-{\mu }_i{)}^2\right)\mathrm{d}X$$

$$=\frac{{\lambda }_d({\tilde{E}}_i)}{{\lambda }_d(E_i)}\left(\widetilde{\mathrm{Var}}(f(X^{(i)}))+({\mu }_i^{\prime }-{\mu }_i{)}^2\right)$$

$$\ge 2^{-(1+d)}\widetilde{\mathrm{Var}}(f(X^{(i)})),$$

$${\lambda }_d({\tilde{E}}_i)={\lambda }_1(J_i)=\frac{1}{2^{dm}}\ge \frac{1}{2^{1+d}n}=2^{-(1+d)}{\lambda }_d(E_i).$$

$$\frac{s_n^2}{n^{1-\frac{2}{d}}}\ge \frac{1}{2^{1+d}{n}^{1-\frac{2}{d}}}\sum _{i=1}^n\widetilde{\mathrm{Var}}(f(X^{(i)}))=:K_1(n).$$

$$\widetilde{\mathrm{Var}}(f(X^{(i)}))=\frac{1}{12\cdot 2^{2m}}\Vert \nabla f(c_i){\Vert }^2+o(2^{-2m}),$$

$$\underset{n\to \infty }{\mathrm{liminf}}{K}_1(n)$$

$$\ge \underset{n\to \infty }{\mathrm{liminf}}\frac{1}{2^{1+d}{n}^{1-2/d}}\frac{1}{12\cdot 2^{2(1+d)/d}{n}^{2/d}}\sum _{i=1}^n\Vert \nabla f(c_i){\Vert }^2$$

$$=\frac{1}{96\cdot 2^{2/d+d}}{\int }_{[0,1{]}^d}\Vert \nabla f(X){\Vert }^2\mathrm{d}X>0.$$

$$\underset{n\to \infty }{\mathrm{liminf}}\frac{s_n^2}{n^{1-\frac{2}{d}}}\ge \underset{n\to \infty }{\mathrm{liminf}}{K}_1(n)\ge \frac{1}{96\cdot 2^{2/d+d}}{\int }_{[0,1{]}^d}\Vert \nabla f(X){\Vert }^2\mathrm{d}X>0.$$

$$\mathbb{E}\left[{\left|f(X^{(i)})-\mathbb{E}[f(X^{(i)})]\right|}^{2+\delta }\right]\le C_{d,\delta }{n}^{-\frac{2+\delta }{d}},$$

$$\underset{n\to \infty }{\lim }\frac{1}{s_n^{2+\delta }}\sum _{i=1}^n\mathbb{E}\left[{\left|f(X^{(i)})-\mathbb{E}[f(X^{(i)})]\right|}^{2+\delta }\right]$$

Full text

Keywords: asymptotic normality; Hilbert’s space filling curve; van der Corput sequence; randomized quasi-Monte Carlo; extensible grid sampling

1 Introduction

Quasi-Monte Carlo (QMC) sampling has gained increasing popularity in numerical integration over the unit cube $[0,1]^{d}$ (see, e.g., L’Ecuyer (2009); Dick and Pillichshammer (2010); Dick et al. (2013)). It is known that functions with finite variation in the sense of Hardy and Krause can be integrated with an error of $O(n^{-1}(\log n)^{d})$, compared to $O(n^{-1/2})$ for ordinary Monte Carlo sampling; see Niederreiter (1992) for details.

In this paper, we consider an alternative numerical integration based on Hilbert’s space filling curve (HSFC) as introduced in He and Owen (2016). An HSFC is a continuous mapping $H(x)$ from $[0,1]$ to $[0,1]^{d}$ for $d\geq 1$. We take the convention that $H(x)=x$ for $d=1$. Formally, the Hilbert curve is defined through the limit of a series of recursive curves. An illustration of the generative process of the HSFC with increasing recursion order for $d=2$ is presented in Figure 1. For detailed definitions and the properties of the HSFC, we refer to He and Owen (2016). Let us consider the problem of estimating an integral over the $d$-dimensional unit cube $[0,1]^{d}$:

[TABLE]

The HSFC-based estimate takes the form

[TABLE]

where $x_{i}$ are some well-chosen points in $[0,1]$. In this paper, we focus on the case of using van der Corput sequence (in base $b\geq 2$) with the nested uniform scrambling of Owen (1995) as the inputs $x_{i}$, for which the estimate is extensible and unbiased. The nested uniform scrambling method is a kind of randomization techniques used commonly in randomized QMC; see Owen (1995) for details and L’Ecuyer and Lemieux (2002) for a survey of various randomized QMC methods. Matoušek (1998) proposed a random linear scrambling method that does not require as much randomness and storage. He and Owen (2016) found convergence rates of the extensible estimate for functions that are Lipschitz continuous or piecewise Lipschitz continuous. More precisely, for Lipschitz continuous functions, they derived a root mean-squared error (RMSE) of $O(n^{-1/2-1/d})$. For the piecewise Lipschitz continuous functions, an RMSE of $O(n^{-1/2-1/(2d)})$ is obtained. Schretter et al. (2016) compared the star discrepancies and RMSEs of using the van der Corput and the golden ratio generator sequences.

Actually, the upper bounds of the RMSE do not tell much about the error distribution. It is often of interest to obtain asymptotically valid confidence interval type guarantees as the usual Monte Carlo sampling. The central limit theorem (CLT) is invoked routinely to compute a confidence interval on the estimate based on a normal approximation. Indeed, Loh (2003) showed that the nested uniform scrambled $(0,m,d)$-net (in base $b$) estimate has an asymptotic normal distribution for smooth enough functions. Recently, by developing on the work by Loh (2003), Basu and Mukherjee (2016) showed that the scrambled geometric net estimate has an asymptotic normal distribution for certain smooth functions defined on products of suitable subsets of $\Re^{d}$. However, in most cases, the randomly-shifted lattice rule (another branch of randomized QMC techniques) estimates may be far from the normally distributed ones; see L’Ecuyer et al. (2010) for discussions and examples.

In this paper we study the asymptotic normality of the HSFC-based estimate (1.2) with sample sizes $n=b^{m}$, $m=0,1,2,\dots$. In practice, we often choose $b=2$ because the base $b=2$ is the same as used to approximate (and define) the Hilbert curve; see Butz (1971) for the algorithm. The main contribution of our paper is two fold. First, for nontrivial functions in $C^{1}([0,1]^{d})$, we establish a lower bound on $\mathrm{Var}(\hat{\mu}_{n})$ which matches the upper bound $O(n^{-1-2/d})$ found in He and Owen (2016). A similar lower bound is established for piecewise smooth functions. Second, we prove that the asymptotic normality of the HSFC-based estimate $\hat{\mu}_{n}$ holds for three classes of functions. In other words, we show that

[TABLE]

in distribution as $n=b^{m}\to\infty$. These results can be applied to stratified sampling on a regular grid with sample sizes $n=m^{d}$, but the HSFC-based estimate we study, does not require the highly composite sample sizes that the grid sampling requires, particularly for large $d$. The main idea to prove the asymptotic normality is based on the Lyapunov CLT (see, e.g., Chung (2001)), which is quite different from the techniques used in Loh (2003) and Basu and Mukherjee (2016). Our proofs do not rely on the upper bounds established in He and Owen (2016).

The rest of the paper is organized as follows. In Section 2, we study the asymptotic normality of stratified sampling on a regular grid with $n=m^{d}$, which can be viewed as a special case of the HSFC-based estimate (1.2) using scrambled van der Corput sequence. In Section 3, we give lower bounds on $\mathrm{Var}(\hat{\mu}_{n})$ and establish asymptotic normality of the estimate $\hat{\mu}_{n}$ for three cases of integrands. In Section 4, we give some empirical verifications on asymptotic normality for the HSFC sampling and other competitive QMC methods. Section 5 concludes this paper.

2 Grid-based Stratified Sampling

In this section, we consider the regular grid sampling with sample sizes $n=m^{d}$. The $d$-dimensional unit cube $[0,1]^{d}$ can be split into $m^{d}$ congruent subcubes with sides of length $1/m$, say, $E_{i}$, $i=1,\dots,n$. The grid-based stratified estimate of the integral (1.1) is given by

[TABLE]

where $U_{i}\sim\mathrm{Unif}(E_{i})$ independently. Denote $\nabla f(X)=(\frac{\partial f(X)}{\partial X_{1}},\dots,\frac{\partial f(X)}{\partial X_{d}})^{\top}$ as the gradient vector of $f(X)$, and let $\|\cdot\|$ be the usual Euclidean norm. The next lemma discusses the variance of $\tilde{\mu}_{n}$, which was proved in Owen (2013). We prove it here also, because we make extensive use of that result.

Lemma 1**.**

Assume that $f(X)\in C^{1}([0,1]^{d})$. Then

[TABLE]

where the limit is taken through values $n=m^{d}$ as $m\to\infty$.

Proof.

Note that $U_{i}$ is uniformly distributed within the cube $E_{i}$ with sides of length $1/m$. Let $c_{i}$ be the center of $E_{i}$. Since $f\in C^{1}([0,1]^{d})$, the first-order Taylor approximation gives

[TABLE]

where $L_{i}=f(c_{i})+\nabla f(c_{i})^{\top}(U_{i}-c_{i})$ and $R_{i}=o(1/m)$. For the linear term $L_{i}$, we have

[TABLE]

since $U_{i}-c_{i}\sim\mathrm{Unif}[-1/(2m),1/(2m)]^{d}$. For the error term $R_{i}$, we have $\mathrm{Var}(R_{i})=o(1/m^{2})$. Also, $\mathrm{Cov}(L_{i},R_{i})=o(1/m^{2})$. As a result,

[TABLE]

Since $m=n^{1/d}$ and

[TABLE]

we conclude that (2.2) holds. ∎

Theorem 2**.**

If $f(X)\in C^{1}([0,1]^{d})$ and $\sigma^{2}=(1/12)\int_{[0,1]^{d}}\|\nabla f(X)\|^{2}\mathrm{d}X>0$, then

[TABLE]

in distribution as $n=m^{d}\to\infty$.

Proof.

Since $f\in C^{1}([0,1]^{d})$, $f$ is Lipschitz continuous. Then, for any $\delta>0$,

[TABLE]

where $C_{d,\delta}>0$ is some constant that only depends on $d$ and $\delta$, and we used the fact that the diameter of $E_{i}$ is $\sqrt{d}n^{-1/d}$.

Let $s_{n}^{2}=\sum_{i=1}^{n}\sigma_{i}^{2}$, where $\sigma_{i}^{2}=\mathrm{Var}(f(U_{i}))$. From Lemma 1, we have

[TABLE]

Therefore, the Lyapunov condition

[TABLE]

is satisfied. Using the Lyapunov CLT, we get (2.4). ∎

To avoid the trivial case of constant functions (which yields an identically zero variance), Theorem 2 assumes that $\int_{[0,1]^{d}}\|\nabla f(x)\|^{2}\mathrm{d}x>0$. Theorem 2 admits the asymptotic normality of $\tilde{\mu}_{n}$. The grid-based stratified estimate has variance $O(n^{-1-2/d})$, compared to Monte Carlo variance $O(n^{-1})$. This actually holds for Lipschitz continuous functions covering the class of functions $C^{1}([0,1]^{d})$ in Theorem 2.

3 HSFC-based Sampling

In this section, we study the HSFC-based estimate given by (1.2), where $x_{i}$ are the first $n=b^{m}$ points of the scrambled van der Corput sequence in base $b\geq 2$. Let $a_{i}$ be the first $n$ points of the van der Corput sequence in base $b$ van der Corput (1935). The integer $i-1\geq 0$ is written in base $b$ as $i-1=\sum_{j=1}^{\infty}a_{ij}b^{j-1}$ for $a_{ij}\in\{0,\dots,b-1\}$. Then $a_{i}$ is then defined by

[TABLE]

The scrambled version of $a_{1},\dots,a_{n}$ is $x_{1},\dots,x_{n}$ written as $x_{i}=\sum_{j=1}^{\infty}x_{ij}b^{-j}$, where $x_{ij}$ are defined through random permutations of the $a_{ij}$. These permutations depend on $a_{ik}$, for $k<j$. More precisely, $x_{i1}=\pi(a_{i1})$, $x_{i2}=\pi_{a_{i1}}(a_{i2})$ and generally for $j\geq 2$

[TABLE]

Each random permutation is uniformly distributed over the $b!$ permutations of $\{0,\dots,b-1\}$, and the permutations are mutually independent.

In this setting, thanks to the nice property of the nested uniform scrambling, the data values in the scrambled sequence can be reordered such that $x_{i}\sim\mathrm{Unif}(I_{i})$ independently with $I_{i}$ for $i=1,\dots,b^{m}$. Let $E_{i}=H(I_{i})$. As used in He and Owen (2016), the estimate (1.2) can be rewritten as

[TABLE]

where $X^{(i)}=H(x_{i})\sim\mathrm{Unif}(E_{i})$. This implies that the HSFC-based sampling is actually a stratified sampling because $\{E_{i}\}_{i=1}^{n}$ is a split of $[0,1]^{d}$. Figure 2 illustrates such splits of $[0,1]^{2}$ when $b=2$. He and Owen (2016) proved the unbiasedness of $\hat{\mu}_{n}$ for any $f\in L^{2}([0,1]^{d})$ and gave some upper bounds for $\mathrm{Var}(\hat{\mu}_{n})$ under certain assumptions on the class of integrands $f$. Their proofs make use of the properties of the HSFC presented in the next lemma, which are also important in studying the asymptotic normality of the HSFC-based sampling. Denote $\lambda_{d}(\cdot)$ as the Lebesgue measure on $\Re^{d}$.

Lemma 3**.**

Let $A=H([p,q])$ for $0\leq p<q\leq 1$. Then $\lambda_{d}(A)=\lambda_{1}([p,q])=q-p$. If $x\sim\mathrm{Unif}([p,q])$, then $H(x)\sim\mathrm{Unif}(A)$. Let $r$ be the diameter of $A$. Then $r\leq 2\sqrt{d+3}(q-p)^{1/d}$.

3.1 Smooth Functions

Loh (2003) and Basu and Mukherjee (2016) focused on smooth functions whose mixed partial gradient satisfies a Hölder condition, which was first studied in Owen (1997). Here we work with a weaker smoothness condition, in the sense that $f(X)\in C^{1}([0,1]^{d})$ as required in Theorem 2.

Theorem 4**.**

Assume that $f(X)\in C^{1}([0,1]^{d})$ and $\sigma^{2}=\int_{[0,1]^{d}}\|\nabla f(X)\|^{2}\mathrm{d}X>0$. Then for all sufficiently large $n$, we have

[TABLE]

Also,

[TABLE]

in distribution as $n=b^{m}\to\infty$.

Proof.

Let $m=\lceil\frac{\log_{2}n+1}{d}\rceil$. Then there exists an interval $J_{i}$ of the form $[(k-1)/2^{dm},k/2^{dm}]$ such that $J_{i}\subset I_{i}$. This is because $\lambda_{1}(I_{i})\geq 2\lambda_{1}(J_{i})$. Let $s_{n}^{2}=\sum_{i=1}^{n}\sigma_{i}^{2}$, where $\sigma_{i}^{2}=\mathrm{Var}(f(X^{(i)}))$. Let $\mu_{i}=\mathrm{E}(f(X^{(i)}))$, and let $\tilde{E}_{i}=H(J_{i})$. Note that $\tilde{E}_{i}\subset E_{i}$ due to $J_{i}\subset I_{i}$. Let $\mu^{\prime}_{i}=\tilde{\mathrm{E}}(f(X^{(i)}))$, where the expectation is taken from $X^{(i)}\sim\mathrm{Unif}(\tilde{\mathrm{E}}_{i})$. Based on some basic algebra, we find

[TABLE]

where $\widetilde{\mathrm{Var}}$ is taken over $X^{(i)}\sim\mathrm{Unif}(\tilde{E}_{i})$, and we used two results by applying Lemma 3 that $\lambda_{d}(E_{i})=1/n$ and

[TABLE]

We thus have

[TABLE]

Notice that $\tilde{E}_{i}$ is a cube with sides of length $2^{-m}$. Following the proof of Theorem 2, we have

[TABLE]

where $c_{i}$ is the center of $\tilde{E}_{i}$. Therefore,

[TABLE]

The inequality (3.6) is due to $m\leq(\log_{2}n+1)/d+1$. The equality (3.7) is due to $c_{i}\in E_{i}$ and $\{E_{1},\dots,E_{n}\}$ is a split of $[0,1]^{d}$. As a result,

[TABLE]

Combing (3.8) with $\mathrm{Var}(\hat{\mu}_{n})=s_{n}^{2}/n^{2}$ establishes the inequality (3.1).

Similar to (2.5), for any $\delta>0$, there exists a constant $C_{d,\delta}$ depending on $d$ and $\delta$ such that

[TABLE]

because the diameter of $E_{i}$ is not larger than $2\sqrt{d+3}n^{-1/d}$ by Lemma 3. Using (3.8) and (3.9), the Lyapunov condition

[TABLE]

is satisfied. Finally, using the Lyapunov CLT, we obtain (3.2). ∎

From the proof of Theorem 4 in He and Owen (2016), we find that

[TABLE]

for any Lipschitz function $f$ with modulus $M$. Theorem 4 gives an asymptotic lower bound of order $n^{-1-2/d}$ for $\mathrm{Var}(\hat{\mu}_{n})$. Therefore, the rate $O(n^{-1-2/d})$ is tight for $\mathrm{Var}(\hat{\mu}_{n})$ if $f\in C^{1}([0,1]^{d})$ and $\int_{[0,1]^{d}}\|\nabla f(X)\|^{2}\mathrm{d}X>0$. To prove the asymptotic normality, we only require the lower bound as shown in the proof of Theorem 4.

3.2 Piecewise Smooth Functions

In this subsection, we focus on piecewise smooth functions of the form $f(X)=g(X)1_{\Omega}(X)$, where $\partial\Omega$ admits a $(d-1)$-dimensional Minkowski content defined below. This kind of functions was also studied in He and Owen (2016).

Definition 5**.**

For a set $\Omega\subset[0,1]^{d}$, define

[TABLE]

where $(A)_{\epsilon}:=\{x+y|x\in A,\|y\|\leq\epsilon\}$. If $\mathcal{M}(\partial\Omega)$ exists and finite, then $\partial\Omega$ is said to admit a $(d-1)$-dimensional Minkowski content.

In the terminology of geometry, $\mathcal{M}(\partial\Omega)$ is known as the surface area of the set $\Omega$. The Minkowski content has a clear intuitive basis, compared to the Hausdorff measure that provides an alternative to quantify the surface area. We should note that the Minkowski content coincides with the Hausdorff measure, up to a constant factor, in regular cases. It is known that the boundary of any convex set in $[0,1]^{d}$ has a $(d-1)$-dimensional Minkowski content since the surface area of a convex set in $[0,1]^{d}$ is bounded by the surface area of the unit cube $[0,1]^{d}$, which is $2d$. More generally, Ambrosio et al. (2008) found that $\partial\Omega$ admits a ($d-1$)-dimensional Minkowski content when $\Omega$ has a Lipschitz boundary.

Let

[TABLE]

be the indices of collections of $E_{i}$ that are interior to $\Omega$ and at the boundary of $\Omega$, respectively. Denote $|A|$ as the cardinality of the set $A$.

Lemma 6**.**

If $\partial\Omega$ admits a $(d-1)$-dimensional Minkowski content, then $|\mathcal{T}_{\mathrm{bdy}}|=O(n^{1-1/d})$.

Proof.

The proof is given in the proof of Theorem 4 in He and Owen (2016). We provide here for completeness.

From Lemma 3, the diameter of $E_{i}$, denoted by $r_{i}$, satisfies $r_{i}\leq 2\sqrt{d+3}n^{-1/d}$. Let $\epsilon=2\sqrt{d+3}n^{-1/d}$. From (3.11), for any fixed $\delta>2\mathcal{M}(\partial\Omega)$, there exists $\epsilon_{0}>0$ such that $\lambda_{d}((\partial\Omega)_{\epsilon})<\delta\epsilon$ whenever $\epsilon<\epsilon_{0}$. Assume that $n>(2\sqrt{d+3}/\epsilon_{0})^{d}$. Thus $r_{i}\leq\epsilon<\epsilon_{0}$. Note that $\cup_{i\in\mathcal{T}_{\mathrm{bdy}}}E_{i}\subset(\partial\Omega)_{\epsilon}$. This leads to

[TABLE]

which completes the proof. ∎

Theorem 7**.**

Let $f(X)=g(X)1_{\Omega}(X)$, where $g(X)\in C^{1}([0,1]^{d})$, $\Omega\subset[0,1]^{d}$ and $\partial\Omega$ admits a $(d-1)$-dimensional Minkowski content. Suppose that $\sigma_{\Omega}^{2}=\int_{\Omega}\|\nabla g(X)\|^{2}\mathrm{d}X>0$. Then for all sufficiently large $n$,

[TABLE]

If $d>2$,

[TABLE]

in distribution as $n=b^{m}\to\infty$.

Proof.

Following the notations in the proof of Theorem 4, we have

[TABLE]

Similar to the proof of Theorem 4 for $g\in C^{1}([0,1]^{d})$, we have

[TABLE]

where $c_{i}\in E_{i}$ as defined there for $i\in\mathcal{T}_{\mathrm{int}}$. From (3.14), we find that

[TABLE]

Note that

[TABLE]

where we picked $c_{i}\in E_{i}\backslash\Omega\neq\emptyset$ for $i\in\mathcal{T}_{\mathrm{bdy}}$ so that the last term of (3.16) is actually zero. Therefore, similar to (3.8), we have

[TABLE]

that establishes (3.12).

On the other hand,

[TABLE]

Again, for any $\delta>0$,

[TABLE]

where $C_{d,\delta}>0$ is some constant that only depends on $d$ and $\delta$. This leads to

[TABLE]

It follows from (3.17) that

[TABLE]

By the continuity of $g$, there is a constant $D$ with $|g(X)|\leq D$ for all $X\in[0,1]^{d}$. Therefore,

[TABLE]

By Lemma 6, we have $|\mathcal{T}_{\mathrm{bdy}}|=O(n^{1-1/d})$. As a result,

[TABLE]

provided that $d>2$ and $\delta>1$. Together with (3.19) and (3.22), the Lyapunov condition is thus verified. So the asymptotic normality is satisfied by applying the Lyapunov CLT again. ∎

He and Owen (2016) gave an upper bound of $O(n^{-1-1/d})$ for $\mathrm{Var}(\hat{\mu}_{n})$ if $f(X)=g(X)1_{\Omega}(X)$, where $g$ is Lipschitz continuous. Theorem 7 provides a lower bound of order $n^{-1-2/d}$. For discontinuous integrands, we cannot get asymptotically matching lower bound to the upper bound because when we take $\Omega=[0,1]^{d}$, the lower bound (3.12) is in line with the smooth case. To establish the asymptotic normality, Theorem 7 requires $d>2$. It is not clear in general whether the asymptotic normality holds for $d=1,2$. If the last term of (3.14) has a lower bound of $O(n^{1-1/d})$, one would have the asymptotic normality for $d=2$. For $d=1$, let’s consider the function $f(X)=g(X)1_{\{X>\theta\}}(X)$ for some $\theta\in[0,1]$. If $\theta$ is a multiple of $b^{-m_{0}}$ for some $m_{0}>0$, then the error over the set $\mathcal{T}_{\mathrm{bdy}}$ vanishes whenever $m\geq m_{0}$. The Lyapunov condition is thus verified by (3.19). As a result, the asymptotic normality holds for this case. If $\theta$ does not have a terminating $b$-adic representation, we may require some additional conditions to ensure the asymptotic normality. Theorem 7 also requires that $\int_{\Omega}\|\nabla g(X)\|^{2}\mathrm{d}X>0$. That condition actually rules out the case in which $f$ is an indicator function. The analysis of indicator functions is presented in the next subsection.

3.3 Indicator Functions

We now consider indicator functions of the form $f(X)=1_{\Omega}(X)$. Recall that $\mathcal{T}_{\mathrm{bdy}}$ denotes the index of the collections of $E_{i}$ that touch the boundary of $\Omega$. In this case, the variance of the estimate reduces to

[TABLE]

where $\mathrm{Var}(1_{\Omega}(X^{(i)}))=n\lambda_{d}(E_{i}\cap\Omega)(1-n\lambda_{d}(E_{i}\cap\Omega))$. Motivated by the proof of Theorem 7, one needs to derive a suitable lower bound for $s_{n}^{2}=n^{2}\mathrm{Var}(\hat{\mu}_{n})$ to apply the Lyapunov CLT. Note that $s_{n}^{2}\leq|\mathcal{T}_{\mathrm{bdy}}|/4$. Assume that $\partial\Omega$ admits a $(d-1)$-dimensional Minkowski content; $s_{n}^{2}$ then has an upper bound of $O(n^{1-1/d})$ since $|\mathcal{T}_{\mathrm{bdy}}|=O(n^{1-1/d})$ by Lemma 6. It is easy to see that if $s_{n}^{2}\geq cn^{1-1/d}$ for some constant $c>0$, the Lyapunov condition is satisfied for any $d>1$. However, it is possible that $\mathcal{T}_{\mathrm{bdy}}=\emptyset$ for strictly increasing sample sizes $n_{k}$, $k=1,\dots,\infty$, if $\Omega$ is a cube. This leads to an identically zero variance and hence $s_{n_{k}}^{2}=0$. Therefore, to study the asymptotic normality for indicator functions, we need the following assumption on $\Omega$, instead of the Minkowski content condition.

Assumption 8**.**

For $\Omega\subset[0,1]^{d}$, there exist a constant $c>0$ and an $N_{0}\geq 1$ such that for any $n\geq N_{0}$,

[TABLE]

Moreover,

[TABLE]

Theorem 9**.**

Let $f(X)=1_{\Omega}(X)$, where $\Omega$ satisfies Assumption 8. Then

[TABLE]

in distribution as $n=b^{m}\to\infty$.

Proof.

By (3.23) and (3.24), we have $s_{n}^{2}\geq c|\mathcal{T}_{\mathrm{bdy}}|$. The Lyapunov condition

[TABLE]

is satisfied for any $\delta>0$, where we used the condition (3.25) and $c>0$. Applying the Lyapunov CLT, we obtain (3.26). ∎

Note that for $d=1$, the condition (3.25) does not hold if $\Omega$ is a union of $k$ disjoint intervals in [0,1], where $k$ is a given positive integer. This is because $|\mathcal{T}_{\mathrm{bdy}}|\leq 2k$ for all possible $n$. Actually, for such cases, the CLT does not hold since the integration error is distributed over (at most) $k+1$ possible values for any $n=b^{m}$; see also L’Ecuyer et al. (2010) for discussions on randomly-shifted lattice rules.

Define $A_{n}(c)=\{i\in\mathcal{T}_{\mathrm{bdy}}|\mathrm{Var}(1_{\Omega}(X^{(i)}))\geq c\}$. Assumption 8 can be weakened slightly to that there exist $c>0$ and $\delta>0$ such that

[TABLE]

If $\partial\Omega$ admits a $(d-1)$-dimensional Minkowski content additionally, it suffices to verify

[TABLE]

or equivalently, $|A_{n}(c)|^{-1}=o(n^{(1/d-1)/(1+\delta)})$. This is because $|\mathcal{T}_{\mathrm{bdy}}|=O(n^{1-1/d})$.

However, it may be hard to verify Assumption 8 for general $\Omega$. As an illustrative example, we next show that the assumption holds for the case $\Omega=\{X=(X_{1},X_{2})\in[0,1]^{2}|X_{1}+X_{2}\geq 1\}$. We restrict our attention to the van der Corput sequence in base $b=2$ so that $n=2^{m}$. In this case, $E_{i}$ is a square with sides of length $1/\sqrt{n}$ when $m$ is even; when $m$ is odd, $E_{i}$ is a rectangle with width $\sqrt{2/n}$ and height $1/\sqrt{2n}$; see Figure 2 for illustrations. We thus have

[TABLE]

Moreover, for all $i\in\mathcal{T}_{\mathrm{bdy}}$, we find that

[TABLE]

Therefore, Assumption 8 is satisfied with $c=3/16$ and $N_{0}=1$ so that the CLT holds for this example. Similarly, it is easy to see that the CLT still holds for the set $\Omega=\{X=(X_{1},\dots,X_{d})\in[0,1]^{d}|\sum_{i=1}^{d}X_{i}\geq d/2\}$ in $d$ dimensions.

4 Numerical Results

In this section, we present some numerical studies to assess the normality of the standardized errors. We also examine the lower bound established in Theorem 4 for smooth functions. We consider the integrals of the following functions:

• •

a smooth function, $f_{1}(X)=12^{d/2}\prod_{i=1}^{d}(X_{i}-\frac{1}{2})$,

• •

a piecewise smooth function, $f_{2}(X)=(X_{1}-X_{2})1_{\{\sum_{i=1}^{d}X_{i}\geq d/2\}}(X)$, and

• •

an indicator function, $f_{3}(X)=1_{\{\sum_{i=1}^{d}X_{i}\geq d/2\}}(X)$.

Note that for any $d\geq 1$, the exact values of these integrals are $\mu=0$, $\mu=0$, and $\mu=1/2$, respectively. The smooth function was studied in Owen (1997), which satisfies the smooth condition required in Loh (2003). The scrambled $(t,m,d)$-net integration of this smooth function has a variance of $O(n^{-3}(\log n)^{d-1})$ Owen (1997), and it enjoys the asymptotic normality when $t=0$, as confirmed by Loh (2003). However, for the last two discontinuous functions, there is no theoretical guarantee in supporting the asymptotic normality for scrambled net quadratures, since the functions do not fit into the class of smooth functions required in Loh (2003).

We make comparisons with randomized Sobol’ points, which use the nested uniform scrambling of Owen (1995) or the linear scrambling of Matoušek (1998). We use the C++ library of T. Kollig and A. Keller (http://www.uni-kl.de/AG-Heinrich/SamplePack.html) to generate the nested uniform scrambled Sobol’ (NUS–Sobol’) points. To generate the linear scrambled Sobol’ (LS–Sobol’) points, we make use of the generator scramble in MATLAB. To calculate Hilbert’s mapping function $H(x)$, we use the C++ source code in Lawder (2000) which is based on the algorithm in Butz (1971). To estimate the variances of these estimators, we use $R$ independent replications $\hat{\mu}_{n}^{(1)},\dots,\hat{\mu}_{n}^{(R)}$ of the sampling schemes. We then estimate the variances by the corresponding empirical variances

[TABLE]

where $\bar{\mu}=(1/R)\sum_{i=1}^{R}\hat{\mu}_{n}^{(i)}$. To see asymptotic normality, we plot the kernel smoothed density of the standardized errors

[TABLE]

using the function ksdensity in MATLAB. In our experiments, we take $R=1000$ and $n=2^{14}=16384$, in order to get good accuracy in the estimation of the target density.

Now consider the smooth function $f_{1}(X)$. We find that $f_{1}\in C^{1}([0,1]^{d})$, and

[TABLE]

Therefore, by Theorem 4, the HSFC-based estimate follows the CLT for $f_{1}$. The lower bound in (3.1) becomes

[TABLE]

Note that $f_{1}(X)$ is a Lipschitz function whose modulus $M$ satisfies

[TABLE]

Together with (3.10), we obtain an upper bound

[TABLE]

Figure 3 shows the natural logarithm of the empirical variances of the HSFC-based estimator for $n=2^{m}$, $m=0,\dots,18$. The true variance of Monte Carlo sampling is $1/n$ for all $d\geq 1$. The lower bound and the upper bound above are also presented. We observe that the empirical variances decay at the rate $n^{-2}$ for $d=2$, and at the rate $n^{-5/4}$ for $d=8$. This supports that the rate $n^{-1-2/d}$ for the HSFC sampling is tight for smooth functions. Figure 4 displays smoothed density estimations of the standardized errors for plain Monte Carlo, LS–Sobol’, NUS–Sobol’, and HSFC. As expected, a nearly normal distribution appears for both the Monte Carlo and HSFC schemes. For the nested uniform scrambling scheme, a nearly normal distribution is also observed for $d=2$. This is because Sobol’ sequence is a $(t,d)$-sequence in base $b=2$ with $t=0$ for $d=2$ and $t>0$ for $d=8$ Dick and Niederreiter (2008). Therefore, the CLT holds for $d=2$, as confirmed by Loh (2003). For $d=8$, on the other hand, the density of the standardized errors does not look like a normal distribution. Even worse, for the linear scrambling scheme, the distribution of the standardized errors is very different from the normal distribution. It looks rather spiky for $d=2$.

For the two discontinuous functions $f_{2}(X)$ and $f_{3}(X)$, the CLT holds for the HSFC sampling (see Sections 3.2 and 3.3 for details). Figures 5 and 6 show smoothed density estimations of the standardized errors for the two functions, respectively. As expected, a nearly normal distribution appears for both the Monte Carlo and HSFC schemes with $d=2,8$. More interestingly, a nearly normal distribution is also observed for the nested uniform scrambling scheme in all cases, although it is not clear whether the CLT holds for scrambled net integrations of discontinuous functions. Similar to the case of the smooth function, the integration error distribution for the linear scrambling scheme is far from the normal distribution, particularly for $d=2$. Comparing to the nested uniform scrambling, the linear scrambling requires less randomness, and therefore its samples may be strongly dependent. That might explain why the CLT does not hold in most cases for randomized QMC with the linear scrambling.

5 Concluding Remarks

Loh (2003) showed that the scrambled net estimate has an asymptotic normal distribution for certain smooth functions. In a very recent work, Basu and Mukherjee (2016) found that the scrambled geometric net estimate has an asymptotic normal distribution for certain smooth functions defined on products of suitable subsets of $\Re^{d}$. The smoothness conditions required in the two papers are more restrictive than the smooth condition required in Section 3.1. The proofs in both Loh (2003) and Basu and Mukherjee (2016) relied on ensuring a suitable lower bound on the variance of the estimate matching up to constants to the upper bound. The proofs in this paper relies on establishing a suitable lower bound, and then make use of the Lyapunov CLT.

We also proved the asymptotic normality of the HSFC-based stratified estimate for certain discontinuous functions. To our best knowledge, it is not clear whether the asymptotic normality of the scrambled net estimate holds for discontinuous functions. He and Wang (2015) provided some upper bounds of scrambled net variances for piecewise smooth functions of the same form $f(X)=g(X)1_{\Omega}(X)$ studied in Section 3.2, but $g$ is of bounded variation in the sense of Hardy and Krause instead. For future research, following the procedures in Loh (2003), it is desirable to establish a matching lower bound for the variance of scrambled net integration of discontinuous functions.

He and Owen (2016) used randomized van der Corput sequence in base $b$ as the input of the HSFC sampling. This makes the sampling scheme extensible. As in Loh (2003) and Basu and Mukherjee (2016), the analysis in this paper is based on the sample size with the pattern $n=b^{m}$, not with arbitrary $n$. This scheme turns out to be a kind of stratified samplings. In contrast to the usual grid sampling, it is extensible and does not require so highly composite sample sizes, particularly for large $d$. The results on the HSFC sampling can also be applied to the usual grid sampling.

Acknowledgments

The authors gratefully thank Professor Art B. Owen for the helpful comments. Zhijian He is supported by the National Science Foundation of China under Grant 71601189. Lingjiong Zhu is supported by the NSF Grant DMS-1613164.