Generation of space-filling uniform designs in unit hypercubes

Journal of Statistical Planning and Inference 142 (2012) 3189–3200

Contents lists available at SciVerse ScienceDirect

Journal of Statistical Planning and Inference journal homepage: www.elsevier.com/locate/jspi

Generation of space-filling uniform designs in unit hypercubes Ismael S. Talke n, John J. Borkowski Department of Mathematical Sciences, Montana State University, Bozeman, MT 59717, USA

a r t i c l e i n f o

abstract

Article history: Received 17 December 2011 Received in revised form 30 March 2012 Accepted 17 June 2012 Available online 26 June 2012

In this paper a new discrepancy measure of uniformity for uniform designs (UDs) in a unit cube is presented. Alternative measures of uniformity based on distance criteria which can be applied to higher dimensions are also discussed. The good lattice point (glp) method was used to construct the uniform designs. Two approaches (generator equivalence and projection) of reducing the computational cost of the glp method are proposed and discussed. Two examples are presented in this paper. & 2012 Elsevier B.V. All rights reserved.

Keywords: Discrepancy Distance-based criteria Good lattice points Space-filling design Uniform design

1. Introduction Experimental design is a statistical methodology used for studying the relationship between the experimental factors on one or more responses of interest. Experiments often have many factors requiring methodology that can identify those factors having large effects. If there are enough resources to run all combinations of the factor levels, then a full factorial design is often used. However, the availability of enough resources to run a full factorial design does not always exist. In such cases fractional factorial designs are often used to significantly reduce the number of experimental runs. We might, however, encounter an experimental situation when a fractional factorial design can no longer reduce the number of runs to an allowable number due to resource limitations. In this situation a design that provides useful experimental data while minimizing the number of runs has to be used. One such type of experimental design is a Uniform Design (UD). Since their introduction, UDs have been widely used in different fields and have proved to be successful for examining the relationship between the experimental (input) factors and the response (output) variables when the underlying model form is unknown (Fang and Lin, 2003; Fang et al., 2000). When using a classical experimental design (such as central composite or Box–Behnken designs), it is assumed that the underlying model form is known but with unknown model parameter values. A design is then chosen with the goal of producing parameter estimates having high precision (Fang and Lin, 2003). In comparison, the UDs are not model dependent. Thus, a UD is robust to potential models because it can still prove useful even when the exact form of the model is not known (Liang et al., 2001; Li et al., 2004; Zhang et al., 1998). A UD also allows the largest possible number of levels for each factor among all experimental designs for a fixed number of design points. Unlike the classical experimental designs that are defined in terms of combinatorial structures, UDs are assessed on measures of the spread of the design points over the experimental domain. Thus, UDs are one type of space-filling designs.

n

Corresponding author. E-mail addresses: [email protected], [email protected] (I.S. Talke), [email protected] (J.J. Borkowski).

0378-3758/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jspi.2012.06.013

3190

I.S. Talke, J.J. Borkowski / Journal of Statistical Planning and Inference 142 (2012) 3189–3200

A good space-filling design is one in which the design points are uniformly scattered throughout the experimental region (Borkowski and Piepel, 2009). Fang (1980) and Wang and Fang (1981) used number-theoretic methods (NTM’s) (also known as quasi-Monte Carlo methods) to generate uniformly scattered sets of points throughout the experimental domain. In this paper, we consider designs having s factors with the unit cube C s ¼ ½0; 1s as the experimental domain. The goal is to find a ‘representative set of n points’ on Cs. To be representative, the set of n design points P ¼ fx1 , . . . ,xn g in Cs should be ‘scattered evenly’ in Cs. Hence, the aim is to find P that minimizes a measure that assesses the uniformity of scatter. The measures considered are discrepancy and distance-based criteria. Extensive reviews of UDs are given in Wang and Fang (1990), Fang et al. (1999, 2006), Fang and Lin (2003), Fang and Wang (1994), and Liang et al. (2001). In Section 2, we present a brief discussion related to several discrepancy measures. A new discrepancy criterion, called modified star discrepancy (MSTRD), is introduced and its benefits for use in low dimensions are justified. Moreover, distance-based measures of uniformity are presented as alternatives to MSTRD in higher dimensions. In Section 3, the good lattice point (glp) method of constructing UDs is presented, and two approaches for reducing the computational cost of the glp method are introduced. Section 4 contains examples and discussion, and concluding remarks are in Section 5. 2. Measures of uniformity As stated in the Introduction, the reason to use a uniform design is to provide design points that are uniformly scattered throughout the experimental region C s ¼ ½0; 1s : In this section, existing measures of uniformity will be reviewed followed by a new measure. Alternative distance-based measures of uniformity will also be reviewed. Thus, the focus is on quantifiable measures of uniformity that allow for direct comparison of competing designs. Different measures of uniformity have been adopted by different authors. For example, Hickernell (1998a,b) modified the star L2 discrepancy measure and obtained the star L2 centered discrepancy and the star L2 wrap-around discrepancy measures of uniformity. Hickernell also provided analytic formulas for computing both these discrepancy measures. Fang and Wang (1994) also used the star discrepancy in measuring uniformity while Borkowski and Piepel (2009) used distance-based measures of uniformity for highly constrained mixture designs. 2.1. The NTM discrepancy measure Let FðxÞ be the multivariate uniform distribution on the unit cube experimental region C s ¼ ½0; 1s , and let F n ðxÞ be the empirical distribution of the design points P ¼ fx1 , . . . ,xn g F n ðxÞ ¼

n 1X Ifx r x1 ,xi2 r x2 , . . . ,xis r xs g, n i ¼ 1 i1

where x ¼ ðx1 , . . . ,xs Þ, xi ¼ ðxi1 , . . . ,xis Þ for i¼1,y,n, and IðAÞ ¼ 1 if A is true, and is zero otherwise. Thus, F n ðxÞ equals the proportion of points in P that lie in the hyperrectangle formed by the origin and point x. The star Lp-discrepancy is defined as Z 1=p p 9F n ðxÞFðxÞ9 dx : ð1Þ Dp ðPÞ ¼ Cs

The star discrepancy is obtained by setting p ¼ 1 in the star Lp-discrepancy Z 1=p p 9F n ðxÞFðxÞ9 dx ¼ sup9F n ðxÞFðxÞ9: D1 ðPÞ ¼ lim p-1

Cs

xEC s

The star L2 centered discrepancy and the star L2 wrap-around discrepancy measures described by Hickernell (1998a,b) are special cases of the star Lp-discrepancy given in (1). Hickernell provided analytic formulas for these two discrepancies that are easy to compute. However, because they are modifications, both can yield biased underestimates of Dp ðP n Þ. The application of discrepancy measures to number theoretic methods for generating designs will now be discussed. For any vector g ¼ ðg1 , g2 , . . . , gs Þ 2 C s , let Nðg,PÞ ¼ nF n ðgÞ ¼ the number of design points that satisfy Ifxi1 r g1 ,xi2 r g2 , . . . ,xis r gs g ¼ 1. Let the hyperrectangle formed by g and the origin be denoted ½0, g. The star discrepancy measure for a set of design points P is given by   Nðg,PÞ  Vð½0, gÞ, Dðn,PÞ ¼ sup ð2Þ n g2C s Q where Vð½0, gÞ ¼ si ¼ 1 gi denotes the volume of ½0, g. In other words, for Dðn,PÞ, we are comparing the proportion of points Nðg,PÞ=n that lie in ½0, g to the ratio of the volumes Vð½0, gÞ=VðC s Þ. Note that Vð½0, gÞ=VðC s Þ ¼ Vð½0, gÞ because VðC s Þ ¼ 1. Then, Dðn,PÞ is simply the supremum of the absolute difference over all possible hyperrectangles ½0, g for all g 2 C s . If Dðn,PÞ is small then the proportion of points within each hyperrectangle ½0, g is nearly proportional to the volume of the hyperrectangle indicating good uniformity of scatter. Or, in other terms, it is a good space-filling design. Conversely, a large discrepancy measure indicates a poor space-filling design. Dðn,PÞ in (2) is equivalent to the Kolomogorov Smirnov statistic in goodness-offit testing. For an overview of discrepancy measures, see Fang and Wang (1994) and Fang et al. (2000, 2006).

I.S. Talke, J.J. Borkowski / Journal of Statistical Planning and Inference 142 (2012) 3189–3200

1

1

0.9

0.9

0.8

0.8

0.7

3191

0.7 γ

0.6

γ

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0 0

0.5

1

0

0.5

1

Fig. 1. Geometric illustration of regions used to calculate (a) the star discrepancy criteria and (b) modified star discrepancy for a 34 point design.

Fig. 1(a) geometrically illustrates how the star discrepancy Dðn,PÞ measures uniformity. Consider the two (inner and outer) rectangles displayed in Fig. 1(a). There are 34 design points in the outer rectangle corresponding to C 2 ¼ ½0; 12 and 11 design points in ½0, g in the inner rectangle of width 0.5 and height 0.6 (i.e., g ¼ ð0:5,0:6ÞÞ. Thus, Nðg,PÞ=n ¼ 11=34 ¼ 0:3235 ¼ the proportion of points that fall in the rectangle ½0, g. Then 90:3235Vð½0, gÞ9 ¼ 90:32350:30009 ¼ 0:0235 is the discrepancy for this particular g. Therefore, to find the star discrepancy (STRD) Dðn,PÞ, we must consider these absolute deviations over all g 2 C 2 , and then take the maximum. In this paper, however, instead of working with only the hyperrectangles ½0, g for g in Cs, we propose to partition the experimental region Cs into 2s hyperrectangles formed by g and each of the 2s vertices of Cs, and then take the maximum of the 2s discrepancies associated with this set of 2s hyperrectangles. Notationally, let Rj be the jth hyperrectangle formed by g and vertex vj (j ¼ 1, . . . ,2s ) of Cs. Let Nðg,P,Rj Þ be the number of design points in Rj and VðRj Þ be the volume of Rj. Then define   Nðg,P,Rj Þ  VðRj Þ: Dðn,P, gÞ ¼ max s  n j ¼ 1,...,2 The proposed modified star discrepancy (MSTRD) is now defined as MSTRDðn,PÞ ¼ supDðn,P, gÞ:

ð3Þ

g2C s

For example, for s¼ 2, C2 is divided into four rectangles for each g, and is illustrated in Fig. 1(b). The discrepancy Dð34,P, gÞ values for the bottom left (j¼1), bottom right (j¼2), top left (j¼3), and top right (j¼4) rectangles are, respectively, 9 11 34 0:5 9 6 8 0:69 ¼0.0235, 9 34 0:5  0:69¼ 0.0353, 9 34 0:4  0:59¼0.0235, and 9 34 0:4  0:59 ¼0.0353. We then take the maximum¼ 0.0353, which is the numerical measure of modified star discrepancy for this g. Therefore, to find the modified star discrepancy MSTRDðn,PÞ, we must consider taking the maximum over all g 2 C 2 . Note that the MSTRD is greater than or equal to the STRD. As shown in this example STRD¼0.0235 is smaller than MSTRD¼0.0353 (or  67% smaller). Thus, using STRD to compare the uniformity of scatted for competing designs can potentially lead to the selection of an inferior design. The MSTRD also provides an improvement to the STRD as a discrepancy measure because it assesses uniformity of scatter for more subspaces of Cs. The concept of dividing the experiential region into 2s hyperrectangles appears, at first glance, to be similar to the approach used by Hickernell in the star L2 centered discrepancy. His approach, however, is restricted to using only the center point ð0:5, . . . ,0:5Þ (i.e., a single g) to form the hyperrectangles. In this study, we remove this single center point restriction and permit the use of random evaluation points to form the hyperrectangles. The MSTRD approach, like that of Hickernell’s method, is invariant under rotation. Determining exact or near-exact MSTRD (or STRD) values can be computationally expensive because we are optimizing over all g 2 C s . Thus, as the dimension of the experimental region Cs increases, the computational cost for calculating a measure of discrepancy increases. In this paper, the estimation of MSTRD and STRD values is restricted to low (s¼2,3) dimensional unit cubes. For higher dimensions, the distance-based measures in Borkowski and Piepel (2009) are recommended alternatives to discrepancy-based measures. 2.2. The distance-based approach Johnson et al. (1990) introduced minimax and maximin designs which are space-filling designs with points dependent on distance measures (Koehler and Owen, 1996). Recently, Borkowski and Piepel (2009) proposed and used an alternative version of the maximin distance criterion as well as two other distance-based criteria for highly constrained mixture

3192

I.S. Talke, J.J. Borkowski / Journal of Statistical Planning and Inference 142 (2012) 3189–3200

designs for selecting the best space-filling designs. In addition to the proposed MSTRD, the three distance-based criteria proposed by Borkowski and Piepel (2009) will be used in this study. The details of the distance-based approach are now presented. Let x1 , . . . ,xn be any set of n points in Cs and which form the rows of the n  s design matrix X 2 3 2 3 x11 x12 : : : x1s x1 6x 7 6x 7 x : : : x 22 2s 7 6 21 6 27 6 7 6 7 6 7 : : : : : 7 X¼6 6 : 7¼6 : 7 6 : 7 6 : 7 : : : : : 4 5 4 5 xn1 xn2 : : : xns xn The distance between a vector x in Cs and matrix X is defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dðx,XÞ ¼ minJxxj J2 ¼ min 2 ðxxj Þðxxj Þ0 j

j

for j ¼ 1; 2, . . . ,n. That is, dðx,XÞ is the distance between a point x in Cs and the nearest design point in X. The following distance-based criteria are the measures of uniformity for assessing the space-filing properties that were used by Borkowski and Piepel (2009): qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1. The root mean squared distance RMSD ðXÞ ¼ E½dðx,XÞ2  is the square root of the expected value of the squared distance. 2. The average distance ADðXÞ ¼ E½dðx,XÞ is the expected value of the distance. 3. The maximum distance MDðXÞ ¼ maxx2C s ðdðx,XÞÞ is the maximum distance between any point x in Cs and the design point of X closest to x. The ADðXÞ and RMSDðXÞ criteria are based on expectations, and thus require numerical integration over the experimental region to determine their exact values, while the MDðXÞ criterion is the solution to a complex optimization problem. Because determination of exact criterion values is computationally demanding for all three criteria, Borkowski and Piepel (2009) approximated each criterion using a large random sample of points in their mixture design regions. For design region Cs, the corresponding approximations are sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PN 2 2 k ¼ 1 ðdðuk ,XÞÞ rmsdðXÞ ¼ N PN adðXÞ ¼

k¼1

dðuk ,XÞ N

mdðXÞ ¼ max ðdðuk ,XÞÞ k ¼ 1,...,N

where u1 , . . . ,uN is a random sample of points in Cs. A good space-filling design will have small rmsdðXÞ, adðXÞ and mdðXÞ values because small values indicate the points in Cs are close to the design X (i.e., there are no large volume subregions in Cs that do not contain any design points from X). The same applies for the MSTRD. Conversely, large values of these criteria indicate that a design is poor with respect to its space-filling properties. 3. The good lattice point (glp) method A number of authors have proposed different methods of constructing UDs. For each method, we will assume the number of design points n is given and the unit cube C s ¼ ½0; 1s is the experimental region. See Fang and Wang (1994) and Hua and Wang (1981) for a review of different construction methods. The first method considered is the good lattice point (glp) method which is a very efficient quasi-Monte Carlo method proposed by Korobov (1959) and later adopted by many authors such as Fang and Wang (1994) and Hua and Wang (1981). Let glpðn,sÞ denote an n-point glp design having s factors and let gcdða,bÞ be the greatest common divisor of integers a and b. The first step is to find the candidate set of positive integers (called generators). That is, find the ordered set of positive integers Hn ¼ fh1 ,h2 , . . . ,hk g such that 1 r hi ohj on for io j, and gcdðn,hi Þ ¼ 1 for i¼1,y,k. k is the maximum number of generators and is determined by Euler’s f function:      1 1 1 1 . . . 1 k ¼ fðnÞ ¼ n 1 p1 p2 pr where n ¼ pt11 pt22    ptrr is the prime number decomposition of n. Next, for each vector ðh01 ,h02 , . . . ,h0s Þ of s distinct elements of Hn, generate an n  s matrix U ¼ ðuij Þ where uij ¼ ih0j modðnÞ for i¼1,y,n and j¼1,y,s where 1r uij rn. Note that each column of U is a permutation of ð1; 2, . . . ,nÞ because the gcdðn,hj Þ ¼ 1 for each j. Then, from U an n  s dimensional design

I.S. Talke, J.J. Borkowski / Journal of Statistical Planning and Inference 142 (2012) 3189–3200

3193

matrix X is formed. The ith row of X ¼ ðxi1 ,xi2 , . . . ,xis Þ where xij ¼ ð2uij 1Þ=2n for i¼1,y,n and j¼ 1,y,s. Each row of X represents a design point in Cs and each column of X is a permutation of ð1=2n,3=2n, . . . ,ð2n1Þ=2nÞ: See Fang and Wang (1994), Zhang et al. (1998), Fang and Lin (2003), and Fang et al. (2006) for details. For example, for n ¼15 points and s ¼2, consider generating vector h ¼ ð1; 2Þ. This vector generates matrix U and glpð15; 2Þ design X given by 2 2 3 3 1 2 0:0333 0:1000 6 2 6 7 7 4 7 6 6 0:1000 0:2333 7 7 6 6 7 6 3 6 0:1667 0:3667 7 6 7 7 6 6 7 7 6 4 6 7 8 7 6 6 0:2333 0:5000 7 7 6 6 7 6 5 10 7 6 0:3000 0:6333 7 7 6 6 7 6 6 12 7 6 0:3667 0:7667 7 7 6 6 7 7 6 6 7 6 7 14 7 6 0:4333 0:9000 7 7 6 6 7 7 6 6 7 1 7 and X ¼ 6 0:5000 0:0333 7: U¼6 8 7 6 6 7 6 9 6 0:5667 0:1667 7 3 7 7 6 6 7 7 6 6 7 6 10 5 7 6 0:6333 0:3000 7 7 6 6 7 6 11 7 7 6 0:7000 0:4333 7 7 6 6 7 7 6 6 7 6 12 9 7 6 0:7667 0:5667 7 7 7 6 6 6 13 11 7 6 0:8333 0:7000 7 7 7 6 6 7 7 6 6 4 14 13 5 4 0:9000 0:8333 5 15

15

0:9667

0:9667

The complete set of X matrices is generated by considering each of ðfðnÞ s Þ possible generating vectors. A measure of uniformity (e.g., discrepancy) is then calculated for each of the ðfðnÞ s Þ matrices. The matrix Xðn,hnÞ with generating vector n n n h ¼ ðh1 , . . . ,hs Þ that minimizes the uniformity measure is then selected as the best UD. For example, for n ¼15 the only prime factors are p1 ¼ 3 and p2 ¼ 5. Thus, fð15Þ ¼ 8 with H15 ¼ f1; 2,4; 7,8; 11,13; 14g. For s¼2 there are ð82Þ ¼ 28 candidate n glpð15; 2Þ designs. Therefore, the aim is to choose generating vector hn such that design Xð15,h Þ has the smallest discrepancy among the 28 possible designs. However, Fang and Wang (1994) showed that fixing h1 ¼ 1 reduces the 7 number of candidate designs to consider. Thus, only ðfðnÞ1 s1 Þ ¼ ð1Þ ¼ 7 candidate designs need to be considered instead of 28 designs. When the number of design points n is a prime p and the dimensionality of the experimental region s are both large (Fang et al., 2006; Fang and Wang, 1994) described the power generator method using primitive roots to significantly reduce the extremely large number of candidate generating vectors. By limiting the search to designs generated by this method, generating vectors that produce the best UDs might be excluded. In addition to the glp method, three other methods of constructing UDs were initially considered. These were the cyclotomic field good point method, the Halton method, and the Hammersley method (Fang and Wang, 1994; Halton and Smith, 1964; Hammersley, 1964; Halton, 1960). Our preliminary comparisons among these four UD construction methods for design sizes n ¼ 10; 11, . . . ,50 clearly indicated that the glp method generated the best designs in Cs for s ¼ 2; 3,4. Thus, this research is focused only on the glp method for design generation. Because of the computational demands of calculating STRD and MSTRD, using the glp method to find the best design can also inflate the computational cost due to a large number ðfðnÞ1 s1 Þ of possible generating vectors. To circumvent the computational requirements when using the glp method, two approaches are proposed and discussed in the next section. 3.1. Computational reduction method I: equivalence Because the MSTRD and the distance-based measures are invariant under rotation (of 901,1801,2701, . . .), and because some generating vectors produce points that are simply a rotation of points from other generating vectors, there is no need to consider sets of generating vectors that yield the same measure of uniformity under rotation. Two NT-nets generated by the glp method are defined to be equivalent if the points in one NT-net can be generated by a column permutation of the other. Therefore, because of the invariance under rotations, two equivalent NT-nets generated by the glp method will have equal values for the uniformity measures considered in this paper. Two computational reduction methods for the glp method for s¼2 and s¼3 will now be presented. For s ¼ 2, consider generating vectors ð1,hi Þ and ðhj ,1Þ where hj 4 hi . If modðhi  hj ,nÞ ¼ 1, then the set of NT-net points generated from ð1,hi Þ is equivalent to the set of NT-Net points generated from ð1,hj Þ. For s ¼ 2, the permutation matrix,  0 1 P2 ¼ 1 0 is used to permute the columns of ðhj ,1Þ into ð1,hj Þ. That is, the set of points in the NT-net generated by ðhj ,1Þ  P2 is the same as those generated by ð1,hi Þ. Hence, only measures of uniformity for NT-nets generated by ð1,hi Þ (and not ð1,hj Þ) need to be computed.

3194

I.S. Talke, J.J. Borkowski / Journal of Statistical Planning and Inference 142 (2012) 3189–3200

For s ¼ 3, consider generating vectors ð1,hi ,hj Þ, ðhk ,1,hl Þ and ðhm ,hq ,1Þ where hj 4 hi and hk 4hi . If the following six conditions: (i) (ii) (iii) (iv) (v) (vi)

modðhi  hk ,nÞ ¼ 1, modðhj  hm ,nÞ ¼ 1, modðhi  hl  hm ,nÞ ¼ 1, modðhl  hq ,nÞ ¼ 1, modðhj  hk  hq ,nÞ ¼ 1, modðhi  hj  hk  hl  hm  hq ,nÞ ¼ 1

are satisfied, then, with respect to a uniformity measure, the set of NT-net points generated from ð1,hi ,hj Þ is equivalent to the set of NT-net points generated by permuting the columns of ðhk ,1,hl Þ and ðhm ,hq ,1Þ. To see this, note the following: 1. If hk 4hl , then the NT-net formed from ð1,hi ,hj Þ is equivalent to the NT-net formed from ðhk ,1,hl Þ. The permutation matrix 2 3 0 0 1 6 7 P 3A ¼ 4 1 0 0 5 0 1 0 is used to permute the columns of ðhk ,1,hl Þ into ð1,hi ,hj Þ. That is, the set of points in the NT-net generated by ðhk ,1,hl Þ  P 3A is the same as those generated by ð1,hi ,hj Þ. Hence, only measures of uniformity for NT-nets generated by ð1,hi ,hj Þ (and not ðhk ,1,hl Þ) need to be computed. If, however, hl 4hk , then interchange the 1st and 3rd rows of the permutation matrix P3A. 2. Similarly, if hm 4hq , then the NT-net formed from ð1,hi ,hj Þ is equivalent to the NT-net formed from ðhm ,hq ,1Þ. The permutation matrix 2 3 0 0 1 6 7 P 3B ¼ 4 0 1 0 5 1 0 0 is used to permute the columns of ðhm ,hq ,1Þ into ð1,hm ,hq Þ. That is, the set of points in the NT-net generated by ðhm ,hq ,1Þ  P 3B is the same as those generated by ð1,hi ,hj Þ. Hence, only measures of uniformity for NT-nets generated by ð1,hi ,hj Þ (and not ðhm ,hq ,1Þ) need to be computed. If, however, hq 4 hm , then interchange the 1st and 2nd rows of the permutation matrix P3B. Hence, using equivalence can significantly reduce the computational demands of finding the best UD using the glp method. For, s¼2 and s¼3, this can lead to a maximum reduction of 1/2 and 2/3, respectively. For example, consider n ¼34 point UDs. Table 1 summarizes the sets of equivalent generating vectors for s¼2 and s¼3. For example, for s¼2, (1,3) and (23,1) are equivalent, and thus, the NT-nets generated from ð1; 3Þ and (1,23) have equal measures of uniformity. Thus, these measure only need to be computed for the NT-net generated by (1,3). Similarly, for s ¼3, (1,3,5), (23,1,13), and (7,21,1) are equivalent, and thus, after using the appropriate permutations of (23,1,13) and (7,21,1), the NT-nets generated from (1,3,5), (1,13,23) and (1,7,21) all have equal measures of uniformity. Thus, these measures only need to be computed for the NT-net generated by (1,3,5). After a check of equivalence, the generators that are in the first column of Table 1(a) and (b) form the reduced candidate set of generators for n ¼34. For s ¼ 2, the computations for n ¼34 are reduced from 15 to 8 generators (  47% reduction), and, for s ¼ 3, the computations are reduced by 23 (  67%) from 105 to 35 generators. For dimension s ¼4, 24 conditions are required to determine sets of equivalent generators, with the number of conditions rapidly increasing as s increases. Even with computational reduction using equivalence, the number of generating vectors to be considered will become very large as the dimension increases and especially when the number of design points is prime. For s ¼4, a Matlab software program was written to search for the equivalent NT-nets. For higher dimensions, an alternative approach to using equivalence to reduce computations is considered and is now presented. 3.2. Computational reduction method II: projection The goal is to provide an alternative in using generator equivalence to reduce the computational burden of finding very good, if not the best, UDs generated by the glp method in a unit cube Cs, especially for higher dimensions. The proposed method is called the projection method and is defined by the following procedure: 1. Find the best glp generating vectors ð1,hi Þ for C2 or, when available, find the best ð1,hi ,hj Þ for C3. These vectors are easy to find because of the minimal computation required for low dimensions. 2. Project the best UDs generated in Step 1 to the next higher dimension. This is accomplished by creating a new set of

I.S. Talke, J.J. Borkowski / Journal of Statistical Planning and Inference 142 (2012) 3189–3200

3195

Table 1 Equivalent generators of NT-nets for s ¼2 and s ¼3 when n¼ 34. (a) Equivalent NT-nets for s ¼2 and n¼ 34 (1,3) (1,5) (1,9) (1,hi) (23,1) (7,1) (19,1) (hj,1) (b) Equivalent NT-nets for s ¼ 3 and n ¼34

(1,11) (31,1)

(1,13) (21,1)

(1,15) (25,1)

(1,27) (29,1)

(1,33) –

ð1,hi ,hj Þ

ðhk ,1,hl Þ

ðhm ,hq ,1Þ

ð1,hi ,hj Þ

ðhk ,1,hl Þ

ðhm ,hq ,1Þ

(1,3,5) (1,3,7) (1,3,9) (1,3,11) (1,3,13) (1,3,15) (1,3,19) (1,3,21) (1,3,25) (1,3,27) (1,3,29) (1,3,31) (1,3,33) (1,5,7) (1,5,9) (1,5,11) (1,5,19) (1,5,21)

(23,1,13) (23,1,25) (23,1,3) (23,1,15) (23,1,27) (23,1,5) (23,1,29) (23,1,7) (23,1,31) (23,1,9) (23,1,21) (23,1,33) (23,1,11) (7,1,15) (7,1,29) (7,1,9) (7,1,31) (7,1,11)

(7,21,1) (5,15,1) (19,23,1) (31,25,1) (21,29,1) (25,7,1) (9,27,1) (13,5,1) (15,11,1) (29,19,1) (27,13,1) (11,33,1) (33,31,1) (5,25,1) (19,27,1) (31,19,1) (9,11,1) (13,31,1)

(1,5,27) (1,5,29) (1,5,31) (1,5,33) (1,9,13) (1,9,15) (1,9,21) (1,9,25) (1,9,31) (1,9,33) (1,11,13) (1,11,25) (1,11,27) (1,11,29) (1,13,21) (1,13,25) (1,15,27)

(7,1,19) (7,1,33) (7,1,13) (7,1,27) (19,1,9) (19,1,13) (19,1,25) (19,1,33) (19,1,11) (19,1,15) (31,1,29) (31,1,27) (31,1,21) (31,1,15) (21,1,33) (21,1,15) (25,1,29)

(29,9,1) (27,33,1) (11,21,1) (33,29,1) (21,19,1) (25,21,1) (13,15,1) (15,33,1) (11,31,1) (33,25,1) (21,27,1) (15,29,1) (29,13,1) (27,25,1) (13,33,1) (15,25,1) (29,27,1)

glp (11;1,7)

glp (11; 1,7, 9)

1.0 0.8 1.0

0.6 y

0.8 0.4 z

0.6 0.4

0.2

y

0.2

1.0 0.8 0.6 0.4 0.2 0.0

0.0 0.0

0.2

0.4

0.6 x

0.8

1.0

0.0 0.0 0.2 0.4 0.6 0.8 1.0 x

Fig. 2. Geometrical illustration of projection from s ¼2 to s ¼3 dimensions.

generating vectors by adding an additional h value to the existing best generating vectors while retaining the sequence of increasing values in the generating vectors. Compute the measure of uniformity for each of these vectors of length s þ 1. Retain the vectors that generated UDs having the smallest measure of uniformity. 3. Iterate the procedure in Step 2 until the desired dimension is reached. 4. Retain the UD having the smallest measure of uniformity. For design size n¼11, Fig. 2 shows the geometric projection of the best glp generator from C2 into C3. The best design for s¼ 2 had generator ð1; 7Þ. Adding 9 to ð1; 7Þ to create generator ð1; 7,9Þ was the optimal choice. This figure illustrates why this approach works well. That is, for any design having good space-filling properties in 3 dimensions, it must also have good space-filling properties in two-dimensional projections. Given designs of size n projection from dimension s to s þ 1 is done by considering all the best generators found in dimension s as candidate generators to be projected to s þ 1: For example, to generate n ¼100 design points for s ¼ 10, start with s¼2. As can be seen in Table 2 the generator (1,23) is found to be the best based on the rmsd and md criteria and (1,13) based on the ad criterion and though not reported in Table 2 the generator (1,37) was the best based on MSTRD for s¼ 2. Hence, (1,23), (1,13) and (1,37) are used as candidate generators that will be projected to s ¼3. In the rows for s¼3 in Table 2 the generator (1,17,23) is best based on the rmsd and ad criteria, (1,7,37) is best based on md criterion, and though

3196

I.S. Talke, J.J. Borkowski / Journal of Statistical Planning and Inference 142 (2012) 3189–3200

Table 2 Best generators for n ¼100 point glp designs using the projection method for s ¼ 3, . . . ,10 and the distance based measures from Scrambled Sobol sequence, Sobol sequence method and the random LHD for s ¼ 2; 3, . . . ,10:. Measures of uniformity s

Generator

rmsd

Generator

ad

Generator

md

2 3 4 5 6 7 8 9 10

(1,23) (1,17,23) (1,7,37,83) (1,7,37,51,83) (1,7,37,51,83,91) (1,7,37,51,69,83,91) (1,7,11,37,51,69,83,91) (1,7,11,37,51,61,69,83,91) (1,7,11,33,37,51,61,69,83,91)

0.0425 0.1139 0.1993 0.2840 0.3698 0.4496 0.5263 0.5985 0.6675

(1,13) (1,17,23) (1,7,37,83) (1,7,37,51,83) (1,7,37,51,83,91) (1,7,37,51,69,83,91) (1,7,11,37,51,69,83,91) (1,7,11,37,51,61,69,83,91) (1,7,11,33,37,51,61,69,83,91)

0.0395 0.1087 0.1923 0.2760 0.3610 0.4406 0.5174 0.5898 0.6590

(1,23) (1,7,37) (1,17,23,37) (1,3,7,37,83) (1,3,7,37,57,83) (1,7,31,37,51,83,91) (1,7,37,51,53,69,83,91) (1,7,11,37,51,59,69,83,91) (1,7,11,37,41,51,61,69,83,91)

0.0930 0.2325 0.4056 0.5427 0.7043 0.8192 0.9091 1.0117 1.1089

s

2 3 4 5 6 7 8 9 10

rmsd

ad

md

Scrambled

Sobol

LHD

Scrambled

Sobol

LHD

Scrambled

Sobol

LHD

0.0522 0.1308 0.2167 0.3072 0.3866 0.4671 0.5428 0.6088 0.6754

0.0532 0.1330 0.2201 0.3075 0.3859 0.4650 0.5360 0.6076 0.6748

0.0547 0.1360 0.2212 0.3090 0.3858 0.4661 0.5388 0.6106 0.6745

0.0468 0.1224 0.2074 0.2967 0.3765 0.4568 0.5328 0.5997 0.6660

0.0476 0.1242 0.2099 0.2966 0.3753 0.4546 0.5263 0.5981 0.6657

0.0484 0.1262 0.2105 0.2983 0.3757 0.4560 0.5291 0.6014 0.6657

0.1259 0.3029 0.4868 0.6526 0.7611 0.9314 0.9832 1.0537 1.1709

0.1327 0.3664 0.4629 0.6412 0.7607 0.8900 0.9541 1.0582 1.1627

0.1610 0.3889 0.5345 0.6394 0.7745 0.8615 0.9931 1.0856 1.1845

not reported generator (1,23,37) was the best based on MSTRD. Thus, (1,7,37), (1,17,23) and (1,23,37) are the candidate generators to be projected to s ¼ 4, and those that appear in the s ¼4 rows were found to be the best. The projection process continues from s¼4 until the desired dimension s¼10 is reached. That is, the best generators in dimension s serve as base generators for dimension s þ1 starting for s ¼ 2; 3, . . . ,9: The bold-faced values in Table 2 are the best generators for n ¼100 design points in dimension s¼10 using the projection method. To investigate the performance of the proposed design from the glp method using the projection method in Table 2, three other designs are also considered. These are the latin hypercube design (LHD) introduced by McKay et al. (1979), the Sobol sequence of Sobol (1967) and the Scrambled Sobol sequence. The Sobol sequence and the Scrambled Sobol sequence of points are generated using the algorithms by Bratley and Fox (1988) and Hong and Hickernell (2003), respectively. The Matlab function sobolset in the Statistics toolbox was used to implement these algorithms. The lhsdesign function from the Statistics Toolbox of Matlab was used to randomly generates LHDs. Table 2 contains the computed three distance based measures for the Scrambled Sobol sequence method, the Sobol sequences and the LHD for n ¼100 and s ¼ 2; 3, . . . ,10. The proposed glp design performs best by all the three distance based measures compared to these three methods used for comparison. This demonstrates how efficient the proposed method can be. Note that the values of the measures of uniformity in Table 2 increase as the dimension increases for all methods. The equivalence approach and the projection method approach are both effective ways to reduce the computational demands for finding very good, if not the best, designs generated using the glp method. The projection method, in particular, is simple and computationally thrifty, and is, therefore, recommended for generating good designs in higher dimensions.

4. Examples and discussion In this section, numerical examples and details are presented for C2 and C3 for design size n ¼34. The best generating vectors for s¼2 and s ¼3 are presented in Tables 3 and 4, respectively. The generating vector in bold-face produces the UD having the smallest value for each measure of uniformity. For example, for s ¼2, (1,13) is the generating vector producing the 34-point UD having the smallest measures of uniformity for all measures except STRD. The number of random g points used for s ¼2 is 4000. Fig. 3 contains scatterplots of points for the eight candidate designs. One can easily see that the values for STRD are all smaller than those for MSTRD, and thus STRD can seriously underestimate true discrepancy, and, more importantly, using STRD as the selection criterion could also lead to selection of an inferior UD. For s¼3 dimensions, the same five measures are used to assess the uniformity of scatter of the design points in C3. For s ¼3, we will have 23 ¼ 8 hyperrectangles created by g and the origin. Then, the number, and subsequently the proportion, of design that fall in each of these eight hyperrectangles is determined. Each of the eight proportions is then compared to the volume of the corresponding hyperrectangle, and the maximum of the eight absolute differences between the

I.S. Talke, J.J. Borkowski / Journal of Statistical Planning and Inference 142 (2012) 3189–3200

3197

Table 3 Measures of uniformity for the n¼ 34 point glp designs. Generator

rmsd

ad

md

STRD

MSTRD

(1,3) (1,5) (1,9) (1,11) (1,13) (1,15) (1,27) (1,33)

0.1098 0.0773 0.0811 0.1025 0.0743 0.0789 0.0757 0.2479

0.0949 0.0707 0.0737 0.0900 0.0691 0.0724 0.0700 0.2065

0.3235 0.2145 0.2286 0.2890 0.1753 0.1853 0.1781 0.6667

0.1113 0.0867 0.0660 0.0841 0.0531 0.0525 0.0635 0.2467

0.1113 0.0867 0.0667 0.0936 0.0531 0.0578 0.0694 0.2467

Table 4 Measures of uniformity for the n¼ 34 point glp designs in C3. Generator

rmsd

ad

md

STRD

MSTRD

(1,3,5) (1,3,7) (1,3,9) (1,3,11) (1,3,13) (1,3,15) (1,3,19) (1,3,21) (1,3,25) (1,3,27) (1,3,29) (1,3,31) (1,3,33) (1,5,7) (1,5,9) (1,5,11) (1,5,19) (1,5,21) (1,5,27) (1,5,29) (1,5,31) (1,5,33) (1,9,13) (1,9,15) (1,9,21) (1,9,25) (1,9,31) (1,9,33) (1,11,13) (1,11,25) (1,11,27) (1,11,29) (1,13,21) (1,13,25) (1,15,27)

0.1958 0.1818 0.1729 0.1730 0.1706 0.1834 0.1837 0.1719 0.1714 0.1795 0.1977 0.2956 0.3006 0.1797 0.1801 0.1845 0.1819 0.1707 0.1809 0.2880 0.1960 0.2867 0.1776 0.1769 0.1785 0.2910 0.1731 0.2860 0.1942 0.1788 0.1697 0.1806 0.2883 0.1779 0.1774

0.1830 0.1712 0.1636 0.1637 0.1617 0.1724 0.1722 0.1627 0.1618 0.1693 0.1851 0.2605 0.2675 0.1685 0.1689 0.1730 0.1713 0.1617 0.1700 0.2500 0.1836 0.2497 0.1673 0.1671 0.1688 0.2540 0.1640 0.2486 0.1822 0.1682 0.1608 0.1696 0.2501 0.1677 0.1671

0.4634 0.3922 0.3742 0.3808 0.3595 0.4044 0.4188 0.3817 0.3610 0.3781 0.4818 0.7160 0.7211 0.4485 0.4452 0.4078 0.3980 0.3873 0.4357 0.7122 0.4558 0.7183 0.3831 0.3701 0.3910 0.7098 0.3496 0.7145 0.4232 0.3895 0.3716 0.3936 0.7112 0.3824 0.3961

0.1192 0.1242 0.1360 0.1054 0.1054 0.1360 0.1054 0.1301 0.1197 0.1192 0.1261 0.2264 0.2392 0.1081 0.0811 0.0943 0.0977 0.0886 0.1049 0.2264 0.0940 0.2392 0.0927 0.0903 0.0801 0.2264 0.0973 0.2392 0.0956 0.1035 0.0987 0.1116 0.2264 0.0732 0.0951

0.1192 0.1271 0.1360 0.1283 0.1158 0.1360 0.1307 0.1390 0.1307 0.1334 0.1261 0.2759 0.2537 0.1081 0.1063 0.1429 0.1268 0.1231 0.1161 0.2629 0.1241 0.2590 0.1032 0.1054 0.0894 0.2477 0.1143 0.2459 0.1013 0.1146 .1185 0.1217 0.2464 0.0980 0.1101

proportions and their corresponding volumes is taken. This process was repeated for 10,000 random g points in C3, and the maximum of the 10,000 maximum differences is the estimated MSTRD value. The five discrepancy and distance-based measures of uniformity for the UDs in C3 with n ¼34 points are presented in Table 4. Recall that from Table 1(b), the glp method originally has 105 candidate generators but this was reduced to only 35 generators using equivalence. These are the generators appearing in Table 4. The best generators for each measure correspond to the bold-faced values in each column. Thus, for the glp method, ð1; 11,27Þ is the best generator for the rmsd and ad criteria, ð1; 9,31Þ for the md criterion, and ð1; 13,25Þ and ð1; 9,21Þ for the STRD and MSTRD discrepancy criteria, respectively. This highlights that different criteria can lead to different ‘best’ designs. The 34 design points generated by (1,11,27) and by ð1; 9,21Þ are displayed in Fig. 4. In Fig. 5, the values of the rmsd and MSTRD criteria for s¼3 and the rmsd criterion for s ¼4 obtained from the equivalence and projection methods are plotted for the best generators for design sizes n ¼ 10; 11, . . . ,50. Fig. 5 indicates that for the rmsd criterion, both methods are equally efficient for generating good UDs in 3 and 4 dimensions. However, for the MSTRD criterion, the equivalence method tends to produce UDs having smaller measures of uniformity than UDs

3198

I.S. Talke, J.J. Borkowski / Journal of Statistical Planning and Inference 142 (2012) 3189–3200

glp (34,1,3)

glp (34,1,5)

glp (34,1,9)

glp (34,1,11)

1

1

1

1

0.5

0.5

0.5

0.5

0

0 0

0.5

1

0 0

glp (34,1,13)

0.5

1

0 0

glp (34,1,15)

0.5

1

0

glp (34,1,27)

1

1

1

0.5

0.5

0.5

0.5

0

0.5

0

1

0

0.5

1

0

0

1

glp (34,1,33)

1

0

0.5

0.5

1

0

0

0.5

1

Fig. 3. Plots of the glp design points for n¼ 34.

Best by rmsd and ad (34;1 11 27)

Best by MSTRD (34;1 9 21)

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0 1

0 1 1 0.5

1 0.5

0.5 0

0

0.5 0

0 3

Fig. 4. Best 34-point glp UDs in C for the rmsd and ad criteria (left) and MSTRD criterion (right).

0.32

Projection Equivalent

0.3

Projection Equivalent

0.42 0.4

0.28

0.25

0.26 0.24 0.22 0.2

0.38 rmsd for s=4

MSTRD for s=3

rmsd for s=3

0.44

Projection Equivalent

0.3

0.2

0.36 0.34 0.32

0.15 0.3 0.28

0.18

0.1 0.26

0.16 10

20

30

40

Design Points

50

10

20

30

40

Design Points

50

10

20

30

40

50

Design Points

Fig. 5. The projection and the equivalent methods for the design sizes 10–50 based on rmsd for dimension 3 and 4 and based on the MSTRD for dimension 3 only.

I.S. Talke, J.J. Borkowski / Journal of Statistical Planning and Inference 142 (2012) 3189–3200

3199

Table 5 Best generating vectors of the glp designs for s ¼2 and s ¼3 for the rmsd and MSTRD criteria. s¼2

s¼2

n

rmsd

MSTRD

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

(1,3) (1,7) (1,5) (1,5) (1,9) (1,11) (1,5) (1,5) (1,5) (1,14) (1,13) (1,13) (1,5) (1,9) (1,19) (1,7) (1,11) (1,16) (1,11) (1,18) (1,7)

(1,3) (1,7) (1,5) (1,5) (1,9) (1,11) (1,7) (1,11) (1,7) (1,14) (1,9) (1,5) (1,13) (1,9) (1,17) (1,11) (1,11) (1,16) (1,11) (1,9) (1,19)

n 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

s¼3

s¼3

rmsd

MSTRD

n

rmsd

MSTRD

(1,12) (1,7) (1,14) (1,13) (1,29) (1,7) (1,8) (1,7) (1,7) (1,17) (1,9) (1,25) (1,25) (1,7) (1,26) (1,13) (1,18) (1,41) (1,9) (1,11)

(1,14) (1,9) (1,14) (1,13) (1,13) (1,11) (1,23) (1,27) (1,25) (1,31) (1,17) (1,19) (1,19) (1,27) (1,19) (1,27) (1,14) (1,41) (1,22) (1,21)

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

(1,3,9) (1,3,5) (1,5,11) (1,4,6) (1,3,5) (1,7,13) (1,3,7) (1,4,10) (1,5,7) (1,3,12) (1,3,13) (1,4,10) (1,5,7) (1,9,20) (1,11,19) (1,8,14) (1,3,17) (1,4,17) (1,5,13) (1,16,20) (1,7,13)

(1,3,7) (1,3,4) (1,5,7) (1,4,5) (1,3,5) (1,7,13) (1,5,9) (1,4,5) (1,5,7) (1,6,8) (1,3,17) (1,4,5) (1,5,13) (1,7,18) (1,11,17) (1,6,11) (1,5,11) (1,5,19) (1,5,13) (1,9,17) (1,7,13)

n 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

rmsd

MSTRD

(1,10,17) (1,3,13) (1,4,23) (1,11,27) (1,11,19) (1,7,17) (1,7,17) (1,7,17) (1,4,14) (1,13,17) (1,8,19) (1,25,31) (1,12,35) (1,7,27) (1,4,17) (1,5,33) (1,4,18) (1,13,19) (1,6,15) (1,7,23)

(1,6,22) (1,5,7) (1,8,14) (1,9,21) (1,8,22) (1,5,7) (1,7,20) (1,9,21) (1,11,23) (1,17,33) (1,9,26) (1,11,25) (1,10,18) (1,25,27) (1,28,37) (1,13,29) (1,15,28) (1,13,19) (1,15,29) (1,11,21)

resulting from the projection method. In general, however, the differences in the measures of uniformity between the two methods is very small. 5. Conclusion In this research, a new discrepancy measure of uniformity (MSTRD) for a unit cube was introduced. MSTRD is invariant under rotation and will produce consistent estimates of discrepancy (unlike STRD which can seriously underestimate true discrepancy). Besides the glp method of constructing UDs, the cyclotomic field, Halton, and Hammersley methods were studied but not presented here because of poor performance in generating good UDs. The glp method, however, is computationally expensive, especially for higher dimensional Cs and when the design size n is prime. Two approaches were proposed to reduce the computational requirements for using the glp method. These methods involve the exploitation of generator equivalence as well as the projection approach to produce UDs in higher dimensions. The limitation of the study is that the proposed new MSTRD discrepancy measure is only feasible for a low dimensional unit cube. For higher dimensions, the use of a distance-based measure (rmsd, ad or md) is recommended for design selection. Current research is now being done to find implementable discrepancy measures and computational reduction methods for experimental regions other than hypercubes. Table 5 contains a summary of the best generating vectors of glp designs for s¼2 and s¼3 for the rmsd and MSTRD criteria. A complete summary of the efficiency measures for s ¼2,3 and 4 and for n ¼ 10; 11, . . . ,50 and the associated generating vectors can be found at http://www.math.montana.edu/  jobo/cr/. References Borkowski, J.J., Piepel, G.F., 2009. Space-filling designs for highly-constrained mixture experiments. Journal of Quality Technology 41 (1), 35–47. Bratley, P., Fox, B.L., 1988. Algorithm 659 Implementing Sobol’s quasirandom sequence generator. ACM Transactions on Mathematical Software 14 (1), 88–100. Fang, K.T., 1980. The uniform design: application of number theoretic methods in experimental design. Acta Mathematicae Applicatae Sinica 3, 363–372. Fang, K.T., Li, R., Sudjianto, A., 2006. Design and Modeling for Computer Experiments. Chapman and Hall, CRC, Boca Raton. Fang, K.T., Lin, D.K.J., 2003. Uniform experimental designs and its application in industry. In: Rao, C.R., Khattree, R. (Eds.), Handbook of Statistics in Industry, vol. 22; 2003. Fang, K.T., Lin, D.K.J., Winker, P., Zhang, Y., 2000. Uniform design: theory and application. Technometrics 42, 237–248. Fang, K.T., Shiu, W., Pan, J., 1999. Uniform design based on latin squares. Statistica Sinica 9, 905–912. Fang, K.T., Wang, Y., 1994. Number-Theoretic Methods in Statistics. Chapman and Hall, London. Halton, J.H., 1960. On efficiency of certain quasi-random sequence of points in evaluating multidimensional integrals. Numerische Mathematik (2), 84– 90. Halton, J.H., Smith, G.B., 1964. Algorithm 247, radical-inverse quasi-random points sequences. Communications of the ACM (7), 701–702. Hammersley, J.M., 1964. Monte-Carlo methods for solving multivariable problems. Annals of the New York Academy of Science 86, 844–874. Hickernell, F.J., 1998a. A generalized discrepancy and quadrature error bound. Mathematics of Computation 67 (221), 299–332. Hickernell, F.J., 1998b. Lattice rules: How well do they measure up?. In: Hellekalek, P., Larcher, G. (Eds.), Random and Quasi-Random Point Sets. Lecture Notes in Statistics, vol. 138; 1998, pp. 109–166. Hong, H.S., Hickernell, F.J., 2003. Algorithm 823: implementing scrambled digital sequences. ACM Transactions on Mathematical Software 29 (2), 95–109.

3200

I.S. Talke, J.J. Borkowski / Journal of Statistical Planning and Inference 142 (2012) 3189–3200

Hua, L.K., Wang, Y., 1981. Application of Number Theory to Numerical Analysis. Springer-Verlag, Science Press, Berlin, Beijing. Johnson, M.E., Moore, L.M., Ylvisaker, D., 1990. Minimax and maxmin distance designs. Journal of Statistical Planning and Inference 26 (2), 131–148. Koehler, J.R., Owen, A.B., 1996. Computer experiments. In: Ghosh, S., Rao, C.R. (Eds.), Handbook of Statistics, vol. 13. , Elsevier Science B.V.. Korobov, N.M., 1959. The approximate computation of multiple integrals. Doklady Akademii Nauk SSSR 124, 1207–1210. Li, R., Lin, D.K.J., Chen, Y., 2004. Uniform design: design, analysis and application. International Journal of Materials and Product Technology 20, 101–114. Liang, Y., Fang, K.T., Xu, Q., 2001. Uniform designs and its application in chemistry and chemical engineering. Chemometrics and Intelligent Laboratory Systems 58, 44–57. McKay, M.D., Beckman, R.J., Conover, W.J., 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), 239–245. Sobol, I.M., 1967. On the distribution of points in a cube and the approximate evaluation of integrals. USSR Computational Mathematics and Mathematical Physics 7, 86–112. Wang, Y., Fang, K.T., 1981. A note on uniform distribution and experimental design. KeXue TongBao 6, 485–489. Wang, Y., Fang, K.T., 1990. Number theoretic methods in applied statistics. Chinese Annals of Mathematics 11B (1), 51–65. Zhang, L., Liang, Y., Jiang, J., Yu, R., Fang, K.T., 1998. Uniform designs applied to nonlinear multivariate calibration by ANN. Analytica Chimica Acta 370 (1), 65–77.