Orthogonal Latin hypercube designs for Fourier-polynomial models

Journal of Statistical Planning and Inference 143 (2013) 307–314

Contents lists available at SciVerse ScienceDirect

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

Orthogonal Latin hypercube designs for Fourier-polynomial models Yuhui Yin, Min-Qian Liu n Department of Statistics, School of Mathematical Sciences and LPMC, Nankai University, Tianjin 300071, China

a r t i c l e i n f o

abstract

Article history: Received 20 October 2011 Received in revised form 30 July 2012 Accepted 3 August 2012 Available online 13 August 2012

Latin hypercube designs (LHDs) are widely used in computer experiments because of their one-dimensional uniformity and other properties. Recently, a number of methods have been proposed to construct LHDs with properties that all linear effects are mutually orthogonal and orthogonal to all second-order effects, i.e., quadratic effects and bilinear interactions. This paper focuses on the construction of LHDs with the above desirable properties under the Fourier-polynomial model. A convenient and flexible algorithm for constructing such orthogonal LHDs is provided. Most of the resulting designs have different run sizes from that of Butler (2001), and thus are new and very suitable for factor screening and building Fourier-polynomial models in computer experiments as discussed in Butler (2001). & 2012 Elsevier B.V. All rights reserved.

Keywords: Computer experiment Latin hypercube Fourier-polynomial model Orthogonal design

1. Introduction Computer experiments are frequently the most effective approach to probing the complex systems. When a simulation model of the system performance in a computer experiment is computationally expensive, a metamodel is often used to approximate the relationship between the system performance and the design parameters. In metamodeling, one basic task that the experimenter has to face is to select a suitable experimental design. The design of computer experiments is different from that of physical experiments because the output of a computer experiment is not subject to random variation. In a computer experiment, none of the traditional principles of blocking, randomization, and replication are of use in solving the design and analysis problems. Latin hypercube designs (LHDs) have been almost exclusively recommended in such conditions because they have good one-dimensional projective properties. An LHD for m factors in n runs can be denoted by an n  m matrix Lðn,mÞ, where the columns of the matrix are required to be permutations of n uniformly spaced numbers. Various efforts, resulted from different perspectives, have been made to optimize LHDs since they were first introduced by McKay et al. (1979). Owen (1992) and Tang (1993) proposed orthogonal array-based LHDs with the feature that they achieve stratification in low dimensions. Owen (1994) and Tang (1998) provided computational algorithms for searching nearly orthogonal LHDs. Park (1994) constructed LHDs that optimize the integrated mean squared error criterion. Morris and Mitchell (1995) investigated the maximin distance LHDs. Ye (1998) presented orthogonal LHDs by a full algebraic construction method, and Cioppa and Lucas (2007) extended Ye’s approach by adding new orthogonal columns to his orthogonal LHDs. Butler (2001) considered orthogonal LHDs with respect to Fourier-polynomial models. Beattie and Lin (2004, 2005) constructed LHDs by rotating the points in p-level full factorial designs, and then Steinberg and Lin (2006) and Pang et al. (2009) developed orthogonal LHDs by means of rotating factorial designs. Joseph and Hung (2008) proposed a multi-objective optimization approach to find good LHDs by

n

Corresponding author. Tel.: þ86 2223504709; fax: þ86 2223506423. E-mail address: [email protected] (M.-Q. Liu).

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

308

Y. Yin, M.-Q. Liu / Journal of Statistical Planning and Inference 143 (2013) 307–314

combining correlation and distance performance measures. Bingham et al. (2009), Lin et al. (2009, 2010) and Sun et al. (2011) proposed construction methods for orthogonal and nearly orthogonal LHDs. Recently, Georgiou (2009), Sun et al. (2009, 2010) and Yang and Liu (2012) provided flexible methods to construct orthogonal LHDs with properties that all linear effects are mutually orthogonal and orthogonal to all second-order effects, i.e., quadratic effects and bilinear interactions, under the polynomial regression model. Note that the LHDs constructed by Ye (1998) and Cioppa and Lucas (2007) also have such properties. In computer experiments, a variety of metamodeling techniques exist, e.g., the response surface methodology, Kriging method, artificial neural network method, multivariate adaptive regression spline and radial basis function approximation (see e.g., Santner et al., 2003; Fang et al., 2006). However, the most popular methods are perhaps the response surface methodology and Kriging method because of their easy explanations. These two methods have their respective advantages and disadvantages. The polynomial regression can be used to narrow the design variables to the most important ones in the situation with a large dimension, and the smoothing method capacity allows quick convergence of noise functions. But the drawback will stand out when applying polynomial regression to model highly nonlinear behaviors. In contrast, the Kriging model is more suitable for fitting the nonlinear systems (Jin et al., 2000). When we have no prior information about the system, the Fourier-polynomial model which was first proposed by Butler (2001) can be a good choice because of its balance between the polynomial model and the spatial model. Justifications for using the Fourier-polynomial model as a good approximation to the spatial model as well as the polynomial model in computer experiments were provided in Section 2 of Butler (2001). Meanwhile, he also showed that the correlations between the linear effects (under a polynomial model) of a newly constructed Lð7,6Þ and the correlations between the Fourier-linear effects (under a Fourier-polynomial model) of this LHD are similar, and are small (i.e. the effects are near-orthogonal) in both cases. The similarity in these two kinds of correlations illustrated the close relationship between linear effects and Fourier-linear effects, and the near-orthogonality is useful for factor screening in both cases. In fact, when the Fourier-polynomial model is used for metamodeling, LHDs with the following orthogonality properties: (a) all Fourier-linear effects are mutually orthogonal, (b) all Fourier-linear effects are orthogonal to all second-order effects, (c) the overall mean effect is orthogonal to all Fourier-linear effects and second-order effects, are suitable for factor screening, which is vital in computer experiments for reducing the dimension of the factor space before carrying out more detailed experimentation (cf. Butler, 2001). The LHDs constructed by Butler (2001) have the above orthogonality properties, where the number of runs, n, is assumed to be a prime. In this paper, we will construct a new class of LHDs with the same properties but with different run sizes. The paper is organized as follows. In Section 2, we introduce the second-order Fourier-polynomial model and provide a theoretical result that will play an important role in the construction. Section 3 presents an approach to constructing the orthogonal LHDs based on orthogonal designs. Some discussion and concluding remarks are given in Section 4. 2. Second-order Fourier-polynomial model Butler (2001) proposed to use the second-order Fourier-polynomial model as a good simulation model in computer experiments. A full second-order Fourier-polynomial model is of the form   pffiffiffi X       m m m 1 X m X pffiffiffi X pðxj 0:5Þ pðxi 0:5Þ 2pðxi 0:5Þ pðxi 0:5Þ þ 2 þ2 cos þ e, y ¼ b0  2 bi cos bii cos bij cos n n n n i¼1 i¼1 i ¼ 1 j ¼ iþ1 ð1Þ where e  Nð0, s2 Þ. Here e is required in order to model higher-order systematic effects. This model may also be expressed in the linear-model form as Y ¼ FðXÞT b þ e,

e  Nð0, s2 In Þ,

where X ¼ ðX ti Þ is the design with n runs and m factors, Y is the n  1 column vector of observations, b ¼ ðb0 , b1 , . . . , bm , b11 , . . . , bmm , b12 , . . . , bm1,m Þ0 denotes the vector of unknown parameters, and In is the n  n identity matrix. Furthermore, FðXÞ ¼ ð1n ,F L ,F Q ,F I Þ is the corresponding design matrix with 1n being the n  1 column vector with all elements unity, and   pffiffiffi pðX ti 0:5Þ , ð2Þ ðF L Þti ¼  2 cos n ðF Q Þti ¼

  pffiffiffi 2pðX ti 0:5Þ , 2 cos n

ðF I Þt,ði,jÞ ¼ 2 cos





pðX ti 0:5Þ n

cos

ð3Þ 

pðX tj 0:5Þ n

 ,

ð4Þ

Y. Yin, M.-Q. Liu / Journal of Statistical Planning and Inference 143 (2013) 307–314

309

for t ¼ 1, . . . ,n, 1 r i,j rm and ði,jÞ ¼ j þmði1Þiði þ1Þ=2 for i o j. Here FL denotes the design matrix of the Fourier-linear effects, FQ denotes the design matrix of the Fourier-quadratic effects and FI denotes the design matrix of the interactions between the Fourier-linear effects. As discussed in Butler (2001), model (1) is able to model the large-scale spatial trends and so may be interpreted as an approximation to spatial models, and this model is also a good approximation to the polynomial model of the form y ¼ b0 þ

m X

bi x i þ

i¼1

m X

bii x2i þ

i¼1

m 1 X

m X

bij xi xj þ e:

i ¼ 1 j ¼ iþ1

An LHD is said to have resolution IV if it satisfies properties (a), (b) and (c) (cf. Butler, 2001). Note that for an LHD, FL and FQ satisfy the constraints 1n 0 F L ¼ 1n 0 F Q ¼ 0, thus a resolution IV LHD requires that (or equivalently, properties (a), (b) and (c) become) (a) F L 0 F L ¼ nIm , (b) F L 0 F Q ¼ 0,F L 0 F I ¼ 0, (c) 1n 0 F I ¼ 0. Note that F L 0 F L ¼ nIm implies 1n 0 F I ¼ 0, thus an LHD with properties (a) and (b) is a resolution IV design. Butler (2001) used the Williams transformation to construct resolution IV LHDs with prime numbers of runs. Before presenting our new construction method for orthogonal LHDs in the subsequent section, we provide the following important result for the construction. Theorem 1. Suppose l(x) and q(x) are two functions satisfying lðxÞ ¼ lðn þ 1xÞ

and

ð5Þ

qðxÞ ¼ qðn þ 1xÞ

for any x 2 R,

ð6Þ

respectively, and X ¼ ðX ti Þ is a design with n runs and m factors. If X has the form 0 1 ! B B B n þ 1 1 0, C X¼ or X ¼ @ 2 m A for some B, ðn þ 1ÞB ðn þ1ÞB

ð7Þ

then X satisfies that n X

lðX ti ÞqðX tj Þ ¼ 0

for 1 r i,j rm

and

ð8Þ

t¼1

n X

lðX ti ÞlðX tj ÞlðX tk Þ ¼ 0

for 1r i,j,k rm:

ð9Þ

t¼1

Proof. Since for the design X in (7), lðX tj Þ ¼ lðX bðn þ 1Þ=2c þ t,j Þ

and

qðX tj Þ ¼ qðX bðn þ 1Þ=2c þ t,j Þ,

for t ¼ 1, . . . ,bn=2c and 1 r j rm, where bzc means the integer part of z, Eqs. (8) and (9) can then be obtained directly.

&

Now, we turn to the Cosine functions that define the design matrices F L ,F Q and FI in (2), (3) and (4), respectively. It can be easily checked that cosfpðx0:5Þ=ng and cosf2pðx0:5Þ=ng satisfy (5) and (6), respectively, i.e.,     pðx0:5Þ pððn þ 1xÞ0:5Þ ¼ cos cos n n     2pðx0:5Þ 2pððn þ 1xÞ0:5Þ and cos ¼ cos : n n Then, from Theorem 1, we have Corollary 1. For a given design X of the form (7), F L 0 F Q ¼ 0 and F L 0 F I ¼ 0.

310

Y. Yin, M.-Q. Liu / Journal of Statistical Planning and Inference 143 (2013) 307–314

Remark 1. Note that when the X in (7) is an LHD with levels from 1 to n, it is called a symmetric LHD (cf. Ye et al., 2000). Here, the X in Theorem 1 or Corollary 1 is not restricted to be a symmetric LHD, it can be any design of the form (7), even not an LHD.

3. New construction method The designs constructed by Butler (2001) have prime numbers of runs. Here we present a new class of LHDs with the above properties but with more flexible run sizes based on a general kind of orthogonal designs. An l  m orthogonal design on fy1 , . . . ,yl g, denoted by ODðl,m; 1, . . . ,1Þ, is an l  m matrix A that satisfies A0 A ¼

l X

ðy2i ÞIm :

i¼1

The entries of A are furthermore, restricted to be of the form 7 yi , i ¼ 1, . . . ,l, where the yi’s are variables that take values from the set of real numbers R, i.e., ðy1 , . . . ,yl Þ 2 Rl . The matrix 0

y1 B y B 2 B B y4 B B B y3 A¼B B y6 B B B y5 B B y 8 @ y7

y2

y4

y3

y6

y5

y8

y1

y3

y4

y5

y6

y7

y3

y1

y2

y8

y7

y6

y4 y5

y2 y8

y1 y7

y7 y1

y8 y2

y5 y4 y3

y6

y7

y8

y2

y1

y7

y6

y5

y4

y3

y1

y8

y5

y6

y3

y4

y2

y7

1

y8 C C C y5 C C C y6 C C y3 C C C y4 C C y2 C A

ð10Þ

y1

is thus an example of an 8  8 orthogonal design (cf. Georgiou, 2009). The construction algorithm works as follows. Algorithm 1. Step1. Step2. Step3.

Give an ODðl,m; 1, . . . ,1Þ on fy1 , . . . ,yl g, denoted by A, and assume that y1 , . . . ,yl are all positive variables. Set D1 ¼ ðA0 ,A0 Þ0 and D2 ¼ ðA0 ,0m ,A0 Þ0 , then two matrices with m columns, n ¼ 2 l and n ¼ 2 l þ1 rows are generated, respectively. Obtain two symmetric LHDs X1 and X2 by applying a suitably selected one-to-one transformation that transforms the entries of D1 and D2 to 1, . . . ,n, respectively.

Remark 2. The transformation used in Algorithm 1 is not unique, for example, it can be any of the following four transformations: 8 if t ¼ yi > : ðn þ 1Þi if t ¼ y i ( (ii) ði1=2Þ þ ðn þ 1Þ=2 if t ¼ yi for t 2 D1 ,x 2 X 1 and xðtÞ ¼ ði1=2Þ þ ðn þ 1Þ=2 if t ¼ yi 8 if t ¼ yi > < i þ ðn þ 1Þ=2 if t ¼ 0 xðtÞ ¼ ðn þ 1Þ=2 for t 2 D2 ,x 2 X 2 , respectively; > : i þ ðn þ1Þ=2 if t ¼ y i (i)

(iii) a one-to-one transformation that first adds ðn þ1Þ to each yi for i ¼ 1, . . . ,l, and then replaces variables y1 , . . . ,yl by a permutation of f1, . . . ,lg, and replaces 0 with ðn þ 1Þ=2 if it exists; (iv) a one-to-one transformation that replaces variables y1 , . . . ,yl by a permutation of f0:5,1:5, . . . ,l0:5g for even n, or a permutation of f1, . . . ,lg for odd n, and then add ðn þ1Þ=2 to all the entries. Note that transformations (i) and (ii) are specific cases of (iii) and (iv), respectively, and (iii) and (iv) may produce the same symmetric LHDs. Here is an illustrative example for Algorithm 1.

Y. Yin, M.-Q. Liu / Journal of Statistical Planning and Inference 143 (2013) 307–314

311

Example 1. Use the orthogonal design A in (10), we can construct an orthogonal LHD with n¼16 runs and m¼ 8 factors by Algorithm 1 and transformation (iv). First, we have 0 1 y2 y4 y3 y6 y5 y8 y7 y1 B y y1 y3 y4 y5 y6 y7 y8 C B C 2 B C B y4 y3 y1 y2 y8 y7 y6 y5 C B C B y y4 y2 y1 y7 y8 y5 y6 C B C 3 B C B y6 y5 y8 y7 y1 y2 y4 y3 C B C B y y6 y7 y8 y2 y1 y3 y4 C B C 5 B C B y8 y7 y6 C y y y y y 5 4 3 1 2 C B C   B y8 y5 y6 y3 y4 y2 y1 C B y7 A C D1 ¼ ¼B B y1 y2 y4 y3 y6 y5 y8 y7 C: A B C B C y1 y3 y4 y5 y6 y7 y8 C B y2 B C B y y3 y1 y2 y8 y7 y6 y5 C B 4 C B C B y3 y4 y2 y1 y7 y8 y5 y6 C B C B y y5 y8 y7 y1 y2 y4 y3 C B 6 C B C B y5 y6 y7 y8 y2 y1 y3 y4 C B C B y y7 y6 y5 y4 y3 y1 y2 C @ 8 A y7 y8 y5 y6 y3 y4 y2 y1 Then applying transformation (iv) with the permutation of f0:5,1:5, . . . ,7:5g being (6.5, 0.5, 2.5, 5.5, 1.5, 4.5, 3.5,7.5) to the entries of D1, i.e., the one-to-one transformation t

y1

y2

y3

y4

y5

y6

y7

y8

y1

y2

y3

y4

y5

y6

y7

y8

x

2

8

6

3

7

4

5

1

15

9

11

14

10

13

12

16

we get the following 0 15 8 B 9 15 B B B 14 11 B B 11 3 B X¼B B 13 10 B B 10 4 B B @ 16 12 12 1

,

LHD for the Fourier-polynomial model: 3

6

4

7

6 15

1

5

2

9

14

11

13

14

7

8

16

13

5

16

8

2

11

3

5

4

10

3

6

2

9

10

10

16

12

10

4

12

1

12

13 7

1 C C C 7 C C 4 C C C: 11 C C 14 C C C 9 A

9

15

5

1

10

13

6

14

8

2

12

16

1

12

15

8

14

6

4

7

16

5

2

9

3

12

16

9

15

6

3

7

13

5

1

8

2

11

13 7

7 4

3 11

11 14

15 9

8 15

1 5

5 16

4 10

10 13

14 6

6 3

2 8

2

It can be easily checked that X is a symmetric LHD satisfying properties (a) and (b). Based on the inherent orthogonal feature of A that the orthogonality is irrelevant to the specific values of yi’s, the foldover structures of D1 and D2, and Theorem 1, we can easily show Theorem 2 below. Theorem 2. The two symmetric LHDs X1 and X2 constructed in Algorithm 1 by applying any of the transformations given in Remark 2 satisfy properties (a) and (b), i.e., they have resolution IV. Note that both X1 and X2 constructed from an ODðm,m; 1, . . . ,1Þ satisfy that m ¼ bn=2c. In fact, the number of factors m in X1 or X2 constructed from an ODðm,m; 1, . . . ,1Þ attains its maximum value among all the corresponding resolution IV LHDs. This conclusion is based on the following important theorem. Theorem 3. If X is an n  m resolution IV LHD with entries f1, . . . ,ng, then m rbn=2c. Proof. For X, let f j ¼ ðf 1j , . . . ,f nj Þ0 be the jth column of the corresponding matrix of Fourier-linear effects FL, j ¼ 1, . . . ,m, V 1 ¼ Spanf1n ,f 1 , . . . ,f m g and V 2 ¼ Spanff 1  f 2 , . . . ,f 1  f m g, where Spanf1n ,f 1 , . . . ,f m g denotes a linear space spanned by vectors 1n ,f 1 , . . . ,f m . Note that f 1  f j for j ¼ 2, . . . ,m are columns of the corresponding FI, and X satisfies properties (a), (b) ? and (c), we have V 2 D V ? means the orthogonal complement space of V. In addition, dimðV 2 Þ ¼ m1, where 1 , where V P dimðVÞ denotes the dimension of V. In fact, if there exist m1 numbers ej for j ¼ 2, . . . ,m, satisfying m j ¼ 2 ej f 1  f j ¼ 0, i.e., 2 3 m X ej f j 5 ¼ 0, f1  4 j¼2

then (i) if n is even, we can know from the Cosine function in (2) that f t1 a0 for t ¼ 1, . . . ,n, thus

Pm

j¼2

ej f j ¼ 0,

312

Y. Yin, M.-Q. Liu / Journal of Statistical Planning and Inference 143 (2013) 307–314

(ii) if n is odd, there is a unique zero element in f1, i.e., the one corresponding to X tj ¼ ðn þ 1Þ=2 for some t. Without loss of P 0 0 generality, suppose f 11 ¼ 0, then m j ¼ 2 ej f j ¼ ðg,0, . . . ,0Þ , furthermore, from 1n f j ¼ 0 we can obtain g¼0, thus we have

Pm

j¼2

ej f j ¼ 0, and the orthogonality of fj’s implies that ej ¼ 0 for j ¼ 2, . . . ,m, so f 1  f 2 , . . . ,f 1  f m are

linearly independent, i.e., dimðV 2 Þ ¼ m1. Then dimðV 1 þ V 2 Þ ¼ dimðV 1 Þ þdimðV 2 Þ ¼ 2m rn, this completes the proof.

&

Based on orthogonal designs ODðm,m; 1, . . . ,1Þ, we can construct resolution IV LHDs which attain the maximum numbers of factors. Next, we provide an approach to generate ODð2c1 ,2c2; 1, . . . ,1Þ’s. Ye (1998) suggested an algorithm for constructing orthogonal LHDs with n ¼ 2c or 2c þ 1 runs and m ¼ 2c2 columns, where c is a positive integer. Based on his construction method, we now generalize Ye’s algorithm to construct a kind of orthogonal designs with n ¼ 2c1 runs and m ¼ 2c2 columns. Algorithm 2. I  ffl{zfflfflfflfflfflffl   ffl}I  Q      Q , k ¼ 1, . . . ,c1, where Let u ¼ ðy1 , . . . ,y2c1 Þ0 be a 2c1  1 vector and Ak ¼ |fflfflfflfflfflffl |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}     1 0 0 1 ck1 k I¼ , Q¼ , 0 1 1 0

Step1.

and  denotes the Kronecker product. Let M ¼ fu,Ak u,Ac1 Aj u,k ¼ 1, . . . ,c1,j ¼ 1, . . . ,c2g. Define hk ¼ q1      qc1 , where qck ¼ ð1,1Þ0 , qi ¼ ð1,1Þ0 , iack, and set H ¼ ð1,hi ,h1  hj ,i ¼ 1, . . . , c1, j ¼ 2, . . . ,c1Þ, where  denotes the elementwise product. Let T ¼ M  H. Then T is an ODð2c1 ,2c2; 1, . . . ,1Þ (cf. Ye, 1998).

Step2. Step3. Step4.

Apply any of the transformations in Remark 2 to the entries of D1 ¼ ðT 0 ,T 0 Þ0 and D2 ¼ ðT 0 ,0m ,T 0 Þ0 , respectively, then we get two symmetric LHDs with n ¼ 2c and 2c þ 1 runs, and m ¼ 2c2 factors, denoted by X1 and X2, respectively. Example 2. For c¼3, we can obtain D1 and D2 as follows: 0 10 y2 y3 y4 y1 y2 y3 y4 y1 B C y1 y4 y3 y2 y1 y4 y3 C B y2 C and D1 ¼ B B y4 y3 y2 y1 y4 y3 y2 y1 C @ A y3 y4 y1 y2 y3 y4 y1 y2 0 10 y2 y3 y4 0 y1 y2 y3 y4 y1 B C y1 y4 y3 0 y2 y1 y4 y3 C B y2 C: D2 ¼ B B y4 y3 y2 y1 0 y4 y3 y2 y1 C @ A y3 y4 y1 y2 0 y3 y4 y1 y2 On applying transformation (iii) with the permutation of f1,2,3,4g being ð4,2,3,1Þ to the entries of D1, and with the permutation of f1,2,3,4g being ð3,1,2,4Þ to the entries of D2, respectively, we have two LHDs with resolution IV: 0

10

4 B7 B X1 ¼ B @8

2

3

1

5

7

6

8

4

8

3

2

5

1

6

2

4

1

3

7

6C C C 5A

3

8

5

2

6

1

4

7

3

1

2

4

5

7

9

8

6

B9 B X2 ¼ B @6

3

6

2

5

1

7

4

8

1

3

5

4

2

9

8C C C: 7A

2

6

7

1

5

8

4

3

9

0

and 10

It can be seen that the one-to-one transformations are given by t

y1

y2

y3

y4

y1

y2

y3

y4

x

5

7

6

8

4

2

3

1

t

y1

y2

y3

y4

0

y1

y2

y3

y4

x

7

9

8

6

5

3

1

2

4

for n ¼ 8

and

for n ¼ 9

respectively:

Y. Yin, M.-Q. Liu / Journal of Statistical Planning and Inference 143 (2013) 307–314

313

Table 1 Existence of orthogonal LHDs with nr 33. Run size n

Maximal number of factors m Butler (2001)

4 5 7 8 9 11 13 16 17 19 23 29 31 32 33

2 3

ODð2c1 ,2c2; 1, . . . ,1Þ

ODðm,m; 1, . . . ,1Þ

2 2

2 2

4 4

4 4

6 6

8 8

5 6 8 9 11 14 15

8 8

4. Discussion and concluding remarks In the previous section, we have introduced a method for constructing orthogonal LHDs. The resulting LHDs are resolution IV designs, i.e., they have the appealing properties (a), (b) and (c), and most of these designs have different run sizes from that of Butler (2001), as we can see in Table 1. In particular, the LHDs constructed by the ODð2c1 ,2c2; 1, . . . ,1Þ’s from Algorithm 2 have m ¼ 2c2 factors, and n ¼ 2c or 2c þ 1 runs, and the LHDs constructed by ODðm,m; 1, . . . ,1Þ attain the maximum numbers of factors. Note that the LHDs with the maximum numbers of factors depend on the existence of ODðm,m; 1, . . . ,1Þ, the three ODs used in Table 1 are ! y1 y2 ODð2,2; 1,1Þ : ðGeramita and Seberry;1979Þ; y2 y1 0

y1 B B y2 ODð4,4; 1,1,1,1Þ : B B y3 @ y4

y2

y4

y1

y3

y4

y2

y3

y1

y3

1

C y4 C C y1 C A y2

ðGeramita and Seberry;1979Þ,

and the ODð8,8; 1, . . . ,1Þ is given in (10). From Table 1, it can be seen that there are some run sizes with which both Butler’s approach and the newly proposed method can be used. For such cases, we suggest examining the correlations between pairs of quadratic effects and pairs of bilinear interactions. By comparing the averages and maximums of absolute correlations, we can choose the best one similarly as Jones and Nachtsheim (2011) did for selecting screening designs. Reading the results in Table 1, the proposed method can also provide designs with the same number of factors and different run sizes. For such cases, we can use the cost-saving designs, i.e., those with larger factor-to-run ratios. Note that for the designs listed in Table 1, some of Butler’s (2001) have larger factor-to-run ratios, and some of the newly constructed designs have larger factor-to-run ratios. For example, for designs with m¼6 factors, the one with n¼13 runs of Butler (2001) has a larger factor-to-run ratio, while for the designs with m ¼ 2,4 and 8 factors, the newly constructed ones with n ¼ 4,8 and 16 runs have larger factor-to-run ratios, respectively. If we need to choose a design for six factors, we can use the design with 13 runs provided by Butler (2001) instead of that with 16 or 17 runs; while if a design for eight factors is needed, the new design with 16 runs constructed by the ODð8,8; 1, . . . ,1Þ can be used. Thus, considering the factor-to-run ratio, the proposed method can produce designs with larger ratios, which supplement the results of Butler (2001). As discussed in Butler (2001), the newly constructed LHDs are also appropriate for factor screening as they require very few experimental runs per factor, and they are very useful for building Fourier-polynomial models.

Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant No. 10971107), the ‘‘131’’ Talents Program of Tianjin, and the Fundamental Research Funds for the Central Universities (Grant No. 65030011). The authors thank the Executive Editors and two referees for their valuable comments.

314

Y. Yin, M.-Q. Liu / Journal of Statistical Planning and Inference 143 (2013) 307–314

References Beattie, S.D., Lin, D.K.J., 2004. Rotated factorial design for computer experiments. Journal of the Chinese Statistical Association 42, 289–308. Beattie, S.D., Lin, D.K.J., 2005. A new class of Latin hypercube for computer experiments. In: Fan, J., Li, G. (Eds.), Contemporary Multivariate Analysis and Experimental Designs in Celebration of Professor Kai-Tai Fang’s 65th Birthday, World Scientific, Singapore, pp. 206–226. Bingham, D., Sitter, R.R., Tang, B., 2009. Orthogonal and nearly orthogonal designs for computer experiments. Biometrika 96, 51–65. Butler, N.A., 2001. Optimal and orthogonal Latin hypercube designs for computer experiments. Biometrika 88, 847–857. Cioppa, T.M., Lucas, T.W., 2007. Efficient nearly orthogonal and space-filling Latin hypercubes. Technometrics 49, 45–55. Fang, K.T., Li, R., Sudjianto, A., 2006. Design and Modeling for Computer Experiments. Chapman & Hall, Boca Raton. Georgiou, S.D., 2009. Orthogonal Latin hypercube designs from generalized orthogonal designs. Journal of Statistical Planning and Inference 139, 1530–1540. Geramita, A.V., Seberry, J., 1979. Orthogonal Designs: Quadratic Forms and Hadamard Matrices. Marcel Dekker, New York, Basel. Jin, R., Chen, W., Simpson, T.W., 2000. Comparative studies of metamodeling techniques under multiple modeling criteria. In: 8th AIAA/NASA/USAF/ ISSMO Symposium on Multidisciplinary Analysis and Optimization, Long Beach, CA, Paper 4801. Jones, B., Nachtsheim, C.J., 2011. A class of three-level designs for definitive screening in the presence of second-order effects. Journal of Quality Technology 43, 1–15. Joseph, V.R., Hung, Y., 2008. Orthogonal-maximin Latin hypercube designs. Statistica Sinica 18, 171–186. Lin, C.D., Bingham, D., Sitter, R.R., Tang, B., 2010. A new and flexible method for constructing designs for computer experiments. Annals of Statistics 38, 1460–1477. Lin, C.D., Mukerjee, R., Tang, B., 2009. Construction of orthogonal and nearly orthogonal Latin hypercubes. Biometrika 96, 243–247. 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, 239–245. Morris, M.D., Mitchell, T.J., 1995. Exploratory designs for computational experiments. Journal of Statistical Planning and Inference 43, 381–402. Owen, A.B., 1992. Orthogonal arrays for computer experiments, integration and visualization. Statistica Sinica 2, 439–452. Owen, A.B., 1994. Controlling correlations in Latin hypercube samples. Journal of the American Statistical Association 89, 1517–1522. Pang, F., Liu, M.Q., Lin, D.K.J., 2009. A construction method for orthogonal Latin hypercube designs with prime power levels. Statistica Sinica 19, 1721–1728. Park, J.S., 1994. Optimal Latin-hypercube designs for computer experiments. Journal of Statistical Planning and Inference 39, 95–111. Santner, T., Williams, B., Notz, W., 2003. The Design and Analysis of Computer Experiments. Springer, New York. Steinberg, D.M., Lin, D.K.J., 2006. A construction method for orthogonal Latin hypercube designs. Biometrika 93, 279–288. Sun, F.S., Liu, M.Q., Lin, D.K.J., 2009. Construction of orthogonal Latin hypercube designs. Biometrika 96, 971–974. Sun, F.S., Liu, M.Q., Lin, D.K.J., 2010. Construction of orthogonal Latin hypercube designs with flexible run sizes. Journal of Statistical Planning and Inference 140, 3236–3242. Sun, F.S., Pang, F., Liu, M.Q., 2011. Construction of column-orthogonal designs for computer experiments. Science China Mathematics 54, 2683–2692. Tang, B., 1993. Orthogonal array-based Latin hypercubes. Journal of the American Statistical Association 88, 1392–1397. Tang, B., 1998. Selecting Latin hypercubes using correlation criteria. Statistica Sinica 8, 965–977. Yang, J.Y., Liu, M.Q., 2012. Construction of orthogonal and nearly orthogonal Latin hypercube designs from orthogonal designs. Statistica Sinica 22, 433–442. Ye, K.Q., 1998. Orthogonal column Latin hypercubes and their application in computer experiments. Journal of the American Statistical Association 93, 1430–1439. Ye, K.Q., Li, W., Sudjianto, A., 2000. Algorithmic construction of optimal symmetric Latin hypercube designs. Journal of Statistical Planning and Inference 90, 145–159.