Digital Nets and Sequences for Quasi-Monte Carlo Methods

The most important property of quasirandom sequences is an equidistribution property. The quality of such uniformity for quasirandom sequences is measured by the discrepancy, which is a distance between the continuous uniform distribution on $[0,1)^{s}$ and the empirical distribution of the $\mathbf{y}_{i}$ for $i=1,\cdots,N$. First we define the empirical distribution of the sequence as

$$F_{N}(\mathbf{y})=\frac{1}{N}\sum_{i=1}^{N}1_{[\mathbf{y}_{i},\infty)}(\mathbf{y}),$$

(1.1)

where $1_{A}$ is the indicator function of the set $A$. We define the uniform distribution on $[0,1)^{s}$ as

$$F_{\text{unif}}(\mathbf{y})=\prod_{j=1}^{s}y_{j},\mbox{ }\mathbf{y}=(y_{1},\cdots,y_{s})^{T}\in[0,1)^{s}.$$

(1.2)

The discrepancy arises from the worst case error analysis of quasi-Monte Carlo quadratures. The Koskma-Hlawka inequality, which provides an important theoretical justification for the quasi-Monte Carlo quadratures, uses the discrepancy to provide an upper bound quadrature error of quasi-Monte Carlo methods:

$$|\hat{I}_{N}-I|\leq D^{*}(P)V_{HK}(f)$$

(1.3)

where $V_{HK}$ is said to have bounded variation $f$ in the sense of Hardy and Krause. If $V_{HK}(f)<\infty$, then $f\in BVHK$. For further understanding of $V_{HK}$ see Niederreiter [Nie92]. $D^{*}(P)$ indicates the star discrepancy which will be defined next.

The star discrepancy is the most widely studied discrepancy and it is defined as follows:

$$D^{*}(P)=\sup_{\mathbf{y}\in[0,1)^{s}}|F_{\text{unif}}(\mathbf{y})-F_{N}(\mathbf{y})|=||F_{\text{unif}}-F_{N}||_{\infty}.$$

(1.4)

The star discrepancy can be thought as a special case $(p=\infty)$ of $L_{p}$-star discrepancy which is defined by

$$D^{*}_{p}(P)=||F_{\text{unif}}-F_{N}||_{p},1\leq p\leq\infty.$$

(1.5)

Another discrepancy, which will be used in this thesis, given by [Hic96], called a generalized $L_{2}$-discrepancy, is

	$\displaystyle\hskip-19.91684ptD^{2}(P)=$
	$\displaystyle\hskip-11.38092pt-1+\frac{1}{N^{2}}\sum_{i,j=0}^{N-1}\prod_{r=1}^{s}\left[-\frac{(-\gamma^{2})^{\alpha}}{(2\alpha)!}B_{2\alpha}(\{y_{ir}-y_{jr}\})+\sum_{k=0}^{\alpha}\frac{\gamma^{2k}}{(k!)^{2}}B_{i}(y_{ir})B_{i}(y_{jr})\right],$		(1.6)

where the notation $\{y\}$ means the fractional part of a number. The positive integer $\alpha$ indicates the degree of smoothness of integrands in the underlying Hilbert space, and the parameter $\gamma$ measures the relative importance given to the uniformity of the points in low dimensional projections of the unit cube versus high dimensional projections of the unit cube. The reproducing kernel leading to this discrepancy is for a Hilbert space of integrands whose partial derivatives up to order $\alpha$ in each coordinate direction are square integrable [HH99]. It is known that the root mean square discrepancy for scrambled $(t,m,s)$-nets decays as $O(N^{-3/2+\epsilon})$ as $N\to\infty$ [HH99, HY01] for $\alpha\geq 2$. $B_{i}(y)$ is the $i$th Bernoulli polynomial (see [AS64]). The first few Bernoulli polynomials, which are used in this thesis, are $B_{0}(y)=1$ and

	$\displaystyle B_{1}(y)$	$\displaystyle=$	$\displaystyle y-\frac{1}{2},$
	$\displaystyle B_{2}(y)$	$\displaystyle=$	$\displaystyle y^{2}-y+\frac{1}{6},$
	$\displaystyle B_{3}(y)$	$\displaystyle=$	$\displaystyle y^{3}-\frac{3}{2}y^{2}+\frac{1}{2}y,$
	$\displaystyle B_{4}(y)$	$\displaystyle=$	$\displaystyle y^{4}-2y^{3}+y^{2}-\frac{1}{30}.$

See [Hic96] for further discussion about the $L_{p}$-star discrepancy and several other discrepancies.

The following two sections introduce two important families of quasi-Monte Carlo methods which are :

• •

  integration lattices [Nie92, Slo85], and

• •

  digital nets and sequences [Nie92].

Rank-1 lattices, also known as good lattice point (glp) sets, were introduced by Korobov[Kor59] and have been widely studied since then. See [SJ94, Nie92] for further details. The formula for the node points of a rank-1 lattice is simply

$$P=\{\{i\mathbf{h}/N\}:i=0,\cdots,N-1\},$$

(1.7)

where $N$ is the number of points, $\mathbf{h}$ is an $s$-dimensional generating vector of integers The formula for a lattice is rather simple. However finding a good generating vector $\mathbf{h}$ that makes the lattice have low discrepancy for the chosen $N$ and $s$ is not trivial. Recently the formula for lattice is extended to an infinite sequence [HHLL01]. This is done by using the radical inverse function, $Q_{b}(i)$. For any integer $b\geq 2$, let any non-negative integer $i$ be represented in base $b$ as $i=\cdots i_{3}i_{2}i_{1}$, where the digits $i_{k}$ take on values between $0$ and $b-1$. The function $Q_{b}(i)$ flips the digits about the decimal point, i.e.,

$$Q_{b}(i)=0.i_{1}i_{2}i_{3}\cdots(\mbox{ base }b)=\sum_{k=1}^{\infty}i_{k}b^{-k}.$$

(1.8)

The sequence $\{Q_{b}(i):i=0,1,\cdots\}$ is called the van der Corput sequence. An infinite sequence of imbedded lattices is defined by replacing $i/N$ by $Q_{b}(i)$, i.e.,

$$y_{i}=\{\{Q_{b}(i)\mathbf{h}\}:i=0,1,\cdots\},$$

(1.9)

where the first $N$ points of this sequence are a lattice whenever $N$ is a power of $b$ [HHLL01].

Digital nets and sequences [Lar98b, Nie92] are one widely used types of low discrepancy point sets.

Let $b\geq 2$ denote a prime number. For any non-negative integer $i=\cdots i_{3}i_{2}i_{1}(\mbox{base }b)$, we define the $\infty\times 1$ vector $\boldsymbol{\psi}(i)$ as the vector of its digits, i.e., $\boldsymbol{\psi}(i)=(i_{1},i_{2},\ldots)^{T}$. For any point $z=0.z_{1}z_{2}\cdots(\mbox{base }b)\in[0,1)$, let $\boldsymbol{\phi}(z)=(z_{1},z_{2},\ldots)^{T}$ denote the $\infty\times 1$ vector of the digits of $z$. Let ${\mathbf{C}}_{1},\ldots,{\mathbf{C}}_{s}$ denote predetermined $\infty\times\infty$ generator matrices. The digital sequence in base $b$ is $\{\mathbf{y}_{0},\mathbf{y}_{1},\mathbf{y}_{2},\ldots\}$, where each $\mathbf{y}_{i}=(y_{i1},\ldots,y_{is})^{T}\in[0,1)^{s}$ is defined by

$$\boldsymbol{\phi}(y_{ij})=\mathbf{C}_{j}\boldsymbol{\psi}(i),\quad j=1,\ldots,s,\ i=0,1,\ldots.$$

(1.10)

Here all arithmetic operations take place in mod $b$. Thus, the right side of (1.10) should not give a vector ending in an infinite trail of $b-1$s.

Digital nets and sequences are the special cases of $(t,m,s)$-nets and $(t,s)$-sequences and further explanation is given in Chapter 4.

The concept of $t$ in $(t,m,s)$-nets and $(t,s)$-sequences provides one way to measure the quality of low discrepancy sequences. See [Nie92] for more detail explanations. Here we provide a brief outline of elementary intervals of $(t,m,s)$-nets and $(t,s)$-sequences.

Let $s\geq 1$ and $b\geq 2$ be integers. An elementary interval in base $b$ is a subinterval of $[0,1)^{s}$ of the form

$$\prod_{r=1}^{s}\left[\frac{a_{r}}{b^{k_{r}}},\frac{a_{r}+1}{b^{k_{r}}}\right)$$

with integer $k_{r}\geq 0$, $0\leq a_{r}<b^{k_{r}}$. Such an elementary interval has volume $b^{-(k_{1}+\cdots+k_{s})}$.

Let $m\geq 0$ be an integer. A finite set of $N=b^{m}$ points from $[0,1)^{s}$ is a $(t,m,s)$-net in base $b$ if every elementary interval in base $b$ with the volume $b^{t-m}$ contains exactly $b^{t}$ points of a set, where $t$ is a nonnegative integer. Smaller $t$ values imply a better equidistribution property for the net. Obviously the best case is when $t=0$, that is every $b$-ary box of volume $1/N$ contains exactly one of the $N$ points of the set.

For a given $t\geq 0$, an infinite sequence of points from $[0,1)^{s}$ is a $(t,s)$ sequence in base $b$ if for all integers $k\geq 0$ and $m\geq t$ the finite sequence $\{y_{kb^{m}+1},\cdots,y_{(k+1)b^{m}}\}$ is a $(t,m,s)$-net in base $b$.

The constructions of $(t,m,s)$-nets and $(t,s)$-sequences in base $b=2$ were introduced by Sobol’ [Sob67]. Later Faure [Fau82] provided $(0,m,s)$-nets and $(0,s)$-sequences for prime bases $b\geq s$. The general definitions are given by Niederreiter [Nie87]. Mullen, Mahalanabis and Niederreiter [Nie95] provide tables of attainable $t$-values for nets.

The advantage of using nets taken from $(t,s)$-sequences is that one can increase $N$ through a sequence of values $N=\lambda b^{m}$, $1\leq\lambda<b$, and all the points used in $\hat{I}_{\lambda b^{m}}$ are also used in $\hat{I}_{(\lambda+1)b^{m}}$. Owen [Owe97a] introduces the $(\lambda,t,m,s)$-nets to describe such sequences.

The initial $\lambda b^{m}$ points of a $(t,s)$-sequence are well distributed. For integers $m,t,b,\lambda$ with $m\geq 0$, $0\geq t\leq m$, and $1\leq\lambda<b$, a sequence ${\mathbf{y}_{i}}$ of $\lambda b^{m}$ points is called a $(\lambda,t,m,s)$-net in base $b$ if every $b$-ary box of volume $b^{t-m}$ contains $\lambda b^{t}$ points of the sequence and no $b$-ary box of volume $b^{t-m-1}$ contains more than $b^{t}$ points of the sequence.

For functions of bounded variation in the sense of Hardy and Krause, the numerical integration by averaging over the points of a $(t,m,s)$-net has an error of order $N^{-1}(\log N)^{s-1}$. Niederreiter [Nie92] discusses more precise error bounds for $(t,m,s)$-nets and $(t,s)$-sequences.

Quasi-Monte Carlo methods can obtain a better convergence rate than Monte Carlo methods, because the points are chosen to be more uniform. However deterministic quasi-Monte Carlo methods have disadvantages compared with Monte Carlo methods. First, quasi-Monte Carlo methods are statistically biased since the points are chosen in certain deterministic way. Therefore the mean of the quadrature estimate is not the desired integral. Second, quasi-Monte Carlo methods do not facilitate simple error estimates while Monte Carlo methods provide straightforward probabilistic error estimates.

Randomizing quasi-Monte Carlo point is one way to overcome such disadvantages while still preserving their higher convergence rates. One may think of randomized quasi-Monte Carlo as a sophisticated variance reduction technique. There are two types of randomization. One is adding the same $s$-dimensional random shift to every point. Let $P$ be the original low discrepancy point set, then $P_{sh}=\{i\mathbf{h}/N+\Delta\}:i=0,\cdots,N-1\}$, where $\Delta$ is a random vector uniformly distributed on $[0,1)^{s}$. This type of randomization, which was introduced in [CP76], is often used with lattice rule because shifted lattice rules retain their lattice structure. However, shifted nets do not necessarily remain as nets. Another type of randomization, which was proposed by Art Owen [Owe95], is a more sophisticated method that randomly permutes the digits of each point. And it often applies to nets because scrambled nets retain their nets properties. However, scrambled lattice rules do not necessary remain lattice rules. The more detail discussions on Owen’s scrambling can be found in Chapter 2.

In this chapter we have introduced the necessary background of quasi-Monte Carlo methods.

Chapter 2 discusses the randomization of the quasi-Monte Carlo methods specifically randomizing digital nets. Detailed descriptions of Owen’s randomization (or scrambling) and Faure-Tezuka’s randomization (or tumbling) are provided. The actual implementation of randomization for digital $(t,s)$-sequences and some numerical results are presented.

Chapter 3 explores the distribution of the discrepancy of scrambled digital $(t,m,s)$-nets which is based on recent theory. We fit the empirical distribution by a sum of chi squared random variables. The distribution of the discrepancy of randomized lattice points is presented also.

Finding better digital nets is the main content of Chapter 4. Here we use optimization methods to find the generator matrices for digital $(t,m,s)$-nets. Evolutionary computation is introduced as a tool for a finding better net.

Chapter 5 presents the applications of the sequences that we have generated. Various problems are explored and the performance of the new sequences is compared with other existing low discrepancy sequences.

Finally the thesis ends with some concluding remarks.

Owen [Owe95, Owe97a, Owe97b, Owe98, Owe99] proposes scrambled $(t,s)$-sequences as a hybrid of Monte Carlo and quasi-Monte Carlo methods. This clever randomizations of quasi-Monte Carlo methods can combine the best of both by yielding higher accuracy with practical error estimates. The following is a detailed description of Owen’s scrambling.

Let $\{\mathbf{x}_{0},\mathbf{x}_{1},\mathbf{x}_{2},\ldots\}$ denote the randomly scrambled sequence proposed by Owen. Let $x_{ijk}$ denote the $k^{\text{th}}$ digit of the $j^{\text{th}}$ component of $\mathbf{x}_{i}$, and similarly for $y_{ijk}$. Then

	$\displaystyle x_{ij1}$	$\displaystyle=$	$\displaystyle\pi_{j}(y_{ij1}),$
	$\displaystyle x_{ij2}$	$\displaystyle=$	$\displaystyle\pi_{y_{ij1}}(y_{ij2}),$
	$\displaystyle x_{ij3}$	$\displaystyle=$	$\displaystyle\pi_{y_{ij1},y_{ij2}}(y_{ij3}),\quad\ldots,$
	$\displaystyle x_{ijk}$	$\displaystyle=$	$\displaystyle\pi_{y_{ij1},y_{ij2},\cdots,y_{ijk-1}}(y_{ijk}),\ldots$

where the $\pi_{a_{1}a_{2}\dots}$ are random permutations of the elements in mod $b$ chosen uniformly and mutually independently. Owen [Owe98] provides a geometrical description to help visualize his scrambling as follows: The rule for choosing $x_{ij}$ is like cutting the unit cube into $b$ equal parts along the $\mathbf{x}_{j}$ axis and then reassembling these parts in random order to reform the cube. The rule for choosing $x_{ij2}$ is like cutting the unit cube into $b^{2}$ equal parts along the $\mathbf{x}_{j}$ axis, taking them as $b$ groups of $b$ consecutive parts, and reassembling the $b$ parts within each groups in random order. The rule for $x_{ijk}$ involves cutting the cube into $b^{k}$ equal parts along the $\mathbf{x}_{j}$ axis, forming $b^{k-1}$ groups of $b$ equal parts, and reassembling the $b$ parts within each group in random order.

Owen [Owe95, Owe97b] states that a randomized net inherits certain equidistribution properties of $(t,m,s)$-nets by proving the following propositions:

Proposition 1  If $\{\mathbf{y}_{i}\}$ is a $(\lambda,t,m,s)$-net in base $b$ then $\{\mathbf{x}_{i}\}$ is a $(\lambda,t,m,s)$-net in base $b$ with probability 1.

Proposition 2  Let $\mathbf{y}$ be a point in $[0,1)^{s}$ and let $\mathbf{x}$ be the scrambled version of $\mathbf{y}$ as described above. Then $\mathbf{x}$ has the uniform distribution on $[0,1)^{s}$.

Faure and Tezuka [FT00] proposed another type of randomizing digital nets and sequences. However the effect of the Faure-Tezuka-randomization can be thought as re-ordering the original sequence, rather than permuting its digits like the Owen-scrambling. Thus we refer to this method as Faure-Tezuka-tumbling as suggested by Art Owen (personal communication). The following describes details of Faure-Tezuka-tumbling.

For any $m,\lambda=0,1,\ldots$, let $i=\lambda b^{m},\ldots,(\lambda+1)b^{m}-1$. Then $\boldsymbol{\psi}(i)$ takes all possible values in its top $m$ rows. By the same token $\mathbf{L}^{T}\boldsymbol{\psi}(i)+\mathbf{e}$ takes on all possible values in its top $m$ rows, but not necessarily in the same order. Therefore, the Faure-Tezuka-tumbled $(t,m,s)$-net, $\{\mathbf{z}_{\lambda b^{m}},\ldots,\mathbf{z}_{(\lambda+1)b^{m}-1}\}$, obtained by replacing $\boldsymbol{\psi}(i)$ of (1.10) by $\mathbf{L}^{T}\boldsymbol{\psi}(i)+\mathbf{e}$ is just the same set as the original $(t,m,s)$-net, $\{\mathbf{y}_{\nu b^{m}},\ldots,\mathbf{y}_{(\nu+1)b^{m}-1}\}$, given by (1.10) for some $\nu$.

There is a cost in scrambling $(t,s)$-sequences. The manipulation of each digit of each point requires more time than for an unscrambled sequence. In addition Owen’s scrambling can be rather tedious because of the bookkeeping required to store the many permutations. Here we present an alternative method that is only slightly less general than Owen’s proposal but is an efficient method of scrambling that minimizes this cost. This is called an Owen-scrambling to recognize that it is done in the spirit of Owen’s original proposal.

The discrepancy used in here is the squared discrepancy (1.1) and it has been scaled by a dimension-dependent constant so that the root mean square scaled discrepancy of a simple random sample is $N^{-1/2}$.

Figures 2.2 shows the histogram of the distribution of the square discrepancy of the particular Owen-scrambled Sobol’ sequence for the different numbers of scrambled digits, $K$. $K$ has been chosen to be $10$, $20$, and $30$ for $s=2$ and $m=10$ with 200 independent replications. The range of the $x$-axis is are $10^{-15},\ldots,10^{-6}$ in this figure. From the figure, there is only little distinction in the distribution of the square discrepancy for $K=20$ and $30$. However a larger number of scrambled digits produces a wider range of discrepancies. From Figure 2.3 the root mean square discrepancies for $K=20$ and $30$ are nearly the same. However the root mean square discrepancy for $K=10$ looses the benefit from scrambling as $m$ increases.

Figures 2.3 and 2.4 plot the discrepancy of the particular Owen-scrambled Sobol’ sequence for different choices of $K$ and the choices are same as before. The choices of dimension are $s=2$ and $5$. These figures show the root mean square discrepancy of 100 different replications. There is almost no difference in the root mean square discrepancy between scrambling $K=20$ or $30$ digits. However for $K=10$ the discrepancy looses its superior performance as $N$ increases, which indicates an insufficient number of scrambled digits. Also notice that choosing the number of scrambled digits is independent from $s$.

Figure 2.5 plots the discrepancies of the Faure sequence, the scrambled Faure sequence, and the tumbled Faure sequence. The number of scrambled digits, $K$, is chosen to be 20, and 100 different replications have been performed. The choice of dimension is $s=5$. Since the Faure sequence takes a base $b$ as a smallest prime number which is equal to or greater than $s$, $b$ assigned to be $5$. Therefore $N$ is chosen to be $\lambda 5^{m}$, where $\lambda=1,\cdots,b-1$ and $m=1,\cdots,4$. From Figure 2.5 scrambled sequences perform the best among all sequences and the tumbled sequence performs better than the original Faure sequence. Also notice that the scrambling and tumbling procedure help to flatten out the humps, which are present in the original sequence when the value of $\lambda\neq 1$.

Figures 2.6–2.8 show the root mean square discrepancies of randomly scrambled $(t,m,s)$-nets in base 2 with $N=2^{m}$ points that have been calculated by using 100 different replications. The choices of dimension are $s=1$ and $10$. The number of scrambled digits, $K$, is chosen to be 31.

Theory predicts that the root mean square discrepancy of nets should decay as $\operatorname{O}(N^{-1})$ for $\alpha=1$, and as $\operatorname{O}(N^{-3/2})$ for $\alpha=2$. Figures 2.6–2.8 display this behavior, although for larger $s$ the asymptotic decay rate is only attained for large enough $N$. However, even for smaller $N$ the discrepancy decays no worse than $\operatorname{O}(N^{-1/2})$, the decay rate for a simple Monte Carlo sample.

The discrepancy of unscrambled nets does not attain the $\operatorname{O}(N^{-3/2})$ decay rate for $\alpha=2$. Although for these experiments all possible digits were scrambled, our experiment suggests that scrambling about $K=m+10$ digits is enough to obtain the benefits of scrambling. Scrambling more than this number does not improve the discrepancy but increases the time required to generate the scrambled points.

The discrepancy measures the uniformity of the distribution of a set of points and it can be interpreted as the maximum possible quadrature error over a unit ball of integrands. Since the discrepancy depends only on the quadrature rule, it is often used to compare different quadrature rules.

Let $K(x,y)$ be the reproducing kernel for some Hilbert space of integrands whose domain is $[0,1)^{s}$. To approximate the integral

$$\int_{[0,1)^{s}}f(x)\ dx$$

one may use the quasi-Monte Carlo quadrature rule

$$\frac{1}{N}\sum_{z\in P}f(z),$$

where $P\subset[0,1)^{s}$ is some well-chosen set of $N$ points. The error of this rule is

$$\operatorname{Err}(f)=\int_{[0,1)^{s}}f(x)\ dx-\frac{1}{N}\sum_{z\in P}f(z),$$

and the square discrepancy is

$$D^{2}(P)=\operatorname{Err}_{x}(\operatorname{Err}_{y}(K(x,y))),$$

where $\operatorname{Err}_{x}$ implies that the error functional is applied to the argument $x$ [Hic99].

The reproducing kernel may be written as the (usually infinite) sum

$$K(x,y)=\sum_{\nu}\phi_{\nu}(x)\phi_{\nu}(y),$$

where $\{\phi_{\nu}\}$ is a basis for the Hilbert space of integrands [Wah90]. This implies that the discrepancy may be written as the sum [HW01]

$$D^{2}(P)=\sum_{\nu}[\operatorname{Err}(\phi_{\nu},P)]^{2}.$$

By Loh’s central limit theorem it is known that each $\operatorname{Err}(\phi_{\nu},P)$ is asymptotically distributed as a normal random variable with mean zero. However, the $\operatorname{Err}(\phi_{\nu},P)$ are in general correlated.

By making an orthogonal transformation one may write

$$D^{2}(P)\approx\sum_{\nu=1}^{\infty}\beta_{\nu}X_{\nu},$$

(3.1)

where the $X_{\nu}$ are independent $\chi^{2}(1)$ random variables, and the $\beta_{1}\geq\beta_{2}\cdots$ are some constants. Note that if $\beta_{1}=\cdots=\beta_{n}=\beta$ and $0=\beta_{n+1}=\beta_{n+2}=\cdots$ then $D^{2}(P)/\beta$ is approximately distributed as $\chi^{2}(n)$.

Since equation (3.1) involves an infinite sum it must be further approximated in order to be computationally feasible. This is done by approximating all but the first $n$ terms by a constant:

$$D^{2}(P)\approx\sum_{\nu=1}^{n}\beta_{\nu}X_{\nu}+c$$

(3.2)

where

$$c=E[\sum_{\nu=n+1}^{\infty}\beta_{\nu}X_{\nu}].$$

Here $E$ denotes the expectation.

Let $Z$ denote the random variable given by the square discrepancy of a scrambled net, and let $W$ denote the random variable given by the right hand side of (3.2). Let $G_{Z}$ and $G_{W}$ denote the probability distribution functions of these two random variables, respectively. Ideally, one would find the $\beta_{\nu}$ and $c$ that minimize some loss function involving $G_{Z}$ and $G_{W}$, but these distributions are not known exactly. Thus, the optimal $\beta_{\nu}$ and $c$ are found by minimizing

$$\sum_{i=1}^{M}[\log(\hat{G}_{Z}^{-1}(p_{i}))-\log(\hat{G}_{W}^{-1}(p_{i})]^{2}.$$

(3.3)

Here,

$$p_{i}=(2i-1)/(2M),\text{ for }i=1,\ldots,M,$$

and

$$\hat{G}_{Z}^{-1}(p_{i})=Z_{(i)},$$

the $i^{\text{th}}$ order statistic from a sample of $M$ scrambled net square discrepancies. Also,

$$\hat{G}_{W}^{-1}(p_{i})=W_{(k_{i})},$$

the $k_{i}^{\text{th}}$ order statistic from a sample of $L$ independent and identically distributed drawings of $W$, where

$$k_{i}=iL/M-(L-M)/(2M),$$

and $L$ is an odd multiple of $M$. The logarithm is used in (3.3) because for a fixed $N$ and $s$ the square discrepancy values can vary by a factor of 100 or more, and it is not good for the larger values to unduly influence the fitted distribution.

Fitting the distribution of the square discrepancy by minimizing (3.3) relies on several approximations. The central limit theorem is invoked to obtain (3.1). The infinite sum is replaced by a finite one in (3.2). The probability distributions of each side of (3.2) are approximated by Monte Carlo sampling, and the two probability distributions are compared at only a finite number of points. In spite of these approximations the fitted distribution matches the observed distribution of the square discrepancy rather well.

The discrepancy considered here is (1.1). With $\alpha=2$ and $\gamma=1$ the discrepancy can be rewritten as

$$D^{2}(P)=-1+\frac{1}{N^{2}}\sum_{\mathbf{y},\mathbf{y}^{\prime}\in P}^{N}\prod_{j=1}^{s}\left\{1+\gamma\left[B_{1}(\mathbf{y}_{j})B_{1}(\mathbf{y}^{\prime}_{j})+\frac{1}{4}B_{2}(\mathbf{y}_{j})B_{2}(\mathbf{y}^{\prime}_{j})\right.\right.\\ \left.\left.-\frac{1}{24}B_{4}(\text{\boldmath$\{$}\mathbf{y}_{j}-\mathbf{y}^{\prime}_{j}\text{\boldmath$\}$})\right]\right\}$$

(3.4)

Figures 3.1 and 3.3 show the empirical distribution (dashed line) of the square discrepancy of scrambled nets and its fitted value (thin line). The fit is based on $M=1000$ replications of the scrambled Sobol’ [BF88, HH01], Niederreiter-Xing, and Faure nets, $L=11000$, and $n=3$ or $4$. The optimization was done by using MATLAB’s fminsearch function. The fits are good in the middle of the distributions, but off in the tails. Moreover the Q-Q plots of the empirical and fitted distributions bear this out.

There are several humps in the empirical distribution of the square discrepancy for $s=2$. Moreover as $N$ increases, the number of humps increases. However, these humps are less noticeable as the dimension increases.

The humps can be explained as follows: The scrambling technique used here, as described in Chapter 2, scrambles the generator matrices of these digital nets and also gives them a digital shift. Since this technique preserves the digital nature of the nets, it can be shown that for $l=0,1,\ldots$ the sum of the digits of points numbered $lb^{2},\ldots,(l+1)b^{2}-1$, is zero modulo $b$. However, for Owen’s original scrambling this is not necessarily the case. For example, in one dimension Owen’s scrambling is equivalent to Latin hypercube sampling. Figures 3.9 and 3.10 show the empirical distribution of the square discrepancy for a scrambled one-dimensional Sobol’ net and for Latin hypercube sampling. Since the base of the Sobol’ sequence is $2$ the difference in the two graphs emerges for $N\geq b^{2}=4$. Note that although the distributions of the square discrepancy for Owen’s original scrambling and the variant used here are somewhat different, the means of the two distributions are the same.

The square discrepancy can be written as a sum of polynomials in $\gamma$:

	$\displaystyle D^{2}(P,K)=$		$\displaystyle(\gamma-\gamma_{0})\tilde{D}^{2}_{1}(P,K)+(\gamma-\gamma_{0})(\gamma-\gamma_{1})\tilde{D}^{2}_{2}(P,K)$
			$\displaystyle+\cdots+(\gamma-\gamma_{0})\cdots(\gamma-\gamma_{s-1})\tilde{D}^{2}_{s}(P,K),$

Then the sum of polynomials can be rearranged as

$$D^{2}(P,K)=\gamma D^{2}_{1}(P,K)+\gamma^{2}D^{2}_{2}(P,K)+\cdots+\gamma^{s}D^{2}_{s}(P,K),$$

where $D^{2}_{j}$ measures the uniformity of all $j$-dimensional projections of $P$ [Hic98]. A computationally efficient method for obtaining all of the $D^{2}_{j}$ is to compute $D^{2}(P,K)$ for $s$ different values of $\gamma_{0},\cdots\gamma_{s-1}$ and then perform polynomial interpolation. Specifically Newton’s divided difference formula has been applied for polynomial interpolation.