- Email: [email protected]
An algorithm for calculating D-optimal designs for polynomial regression through a fixed point
Journal of Statistical Planning and Inference 142 (2012) 935–943
Contents lists available at SciVerse ScienceDirect
Journal of Statistical Planning and Inference journal homepage: www.elsevier.com/locate/jspi
An algorithm for calculating D-optimal designs for polynomial regression through a fixed point Hiroto Sekido Department of Applied Mathematics and Physics, Graduate School of Informatics, Kyoto University, Kyoto 606-8501, Japan
a r t i c l e i n f o
abstract
Article history: Received 11 January 2010 Received in revised form 10 September 2011 Accepted 31 October 2011 Available online 10 November 2011
Optimal designs are required to make efficient statistical experiments. By using canonical moments, in 1980, Studden found Ds-optimal designs for polynomial regression models. On the other hand, integrable systems are dynamical systems whose solutions can be written down concretely. In this paper, polynomial regression models through a fixed point are discussed. In order to calculate D-optimal designs for these models, a useful relationship between canonical moments and discrete integrable systems is introduced. By using canonical moments and discrete integrable systems, a new algorithm for calculating D-optimal designs for these models is proposed. & 2011 Elsevier B.V. All rights reserved.
Keywords: D-optimal design Canonical moment Polynomial regression model Integrable system Nonautonomous discrete Toda equation
1. Introduction For a given probability measure on the Borel sets of a finite interval, we can define canonical moments in addition to ordinary moments. Canonical moments are defined as normalized moments in some sense. For example, canonical moments pk satisfy 0 rpk r 1. Other properties of canonical moments were studied by Skibinsky (1967, 1968, 1969). Applications of canonical moments are explained in Dette and Studden (1997). One of the applications is a design of experiments. In this paper, we consider the D-optimal designs which corresponds to an objective function of determinant form. D-optimal designs and Ds-optimal designs for polynomial regression models were studied by Studden (1980). In the calculation by Studden (1980), an important point is to use canonical moments instead of ordinary moments, and D-optimal designs can be identified by its canonical moments. As we see in Section 2, the objective function is written down in term of canonical moments. Besides this, D-optimal designs for various models have been calculated by Biedermann et al. (2009), Dette (1992), Dette and Franke (2000, 2001), Dette et al. (2007), Fang (2002, 2006), Fang and Wiens (2003), Huang et al. (1995), Lau (1988), Lau and Studden (1985), Lim and Studden (1988), Studden (1982b, 1982a, 1989), Zen and Tsai (2004). For example, D-optimal designs for weighted polynomial regression with a weight function xa ð1xÞb were found by Studden (1982a) by using canonical moments. However, in most cases, an explicit form of the D-optimal design is unknown. We here consider polynomial regression models with some prior information. For example, D-optimal designs for this kind of models has been studied by Dette (1992), Fang (2002), Huang et al. (1995). Dette (1992) calculated D-optimal designs for polynomial regression models with only odd (or even) degree terms. Huang et al. (1995) considered polynomial regression without intercept. Fang (2002) considered polynomial regression models through origin.
E-mail address: [email protected] 0378-3758/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2011.10.014
936
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
On the other hand, the term integrable systems is used for nonlinear dynamical systems whose solutions can be written down concretely. For Hamiltonian systems with finite degree of freedom, their integrability is defined in the Liouville– Arnold theorem (Arnold, 1978). However, even now, there is no mathematical definition of integrability for nonlinear systems with infinite degree of freedom. This is why we call an explicitly solvable nonlinear system as an integrable system. For example, the Lotka–Volterra equation: dU k ¼ U k ðU k þ 1 U k1 Þ, dt
ð1Þ
known as a model equation for a food chain is one of integrable systems. Discrete integrable systems are discrete analogues of integrable systems. For example, the discrete Lotka–Volterra equation (Spiridonov and Zhedanov, 1997): þ 1Þ þ 1Þ ð1þ duðn Þ ¼ uðnÞ ð1 þ dukðnÞþ 1 Þ uðn k k1 k
ð2Þ ukðnÞ
indicates the value of uk at the time nd, is one of discrete integrable systems, where the integer n is a discrete time and d is a discrete time step. Note that the discrete Lotka–Volterra equation (2) reduces in the continuous limit d-0 with nd ¼ T to the Lotka–Volterra equation (1). Both (1) and (2) have explicit solutions of determinant form. Here, we approach to a relationship between canonical moments and determinant solutions of discrete integrable systems. This is a new viewpoint in optimal designs. In this paper, we propose an algorithm for calculating D-optimal design for polynomial regression through a given fixed point in x ¼ l, where l is an arbitrary real number. We here consider the design space is X ¼ ½0; 1. When l ¼ 0 and l ¼ 1, the corresponding D-optimal designs are equivalent to D-optimal designs for weighted polynomial regression with a weight function x2 and ð1xÞ2 , respectively. These are special cases of D-optimal designs calculated in Studden (1982a). Fang (2002) proposed the D-optimal designs for polynomial regression through origin with the design space X ¼ ½1; 1. It contains the case of l ¼ 1=2. While in this paper, we consider polynomial regression through an arbitrary fixed point. See Huang et al. (1995) for the motivation in terms of regression models. The resulting algorithm for calculating D-optimal designs can be regarded as a new application of discrete integrable systems. Integrable systems have been applied to numerical analysis. Typical examples are matrix eigenvalue algorithms (Nakamura, 1988; Symes, 1982), and algorithms for computing matrix singular values (Tsujimoto et al., 2001) in terms of the discrete Lotka–Volterra equation (2). Section 3 shows that canonical moments can be rephrased as determinant solutions of (2). This is the starting point of our approach. Therefore, discrete integrable systems, which are generalizations of (2), help us to write down the objective function in terms of canonical moments for polynomial regression through a fixed point. In Section 2, we introduce the notion of D-optimal designs and give a short review on the calculation of D-optimal designs by Studden (1980). Section 3 indicates a relationship between canonical moments and discrete integrable systems at first. Then, an algorithm for calculating D-optimal designs for polynomial regression through a fixed point is proposed. In Section 4, we calculate D-optimal designs for a second degree polynomial regression through a fixed point in terms of the nonautonomous discrete Toda equation, an integrable system. Section 5 is devoted to conclusions. 2. Polynomial regression models and its D-optimal designs: a preliminary In this section, we give a brief review on the D-optimal designs for polynomial regression models. Let us call polynomial regression models PRMs, for short. After that, we introduce canonical moments and give a review on the calculation of D-optimal designs by Studden (1980). Consider common (m1)th degree PRMs: Y ¼ ym1 xm1 þ ym2 xm2 þ þ y0 þ e, where yk is an unknown parameter, and e is a random error term. Let P ½0;1 denote the set of all probability measures on the Borel sets of the interval ½0; 1. For given m 2 P ½0;1 , let ck be the kth moment given by Z 1 xk dmðxÞ, k ¼ 0; 1,2, . . . : ck ¼ 0
The D-optimal designs are defined as probability measures in P ½0;1 which maximize the determinant of the Fisher information matrix, that is, D-optimal designs m 2 P ½0;1 for PRMs are defined as the optimal solution of the optimization problem maximize
Hm ðmÞ,
subject to
m 2 P ½0;1 ,
ð3Þ
where let Hk ðmÞ be the Hankel determinants of moments as m1
Hm ðmÞ ¼ 9ci þ j 9i,j ¼ 0 : See Pukelsheim (1993) for more details about optimal designs.
ð4Þ
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
937
The objective function Hm of (3) is explicitly expressed in terms of the moments ck. However, in terms of moments ck, the set P ½0;1 of feasible solutions of (3) is expressed in a very complicated form. Therefore, we use canonical moments instead of the ordinary moments ck. This method was proposed by Studden (1980) at first, and is explained in the book by Dette and Studden (1997). Now, we define canonical moments. For a given probability measure m 2 P ½0;1 , let ckþ denote the maximum of the k-th moment over the set of all measure having moments c0 ,c1 , . . . ,ck1 . Similarly let c k denote the corresponding minimum. Canonical moments are defined by pk ¼
ck c k , ckþ c k
k ¼ 1; 2, . . . ,N,
ð5Þ
þ where N is the minimum of j which satisfy cjþþ 1 ¼ c j þ 1 . If c j 4 cj for an arbitrary positive integer j, let N be 1. Canonical moments have the property that ( pk 2 ð0; 1Þ, k ¼ 1; 2, . . . ,N1,
pN 2 f0; 1g:
ð6Þ
Conversely, for a given arbitrary sequence fpk gN k ¼ 1 which satisfies (6), there is a unique m 2 P ½0;1 which has the canonical moments fpk gN k ¼ 1. It is to be noted that canonical moments have a Hankel determinant expression. For introduction, we generalize the ðnÞ definition (4) of Hankel determinants Hk ðmÞ. Let HðnÞ ,H k be the Hankel determinants defined by k k1
HkðnÞ ¼ 9ci þ j þ n 9i,j ¼ 0 ,
k1
ðnÞ
H k ¼ 9ci þ j þ n ci þ j þ n þ 1 9i,j ¼ 0 ,
ð7Þ
ðnÞ where k ¼ 1; 2, . . . , n ¼ 0; 1,2, . . . . We here set that HðnÞ 0 ¼ 1, and that H k ¼ 0 if a matrix size k is negative. Obviously, we ð0Þ have Hk ¼ Hk . The canonical moments pk have the following Hankel determinant expression: ð1Þ
ð0Þ
p2k1 ¼
Hð1Þ H k1 k ð1Þ Hkð0Þ H k1
,
p2k ¼
Hð0Þ H k þ 1 k1 ð0Þ Hð1Þ Hk k
,
k ¼ 1; 2, . . . :
We introduce zk by the transformation of canonical moments pk
z0 ¼ 0, z1 ¼ p1 , zk ¼ ð1pk1 Þpk , k ¼ 2; 3, . . . ,N:
ð8Þ
Then zk also have the Hankel determinant expression:
z2k1 ¼
Hð1Þ Hð0Þ k k1 Hð0Þ Hð1Þ k k1
,
z2k ¼
Hkð0Þþ 1 Hð1Þ k1 Hð1Þ Hð0Þ k k
,
k ¼ 1; 2, . . . :
From this expression, Hankel determinants Hm can be expressed in a product of zk or of canonical moments pk, Hm ¼
m1 Y
ðz2k1 z2k Þmk
k¼1
0 ¼@
m1 Y
1 A ð1p2j Þmj1 pmj 2j
j¼1
m1 Y
ðð1p2j1 Þp2j1 Þmj ,
ð9Þ
j¼1
where 2 m2 r N. The expression (9) enable us to find the canonical moments pk which maximize the objective function Hm of (3) by a simple calculation. Consequently, we obtain the canonical moments pk of D-optimal designs for PRMs p2j1 ¼
1 , 2
p2j ¼
mj , 2ðmjÞ1
j ¼ 1; 2, . . . ,m1:
This approach for D-optimal designs is described in the book by Dette and Studden (1997). 3. D-optimal designs for the polynomial regression models through a fixed point In this section, we consider the mth degree polynomial regression models through a known point. In this case, we know a value of Eðy9x ¼ lÞ as a prior information, where l is an arbitrary given real number. Let us denote an mth degree polynomial regression model through a fixed point in x ¼ l as PRMm ðlÞ. In Section 3.2, PRMm ðlÞ and the corresponding D-optimal designs will be formulated as D-optimal design for a linear regression model. By Studden (1980), the objective function Hm ðmÞ of (3) is expressed in terms of canonical moments. In Section 3.3, we propose a method to express the objective function corresponding to PRMm ðlÞ in terms of canonical moments.
938
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
3.1. Canonical moments and discrete integrable systems We indicate a relationship between canonical moments and discrete integrable systems, before considering a formulation of the PRMm ðlÞ. In order to obtain a D-optimal design for PRMm ðlÞ, we can get clues from discrete integrable systems. Hankel determinants are related to both integrable systems and D-optimal designs, furthermore (some generalized) canonical moments pk and the variables zk satisfy a kind of discrete integrable systems. be defined by Proposition 3.1. Let pðnÞ k ðn þ 1Þ
ðnÞ ¼ p2k
HðnÞ H k þ 1 k1
ðnÞ
þ 1Þ Hðn Hk k
ðnÞ
,
pðnÞ ¼ 2k þ 1
þ 1Þ Hðn Hk kþ1
ðn þ 1Þ
HðnÞ H kþ1 k
,
k,n ¼ 0; 1,2, . . . :
satisfy the analogue of the discrete Lotka–Volterra equation: Then pðnÞ k þ 1Þ þ 1Þ ðnÞ ð1pðn Þð1 þ pðn Þ ¼ ð1pðnÞ Þð1þ p2k Þ, 2k þ 1 2k 2k þ 1 þ2 þ 1Þ ðn þ 1Þ ðnÞ p ¼ p2k pðnÞ , pðn 2k þ 2 2k þ 1 þ 2 2k þ 3
Let
zðnÞ k
k,n ¼ 0; 1,2, . . . :
ð10Þ
be defined by
zðnÞ 2k ¼
þ 1Þ HkðnÞþ 1 Hðn k1 þ 1Þ ðnÞ Hðn Hk k
,
ðnÞ z2k þ1 ¼
þ 1Þ ðnÞ Hðn Hk kþ1
HðnÞ Hðn þ 1Þ kþ1 k
,
k,n ¼ 0; 1,2, . . . :
ðnÞ
Then zk satisfy the discrete Toda equation (Hirota, 1997; Spiridonov and Zhedanov, 1995) þ 1Þ ðn þ 1Þ ðnÞ ðnÞ ðn þ 1Þ ðn þ 1Þ ðnÞ ðnÞ zðn þ z2k þ 1 ¼ z2k þ 1 þ z2k þ 2 , z2k þ 1 z2k þ 2 ¼ z2k þ 2 z2k þ 3 , 2k
ð11Þ
where k,n ¼ 0; 1, . . . . ð0Þ
¼ pk , and that zk ¼ zk . This proposition can be prove by using the determinant solution of integrable It is clear that pð0Þ k systems and Sylvester’s determinant identity. Note that (10) is obtained by substituting ðnÞ ¼ pðnÞ , u2k 2k
uðnÞ ¼ 1pðnÞ , 2k þ 1 2k þ 1
d ¼ 1
into (2). Here, Proposition 3.1 indicates that the canonical moments pk are given from the determinant solution of the integrable system (10). It is also to be noted that the discrete Toda equation (11) is equivalent to the qd-algorithm. The qdalgorithm is introduced in Section 1.5 of the book by Dette and Studden (1997). 3.2. Formulation of the PRMm ðlÞ In order to deal PRMm ðlÞ, we introduce a sequence of linear combination of ck and introduce its Hankel determinants. For given moments ck , k ¼ 0; 1,2, . . . and an arbitrary real number l, we define cs,t by the recurrence formula: k cskþ 1,t ¼ cks,tþ 1 ,
þ1 cs,t ¼ cs,t lcs,t , k kþ1 k
c0;0 ¼ ck , k
ð12Þ
where s,t,k are nonnegative integers. We note that cs,t are well-defined. The equations (12) denote the Christoffel k transformation in terms of moments. The Christoffel transformation is described in Chihara (1978), Chapter I. 7. Let Hs,t k be Hankel determinants defined by using (12) as k1
s,t Hs,t k ¼ 9ci þ j 9i,j ¼ 0 ,
where k ¼ 1; 2, . . . , n ¼ 0; 1,2, . . . . We here set that H0s,t ¼ 1, and that Hks,t ¼ 0 if a matrix size k is negative. We can easily find ðnÞ ðnÞ out that Hn,0 k ¼ H k by comparing with the notations (7) of H k . Let us write a linear regression model as 0 1 f 0 ðxÞ B C B f 1 ðxÞ C T C þ e, Y ¼ y f ðxÞ þ e ¼ ðy0 y1 ym1 ÞB B C ^ @ A f m1 ðxÞ where f ðxÞ ¼ ðf 0 ðxÞ,f 1 ðxÞ, . . . ,f m1 ðxÞÞT is a known vector of basis functions. The polynomial regression models PRMm ðlÞ through a fixed point do not correspond to a linear regression model uniquely. However, we can define D-optimal designs for PRMm ðlÞ, since D-optimal designs depend only on the linear space spanned by the basis functions. (See Dette and Studden, 1997, Theorem 5.5.1.) Now, we give a definition of D-optimal design for PRMm ðlÞ. Let the original polynomial regression model be described as Y¼
m X k¼0
yk xk þ e:
ð13Þ
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
939
In the case of PRMm ðlÞ, since we have EðY9x ¼ lÞ, substitute
y0 ¼ EðY9x ¼ lÞy1 ly2 l2 ym lm into (13) to obtain YEðY9x ¼ lÞ ¼
m X
yk ðxk lk Þ þ e:
ð14Þ
k¼1
When we obtain yk by observation, we can calculate the value of yk EðY9x ¼ lÞ easily. Hence, by rewriting YEðY9x ¼ lÞ as Y, we can consider that the model (14) is essentially equivalent to the linear regression model: Y¼
m X
yk ðxk lk Þ þ e:
k¼1
Therefore, PRMm ðlÞ have the vector of basis functions: 2
m
gðxÞ ¼ ðxl,x2 l , . . . ,xm l ÞT : Let A be the following 0 1 0 B l 1 B B A¼B B 0 l B ^ @ ^ 0
0
m m nonsingular matrix 1 0 0 0 0C C C 1 & ^C C: C & & 0A l 1
Then, we can replace g(x) with f ðxÞ ¼ AgðxÞ ¼ ðxlÞð1,x, . . . ,xm1 ÞT : Since the Fisher information matrix is written as
R1 0
f ðxÞf ðxÞT dmðxÞ, we obtain the following lemma.
Lemma 3.2. The D-optimal design m 2 P ½0;1 for PRMm ðlÞ is described as the optimal solution of the optimization problem maximize
H0;2 m ðmÞ,
subjectto
m 2 P ½0;1 :
ð15Þ
Note that the D-optimal design for PRMm ðlÞ is equivalent to the D-optimal design for an (m1)th degree weighted polynomial regression model with weight wðxÞ ¼ ðxlÞ2 , or for an (m1)th degree heteroscedastic polynomial regression model with variance function ðxlÞ2 . The case of l r0 or l Z1 has already been studied. See Dette and Trampisch (2010). 3.3. Construction of the D-optimal designs for PRMm ðlÞ s,t
In order to obtain the D-optimal design for PRMm ðlÞ, we define values corresponding to Hankel determinant Hs,t k . Let zk be defined by
zs,t 2m¼
s þ 1,t Hs,t m þ 1 H m1
Hsmþ 1,t Hs,t m
,
zs,t 2 mþ1 ¼
s,t Hsmþþ1,t 1 Hm s þ 1,t Hs,t m þ 1 Hm
,
m Z 0,
ð16Þ
,
m Z 0,
ð17Þ
s,t
and let xk be defined by
xs,t 2m¼
s,t þ 1 Hs,t m þ 1 H m1 þ 1 s,t Hs,t Hm m
,
x2s,tm þ 1 ¼
þ 1 s,t Hs,t m þ 1 Hm s,t þ 1 Hs,t m þ 1 Hm
n,0
ðnÞ
s,t
s,t
where s,t is nonnegative integers. We can easily show that zk ¼ zk . There are some useful identities among zk and xk which are called the nonautonomous discrete Toda (nd-Toda) equation: s,t zs2þm1,t þ zs2þm1,tþ 1 ¼ z2s,tm þ 1 þ zs,t zs2þm1,tþ 1 zs2þm1,tþ 2 ¼ zs,t 2 m þ 2, 2 m þ 2 z2 m þ 3 ,
ð18aÞ
þ1 s,t þ 1 s,t s,t þ1 s,t þ 1 s,t s,t xs,t xs,t 2 m þ x2 m þ 1 ¼ x2 m þ 1 þ x2 m þ 2 , 2 m þ 1 x2 m þ 2 ¼ x2 m þ 2 x2 m þ 3 ,
ð18bÞ
þ1 s,t þ 1 s,t s,t þ1 s,t þ 1 s,t s,t zs,t zs,t 2 m þ z2 m þ 1 ¼ x2 m þ 1 þ x2 m þ 2 þ l, 2 m þ 1 z2 m þ 2 ¼ x2 m þ 2 x2 m þ 3 ,
ð18cÞ
s,t s,t xs2þm1,t þ xs2þm1,tþ 1 þ l ¼ zs,t xs2þm1,tþ 1 xs2þm1,tþ 2 ¼ zs,t 2 m þ 1 þ z2 m þ 2 , 2 m þ 2 z2 m þ 3 ,
ð18dÞ
where m,s,t are nonnegative integers. Note that (18) reduces to (11) by substituting l ¼ 0. See Kajiwara and Mukaihira (2005) for more details about the nd-Toda equation and its determinant solutions.
940
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
Let us obtain an expression of the objective function H0;2 m ðmÞ in terms of canonical moments pk by the following three s,t s,t steps. First, we can obtain an expression of the objective function H0;2 m in terms of zk or of xk by the following theorem. s,t
s,t
Theorem 3.3. The Hankel determinants Hs,t m have the following expressions in terms of zk or of xk : s,t m Hs,t m ¼ ðc0 Þ
m1 Y
s,t
s,t
ðz2k1 z2k Þmk ¼ ðc0s,t Þm
k¼1
m1 Y
s,t
s,t
ðx2k1 x2k Þmk :
ð19Þ
k¼1
The proof is given in Appendix. From Theorem 3.3, we find out that the objective function is described as m1 Y
0;2 m H0;2 m ¼ ðc0 Þ
0;2
0;2
ðx2k1 x2k Þmk :
k¼1 0;0
Next step, let us obtain an expression of the objective function in terms of zk . Note that the canonical moments pk are 0;0
s,t
s,t
similar to zk . To this aim, we can use the nd-Toda equation (18). The relationship among the variables zk and xk is often illustrated as in Figs. 1 and 2. We can find out that Eqs. (18a) and (18b) describe mutual relations among four vertices in 0;2
each rhombus appeared in Figs. 1 and 2, respectively. Figs. 1 and 2 imply that the variables xk of
z0;0 k
via
x0;1 k ,
x0;0 k
and
xk0,1 .
Here, we consider that the variables
x0,1 k
can be expressed in terms
defined by (17) with some constant c1 . Naturally,
the value of objective function H0;2 m does not depend on c1 . s,t s,t The remaining work in this step is to express c0;2 0 in terms of zk and of xk . We can use the
zst 1 ¼
cs0þ 1t , cst 0
xs,t 1 ¼
þ1 cs,t 0
ð20Þ
cs,t 0 s,t
s,t
that can be obtained by (16) and (17). From Eq. (20), we can express c0;2 0 in terms of zk and of xk , namely, 0;0 0;1
c0;2 0 ¼ x1 x1 : In the last step, let us obtain an expression of the objective function H0;2 m in terms of the canonical moments pk. This is 0;0 the easiest step. We can represent zk ¼ zk by pk by using (8).
s,t
Fig. 1. The relationship among zk .
s,t
Fig. 2. The relationship among xk .
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
941
By putting it all together, the algorithm for calculating D-optimal designs for PRMm ðlÞ is described as follows. The algorithm for calculating D-optimal designs for PRMm ðlÞ Step Step Step Step
1. 2. 3. 4.
s,t
s,t
By using Theorem 3.3, describe the objective function H0;2 m in terms of zk and xk . 0;0 By using the nd-Toda (18), describe the objective function H0;2 m in terms of zk . 0;2 By using (8), describe the objective function Hm in terms of canonical moments pk. Find canonical moments which maximize the objective function H0;2 m .
The D-optimal designs can be obtained by using Theorem 3.4.1 in the book by Dette and Studden (1997) after canonical moments are given. Here, if N is less than 2m, then the canonical moments: pk ,
k ¼ N þ1,N þ 2, . . . ,2m
ð21Þ
and
zk , k ¼ N þ 1,N þ 2, . . . ,2m
ð22Þ
cannot be defined from (5) and (8), respectively. However, it is no matter, because the objective function H0;2 m does not depend on (21) and (22). There are some notes about above algorithm. An expression of the objective function can be obtained by using a computer algebra system, and we can find canonical moments which maximize the objective function numerically for fixed l. There are the relations: þ1 s,t þ 1 s,t s,t þ1 s,t þ 1 s,t s,t zs,t zs,t m Z0, 2 m þ z2 m þ 1 ¼ x2 m þ x2 m þ 1 þ l, 2 m þ 1 z2 m þ 2 ¼ x2 m þ 1 x2 m þ 2 ,
ð23Þ
that are the variant of the nd-Toda equation (18). These relations (23) save the cost of the calculation in the second step.
4. Examples Let us calculate D-optimal designs for PRM2 ðlÞ. In this case, the objective function is H0;2 2 . Perform first three steps in proposed algorithm, we can obtain H0;2 2 as a function of p1 ,p2 ,p3 ,p4 and l by using a computer algebra system. In this case, we can obtain D-optimal design analytically by using a computer algebra system. The D-optimal design for PRM2 ðlÞ is as follows: pp ffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffi ffi . 21Þ 2 r l r1=2: (1) Case l r 1 or ð1 The canonical moments pk of the D-optimal design mn are p1 ¼ 1=2,
p2 ¼ 1,
and the D-optimal design mn has the two support points such as Probmn ð0Þ ¼ 1=2,
Probmn ð1Þ ¼ 1=2,
where Probm ðxÞ denote the probability mass function. (2) Case 1r l r1=6: The canonical moments pk of the D-optimal design mn are p1 ¼
3þl , 4
p2 ¼
1l , 3þl
p3 ¼ 1,
and the D-optimal design mn has the two support points such as Probmn ðð1 þ lÞ=2Þ ¼ 1=2,
Probmn ð1Þ ¼ 1=2:
pffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi ffi . 21Þ 2: (3) Case 1=6r l r ð1 The canonical moments pk of the D-optimal design mn are p1 ¼
lð4l3 4l þ 1Þ , 2 ð2l1Þð8l 8l þ 1Þ 4
p2 ¼
3
2
4ðl1Þð20l 40l þ30l 10l þ1Þl 2
3
3
ð1 þ 8l12l þ 4l Þð14l þ 4l Þ
,
942
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
1 support
0.8 0.6 0.4 0.2 0 -1
-0.5
0 lambda
0.5
1
Fig. 3. A relationship between D-optimal designs and l ð1:2 r l r 1:2Þ.
1 support
0.8 0.6 0.4 0.2 0 0.15
0.155
0.16
0.165 lambda
0.17
0.175
0.18
Fig. 4. A relationship between D-optimal designs and l ð0:15 r l r 0:18Þ. 2
p3 ¼
2
3
l ð6l5Þð1 þ 8l12l þ 4l Þ 4
3
2
ð2l1Þð20l 40l þ 30l 10l þ 1Þ
,
p4 ¼ 1, and the D-optimal design mn has the three support points such as Probmn ð0Þ ¼
Probmn
4ð6l1Þðl1Þ4 2
ð8l 8l þ1Þð4l3Þð2l1Þ
,
4 3 2 lð4l3Þ 8l 16l þ 16l 8l þ 1 ¼ , 2 2l1 ð8l 8l þ 1Þð4l1Þð4l3Þ
Probmn ð1Þ ¼
4ð6l5Þl
4
2
ð2l1Þð4l1Þð8l 8l þ 1Þ
:
(4) Case l 41=2: The D-optimal design mn is the measure corresponding to the D-optimal design for PRMs ð1lÞ under the transformation y ¼ 1x. Figs. 3 and 4 show the relationship between D-optimal designs for PRMm ðlÞ and l. The horizontal axis indicates and the vertical axis indicates indicate the value of l and support points of the D-optimal designs for PRMm ðlÞ. The thickness of the lines, here, corresponds to the respective weights. 5. Conclusions In this paper, we first indicate a relationship between canonical moments and discrete integrable systems. By using this relationship, we propose a new algorithm for calculating an expression of the Hankel determinant H0;2 in terms of the k canonical moments pk, where H0;2 is the objective function of the optimization problem (15). This algorithm enable us to k calculate D-optimal designs for polynomial regression model through a fixed point. In this paper, we use the nd-Toda equation that is a kind of discrete integrable systems. Here, discrete integrable systems describe rhombus diagrams such as Figs. 1 and 2. Therefore they help us to calculate D-optimal designs algorithmically. There are many types of integrable systems that have useful properties such as determinant solutions. We wish that the relationship between canonical moments and discrete integrable systems found in this paper will contribute to advance statistical experiments further.
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
943
Appendix A. Proof of Theorem 3.3 From (16), we obtain s,t zs,t 2 m z2 m1
s þ 1,t s þ 1,t s,t Hs,t Hm1 m þ 1 H m1 H m
¼
Hsmþ 1,t Hs,t m
s þ 1,t Hs,t m H m1
¼
Hs,t mþ1 Hs,t m
,
Hs,t m : Hs,t m1
ð24Þ
Eq. (24) leads Hs,t mþ1 Hs,t m
s,t
s,t
¼ z2 m z2 m1
m Y Hs,t m ¼ cs,t ðz2k1 z2k Þ: 0 s,t Hm1 k¼1
This implies that the first equation of (19): s,t m Hs,t m ¼ ðc0 Þ
m1 Y
s,t
s,t
ðz2k1 z2k Þmk :
k¼1
From (16) and (17), s,t xs,t 2 m x2 m1 ¼
s,t Hs,t m þ 1 H m1 s,t Hs,t m Hm
s,t
s,t
¼ z2 m z2 m1
that implies m ðcs,t 0 Þ
m1 Y k¼1
s,t
s,t
m ðz2k1 z2k Þmk ¼ ðcs,t 0 Þ
m1 Y
s,t
s,t
ðx2k1 x2k Þmk :
k¼1
References Arnold, V.I., 1978. Mathematical Methods of Classical Mechanics. Springer-Verlag, New York, Heidelberg. Biedermann, S., Dette, H., Hoffmann, P., 2009. Constrained optimal discrimination designs for Fourier regression models. Annals of the Institute of Statistical Mathematics 61, 143–157. Chihara, T.S., 1978. An Introduction to Orthogonal Polynomials. Gordon and Breach, New York. Dette, H., 1992. Optimal designs for a class of polynomial regression of odd or even degree. Annals of Statistics 20, 238–259. Dette, H., Franke, T., 2000. Constrained D- and D1-optimal designs for polynomial regression. Annals of Statistics 28 (6), 1702–1727. Dette, H., Franke, T., 2001. Robust designs for polynomial regression by maximizing a minimum of D- and D1-efficiencies. Annals of Statistics 29 (4), 1024–1049. Dette, H., Haines, L.M., Imhof, L.A., 2007. Maximin and Bayesian optimal designs for regression models. Statistica Sinica 17, 463–480. Dette, H., Studden, W., 1997. The Theory of Canonical Moments with Applications in Statistics, Probability and Analysis. Wiley, New York. Dette, H., Trampisch, M., 2010. A general approach to D-optimal designs for weighted univariate polynomial regression models. Journal of Korean Statistical Society 39 (1), 1–26. Fang, Z., 2002. D-optimal designs for polynomial regression models through origin. Statistics & Probability Letters 57, 343–351. Fang, Z., 2006. Some robust designs for polynomial regression models. Canadian Journal of Statistics 34, 623–638. Fang, Z., Wiens, D.P., 2003. Robust regression designs for approximate polynomial models. Journal of Statistical Planning and Inference 117 (2), 305–321. Hirota, R., 1997. Conserved quantities of random-time Toda equation. Journal of the Physical Society of Japan 66, 283–284. Huang, M.N.L., Chang, F.C., Wong, W.K., 1995. D-optimal designs for polynomial regression without intercept. Statistica Sinica 5, 441–458. Kajiwara, K., Mukaihira, A., 2005. Soliton solutions for the non-autonomous discrete-time Toda lattice equation. Journal of Physics A: Mathematical and General 38 (28), 6363–6370. Lau, T.S., 1988. D-optimal designs on the unit q-ball. Journal of Statistical Planning and Inference 19, 299–315. Lau, T.S., Studden, W.J., 1985. Optimal designs for trigonometric and polynomial regression using canonical moments. Annals of Statistics 13, 385–394. Lim, Y.B., Studden, W.J., 1988. Efficient Ds-optimal designs for multivariate polynomial regression on the q-cube. Annals of Statistics 16, 1225–1240. Nakamura, Y., 1988. A tau-function of the finite nonperiodic Toda lattice. Physical Letters A 195 (5–6), 346–350. Pukelsheim, F., 1993. Optimal Design of Experiments. Wiley, New York. Skibinsky, M., 1967. The range of the (nþ 1)th moment for distributions on [0, 1]. Journal of Applied Probability 4 (3), 543–552. Skibinsky, M., 1968. Extreme nth moments for distributions on ½0; 1 and the inverse of a moment space map. Journal of Applied Probability 5 (3), 693–701. Skibinsky, M., 1969. Some striking properties of binomial and beta moments. Annals of Mathematical Statistics 40 (5), 1753–1764. Spiridonov, V., Zhedanov, A., 1995. Discrete Darboux transformations, the discrete-time Toda lattice, and the Askey–Wilson polynomials. Methods and Applications of Analysis 2, 369–398. Spiridonov, V., Zhedanov, A., 1997. Discrete-time Volterra chain and classical orthogonal polynomials. Journal of Physics A: Mathematical and General 30, 8727–8737. Studden, W.J., 1980. Ds-optimal designs for polynomial regression using continued fractions. Annals of Statistics 8 (5), 1132–1141. Studden, W.J., 1982a. Optimal designs for weighted polynomial regression using canonical moments. Statistical Decision Theory and Related Topics III 2, 335–350. Studden, W.J., 1982b. Some robust-type D-optimal design in polynomial regression. Journal of American Statistical Association 77, 916–921. Studden, W.J., 1989. Note on some Fp optimal design for polynomial regression. Annals of Statistics 17, 618–623. Symes, W.W., 1982. The QR algorithm and scattering for the finite nonperiodic Toda lattice. Physica D 4 (2), 275–280. Tsujimoto, S., Nakamura, Y., Iwasaki, M., 2001. The discrete Lotka–Volterra system computes singular values. Inverse Problems 17 (1), 53–58. Zen, M.M., Tsai, M.H., 2004. Criterion-robust optimal designs for model discrimination and parameter estimation in Fourier regression models. Journal of Statistical Planning and Inference 124 (2), 475–487.
Contents lists available at SciVerse ScienceDirect
Journal of Statistical Planning and Inference journal homepage: www.elsevier.com/locate/jspi
An algorithm for calculating D-optimal designs for polynomial regression through a fixed point Hiroto Sekido Department of Applied Mathematics and Physics, Graduate School of Informatics, Kyoto University, Kyoto 606-8501, Japan
a r t i c l e i n f o
abstract
Article history: Received 11 January 2010 Received in revised form 10 September 2011 Accepted 31 October 2011 Available online 10 November 2011
Optimal designs are required to make efficient statistical experiments. By using canonical moments, in 1980, Studden found Ds-optimal designs for polynomial regression models. On the other hand, integrable systems are dynamical systems whose solutions can be written down concretely. In this paper, polynomial regression models through a fixed point are discussed. In order to calculate D-optimal designs for these models, a useful relationship between canonical moments and discrete integrable systems is introduced. By using canonical moments and discrete integrable systems, a new algorithm for calculating D-optimal designs for these models is proposed. & 2011 Elsevier B.V. All rights reserved.
Keywords: D-optimal design Canonical moment Polynomial regression model Integrable system Nonautonomous discrete Toda equation
1. Introduction For a given probability measure on the Borel sets of a finite interval, we can define canonical moments in addition to ordinary moments. Canonical moments are defined as normalized moments in some sense. For example, canonical moments pk satisfy 0 rpk r 1. Other properties of canonical moments were studied by Skibinsky (1967, 1968, 1969). Applications of canonical moments are explained in Dette and Studden (1997). One of the applications is a design of experiments. In this paper, we consider the D-optimal designs which corresponds to an objective function of determinant form. D-optimal designs and Ds-optimal designs for polynomial regression models were studied by Studden (1980). In the calculation by Studden (1980), an important point is to use canonical moments instead of ordinary moments, and D-optimal designs can be identified by its canonical moments. As we see in Section 2, the objective function is written down in term of canonical moments. Besides this, D-optimal designs for various models have been calculated by Biedermann et al. (2009), Dette (1992), Dette and Franke (2000, 2001), Dette et al. (2007), Fang (2002, 2006), Fang and Wiens (2003), Huang et al. (1995), Lau (1988), Lau and Studden (1985), Lim and Studden (1988), Studden (1982b, 1982a, 1989), Zen and Tsai (2004). For example, D-optimal designs for weighted polynomial regression with a weight function xa ð1xÞb were found by Studden (1982a) by using canonical moments. However, in most cases, an explicit form of the D-optimal design is unknown. We here consider polynomial regression models with some prior information. For example, D-optimal designs for this kind of models has been studied by Dette (1992), Fang (2002), Huang et al. (1995). Dette (1992) calculated D-optimal designs for polynomial regression models with only odd (or even) degree terms. Huang et al. (1995) considered polynomial regression without intercept. Fang (2002) considered polynomial regression models through origin.
E-mail address: [email protected] 0378-3758/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2011.10.014
936
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
On the other hand, the term integrable systems is used for nonlinear dynamical systems whose solutions can be written down concretely. For Hamiltonian systems with finite degree of freedom, their integrability is defined in the Liouville– Arnold theorem (Arnold, 1978). However, even now, there is no mathematical definition of integrability for nonlinear systems with infinite degree of freedom. This is why we call an explicitly solvable nonlinear system as an integrable system. For example, the Lotka–Volterra equation: dU k ¼ U k ðU k þ 1 U k1 Þ, dt
ð1Þ
known as a model equation for a food chain is one of integrable systems. Discrete integrable systems are discrete analogues of integrable systems. For example, the discrete Lotka–Volterra equation (Spiridonov and Zhedanov, 1997): þ 1Þ þ 1Þ ð1þ duðn Þ ¼ uðnÞ ð1 þ dukðnÞþ 1 Þ uðn k k1 k
ð2Þ ukðnÞ
indicates the value of uk at the time nd, is one of discrete integrable systems, where the integer n is a discrete time and d is a discrete time step. Note that the discrete Lotka–Volterra equation (2) reduces in the continuous limit d-0 with nd ¼ T to the Lotka–Volterra equation (1). Both (1) and (2) have explicit solutions of determinant form. Here, we approach to a relationship between canonical moments and determinant solutions of discrete integrable systems. This is a new viewpoint in optimal designs. In this paper, we propose an algorithm for calculating D-optimal design for polynomial regression through a given fixed point in x ¼ l, where l is an arbitrary real number. We here consider the design space is X ¼ ½0; 1. When l ¼ 0 and l ¼ 1, the corresponding D-optimal designs are equivalent to D-optimal designs for weighted polynomial regression with a weight function x2 and ð1xÞ2 , respectively. These are special cases of D-optimal designs calculated in Studden (1982a). Fang (2002) proposed the D-optimal designs for polynomial regression through origin with the design space X ¼ ½1; 1. It contains the case of l ¼ 1=2. While in this paper, we consider polynomial regression through an arbitrary fixed point. See Huang et al. (1995) for the motivation in terms of regression models. The resulting algorithm for calculating D-optimal designs can be regarded as a new application of discrete integrable systems. Integrable systems have been applied to numerical analysis. Typical examples are matrix eigenvalue algorithms (Nakamura, 1988; Symes, 1982), and algorithms for computing matrix singular values (Tsujimoto et al., 2001) in terms of the discrete Lotka–Volterra equation (2). Section 3 shows that canonical moments can be rephrased as determinant solutions of (2). This is the starting point of our approach. Therefore, discrete integrable systems, which are generalizations of (2), help us to write down the objective function in terms of canonical moments for polynomial regression through a fixed point. In Section 2, we introduce the notion of D-optimal designs and give a short review on the calculation of D-optimal designs by Studden (1980). Section 3 indicates a relationship between canonical moments and discrete integrable systems at first. Then, an algorithm for calculating D-optimal designs for polynomial regression through a fixed point is proposed. In Section 4, we calculate D-optimal designs for a second degree polynomial regression through a fixed point in terms of the nonautonomous discrete Toda equation, an integrable system. Section 5 is devoted to conclusions. 2. Polynomial regression models and its D-optimal designs: a preliminary In this section, we give a brief review on the D-optimal designs for polynomial regression models. Let us call polynomial regression models PRMs, for short. After that, we introduce canonical moments and give a review on the calculation of D-optimal designs by Studden (1980). Consider common (m1)th degree PRMs: Y ¼ ym1 xm1 þ ym2 xm2 þ þ y0 þ e, where yk is an unknown parameter, and e is a random error term. Let P ½0;1 denote the set of all probability measures on the Borel sets of the interval ½0; 1. For given m 2 P ½0;1 , let ck be the kth moment given by Z 1 xk dmðxÞ, k ¼ 0; 1,2, . . . : ck ¼ 0
The D-optimal designs are defined as probability measures in P ½0;1 which maximize the determinant of the Fisher information matrix, that is, D-optimal designs m 2 P ½0;1 for PRMs are defined as the optimal solution of the optimization problem maximize
Hm ðmÞ,
subject to
m 2 P ½0;1 ,
ð3Þ
where let Hk ðmÞ be the Hankel determinants of moments as m1
Hm ðmÞ ¼ 9ci þ j 9i,j ¼ 0 : See Pukelsheim (1993) for more details about optimal designs.
ð4Þ
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
937
The objective function Hm of (3) is explicitly expressed in terms of the moments ck. However, in terms of moments ck, the set P ½0;1 of feasible solutions of (3) is expressed in a very complicated form. Therefore, we use canonical moments instead of the ordinary moments ck. This method was proposed by Studden (1980) at first, and is explained in the book by Dette and Studden (1997). Now, we define canonical moments. For a given probability measure m 2 P ½0;1 , let ckþ denote the maximum of the k-th moment over the set of all measure having moments c0 ,c1 , . . . ,ck1 . Similarly let c k denote the corresponding minimum. Canonical moments are defined by pk ¼
ck c k , ckþ c k
k ¼ 1; 2, . . . ,N,
ð5Þ
þ where N is the minimum of j which satisfy cjþþ 1 ¼ c j þ 1 . If c j 4 cj for an arbitrary positive integer j, let N be 1. Canonical moments have the property that ( pk 2 ð0; 1Þ, k ¼ 1; 2, . . . ,N1,
pN 2 f0; 1g:
ð6Þ
Conversely, for a given arbitrary sequence fpk gN k ¼ 1 which satisfies (6), there is a unique m 2 P ½0;1 which has the canonical moments fpk gN k ¼ 1. It is to be noted that canonical moments have a Hankel determinant expression. For introduction, we generalize the ðnÞ definition (4) of Hankel determinants Hk ðmÞ. Let HðnÞ ,H k be the Hankel determinants defined by k k1
HkðnÞ ¼ 9ci þ j þ n 9i,j ¼ 0 ,
k1
ðnÞ
H k ¼ 9ci þ j þ n ci þ j þ n þ 1 9i,j ¼ 0 ,
ð7Þ
ðnÞ where k ¼ 1; 2, . . . , n ¼ 0; 1,2, . . . . We here set that HðnÞ 0 ¼ 1, and that H k ¼ 0 if a matrix size k is negative. Obviously, we ð0Þ have Hk ¼ Hk . The canonical moments pk have the following Hankel determinant expression: ð1Þ
ð0Þ
p2k1 ¼
Hð1Þ H k1 k ð1Þ Hkð0Þ H k1
,
p2k ¼
Hð0Þ H k þ 1 k1 ð0Þ Hð1Þ Hk k
,
k ¼ 1; 2, . . . :
We introduce zk by the transformation of canonical moments pk
z0 ¼ 0, z1 ¼ p1 , zk ¼ ð1pk1 Þpk , k ¼ 2; 3, . . . ,N:
ð8Þ
Then zk also have the Hankel determinant expression:
z2k1 ¼
Hð1Þ Hð0Þ k k1 Hð0Þ Hð1Þ k k1
,
z2k ¼
Hkð0Þþ 1 Hð1Þ k1 Hð1Þ Hð0Þ k k
,
k ¼ 1; 2, . . . :
From this expression, Hankel determinants Hm can be expressed in a product of zk or of canonical moments pk, Hm ¼
m1 Y
ðz2k1 z2k Þmk
k¼1
0 ¼@
m1 Y
1 A ð1p2j Þmj1 pmj 2j
j¼1
m1 Y
ðð1p2j1 Þp2j1 Þmj ,
ð9Þ
j¼1
where 2 m2 r N. The expression (9) enable us to find the canonical moments pk which maximize the objective function Hm of (3) by a simple calculation. Consequently, we obtain the canonical moments pk of D-optimal designs for PRMs p2j1 ¼
1 , 2
p2j ¼
mj , 2ðmjÞ1
j ¼ 1; 2, . . . ,m1:
This approach for D-optimal designs is described in the book by Dette and Studden (1997). 3. D-optimal designs for the polynomial regression models through a fixed point In this section, we consider the mth degree polynomial regression models through a known point. In this case, we know a value of Eðy9x ¼ lÞ as a prior information, where l is an arbitrary given real number. Let us denote an mth degree polynomial regression model through a fixed point in x ¼ l as PRMm ðlÞ. In Section 3.2, PRMm ðlÞ and the corresponding D-optimal designs will be formulated as D-optimal design for a linear regression model. By Studden (1980), the objective function Hm ðmÞ of (3) is expressed in terms of canonical moments. In Section 3.3, we propose a method to express the objective function corresponding to PRMm ðlÞ in terms of canonical moments.
938
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
3.1. Canonical moments and discrete integrable systems We indicate a relationship between canonical moments and discrete integrable systems, before considering a formulation of the PRMm ðlÞ. In order to obtain a D-optimal design for PRMm ðlÞ, we can get clues from discrete integrable systems. Hankel determinants are related to both integrable systems and D-optimal designs, furthermore (some generalized) canonical moments pk and the variables zk satisfy a kind of discrete integrable systems. be defined by Proposition 3.1. Let pðnÞ k ðn þ 1Þ
ðnÞ ¼ p2k
HðnÞ H k þ 1 k1
ðnÞ
þ 1Þ Hðn Hk k
ðnÞ
,
pðnÞ ¼ 2k þ 1
þ 1Þ Hðn Hk kþ1
ðn þ 1Þ
HðnÞ H kþ1 k
,
k,n ¼ 0; 1,2, . . . :
satisfy the analogue of the discrete Lotka–Volterra equation: Then pðnÞ k þ 1Þ þ 1Þ ðnÞ ð1pðn Þð1 þ pðn Þ ¼ ð1pðnÞ Þð1þ p2k Þ, 2k þ 1 2k 2k þ 1 þ2 þ 1Þ ðn þ 1Þ ðnÞ p ¼ p2k pðnÞ , pðn 2k þ 2 2k þ 1 þ 2 2k þ 3
Let
zðnÞ k
k,n ¼ 0; 1,2, . . . :
ð10Þ
be defined by
zðnÞ 2k ¼
þ 1Þ HkðnÞþ 1 Hðn k1 þ 1Þ ðnÞ Hðn Hk k
,
ðnÞ z2k þ1 ¼
þ 1Þ ðnÞ Hðn Hk kþ1
HðnÞ Hðn þ 1Þ kþ1 k
,
k,n ¼ 0; 1,2, . . . :
ðnÞ
Then zk satisfy the discrete Toda equation (Hirota, 1997; Spiridonov and Zhedanov, 1995) þ 1Þ ðn þ 1Þ ðnÞ ðnÞ ðn þ 1Þ ðn þ 1Þ ðnÞ ðnÞ zðn þ z2k þ 1 ¼ z2k þ 1 þ z2k þ 2 , z2k þ 1 z2k þ 2 ¼ z2k þ 2 z2k þ 3 , 2k
ð11Þ
where k,n ¼ 0; 1, . . . . ð0Þ
¼ pk , and that zk ¼ zk . This proposition can be prove by using the determinant solution of integrable It is clear that pð0Þ k systems and Sylvester’s determinant identity. Note that (10) is obtained by substituting ðnÞ ¼ pðnÞ , u2k 2k
uðnÞ ¼ 1pðnÞ , 2k þ 1 2k þ 1
d ¼ 1
into (2). Here, Proposition 3.1 indicates that the canonical moments pk are given from the determinant solution of the integrable system (10). It is also to be noted that the discrete Toda equation (11) is equivalent to the qd-algorithm. The qdalgorithm is introduced in Section 1.5 of the book by Dette and Studden (1997). 3.2. Formulation of the PRMm ðlÞ In order to deal PRMm ðlÞ, we introduce a sequence of linear combination of ck and introduce its Hankel determinants. For given moments ck , k ¼ 0; 1,2, . . . and an arbitrary real number l, we define cs,t by the recurrence formula: k cskþ 1,t ¼ cks,tþ 1 ,
þ1 cs,t ¼ cs,t lcs,t , k kþ1 k
c0;0 ¼ ck , k
ð12Þ
where s,t,k are nonnegative integers. We note that cs,t are well-defined. The equations (12) denote the Christoffel k transformation in terms of moments. The Christoffel transformation is described in Chihara (1978), Chapter I. 7. Let Hs,t k be Hankel determinants defined by using (12) as k1
s,t Hs,t k ¼ 9ci þ j 9i,j ¼ 0 ,
where k ¼ 1; 2, . . . , n ¼ 0; 1,2, . . . . We here set that H0s,t ¼ 1, and that Hks,t ¼ 0 if a matrix size k is negative. We can easily find ðnÞ ðnÞ out that Hn,0 k ¼ H k by comparing with the notations (7) of H k . Let us write a linear regression model as 0 1 f 0 ðxÞ B C B f 1 ðxÞ C T C þ e, Y ¼ y f ðxÞ þ e ¼ ðy0 y1 ym1 ÞB B C ^ @ A f m1 ðxÞ where f ðxÞ ¼ ðf 0 ðxÞ,f 1 ðxÞ, . . . ,f m1 ðxÞÞT is a known vector of basis functions. The polynomial regression models PRMm ðlÞ through a fixed point do not correspond to a linear regression model uniquely. However, we can define D-optimal designs for PRMm ðlÞ, since D-optimal designs depend only on the linear space spanned by the basis functions. (See Dette and Studden, 1997, Theorem 5.5.1.) Now, we give a definition of D-optimal design for PRMm ðlÞ. Let the original polynomial regression model be described as Y¼
m X k¼0
yk xk þ e:
ð13Þ
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
939
In the case of PRMm ðlÞ, since we have EðY9x ¼ lÞ, substitute
y0 ¼ EðY9x ¼ lÞy1 ly2 l2 ym lm into (13) to obtain YEðY9x ¼ lÞ ¼
m X
yk ðxk lk Þ þ e:
ð14Þ
k¼1
When we obtain yk by observation, we can calculate the value of yk EðY9x ¼ lÞ easily. Hence, by rewriting YEðY9x ¼ lÞ as Y, we can consider that the model (14) is essentially equivalent to the linear regression model: Y¼
m X
yk ðxk lk Þ þ e:
k¼1
Therefore, PRMm ðlÞ have the vector of basis functions: 2
m
gðxÞ ¼ ðxl,x2 l , . . . ,xm l ÞT : Let A be the following 0 1 0 B l 1 B B A¼B B 0 l B ^ @ ^ 0
0
m m nonsingular matrix 1 0 0 0 0C C C 1 & ^C C: C & & 0A l 1
Then, we can replace g(x) with f ðxÞ ¼ AgðxÞ ¼ ðxlÞð1,x, . . . ,xm1 ÞT : Since the Fisher information matrix is written as
R1 0
f ðxÞf ðxÞT dmðxÞ, we obtain the following lemma.
Lemma 3.2. The D-optimal design m 2 P ½0;1 for PRMm ðlÞ is described as the optimal solution of the optimization problem maximize
H0;2 m ðmÞ,
subjectto
m 2 P ½0;1 :
ð15Þ
Note that the D-optimal design for PRMm ðlÞ is equivalent to the D-optimal design for an (m1)th degree weighted polynomial regression model with weight wðxÞ ¼ ðxlÞ2 , or for an (m1)th degree heteroscedastic polynomial regression model with variance function ðxlÞ2 . The case of l r0 or l Z1 has already been studied. See Dette and Trampisch (2010). 3.3. Construction of the D-optimal designs for PRMm ðlÞ s,t
In order to obtain the D-optimal design for PRMm ðlÞ, we define values corresponding to Hankel determinant Hs,t k . Let zk be defined by
zs,t 2m¼
s þ 1,t Hs,t m þ 1 H m1
Hsmþ 1,t Hs,t m
,
zs,t 2 mþ1 ¼
s,t Hsmþþ1,t 1 Hm s þ 1,t Hs,t m þ 1 Hm
,
m Z 0,
ð16Þ
,
m Z 0,
ð17Þ
s,t
and let xk be defined by
xs,t 2m¼
s,t þ 1 Hs,t m þ 1 H m1 þ 1 s,t Hs,t Hm m
,
x2s,tm þ 1 ¼
þ 1 s,t Hs,t m þ 1 Hm s,t þ 1 Hs,t m þ 1 Hm
n,0
ðnÞ
s,t
s,t
where s,t is nonnegative integers. We can easily show that zk ¼ zk . There are some useful identities among zk and xk which are called the nonautonomous discrete Toda (nd-Toda) equation: s,t zs2þm1,t þ zs2þm1,tþ 1 ¼ z2s,tm þ 1 þ zs,t zs2þm1,tþ 1 zs2þm1,tþ 2 ¼ zs,t 2 m þ 2, 2 m þ 2 z2 m þ 3 ,
ð18aÞ
þ1 s,t þ 1 s,t s,t þ1 s,t þ 1 s,t s,t xs,t xs,t 2 m þ x2 m þ 1 ¼ x2 m þ 1 þ x2 m þ 2 , 2 m þ 1 x2 m þ 2 ¼ x2 m þ 2 x2 m þ 3 ,
ð18bÞ
þ1 s,t þ 1 s,t s,t þ1 s,t þ 1 s,t s,t zs,t zs,t 2 m þ z2 m þ 1 ¼ x2 m þ 1 þ x2 m þ 2 þ l, 2 m þ 1 z2 m þ 2 ¼ x2 m þ 2 x2 m þ 3 ,
ð18cÞ
s,t s,t xs2þm1,t þ xs2þm1,tþ 1 þ l ¼ zs,t xs2þm1,tþ 1 xs2þm1,tþ 2 ¼ zs,t 2 m þ 1 þ z2 m þ 2 , 2 m þ 2 z2 m þ 3 ,
ð18dÞ
where m,s,t are nonnegative integers. Note that (18) reduces to (11) by substituting l ¼ 0. See Kajiwara and Mukaihira (2005) for more details about the nd-Toda equation and its determinant solutions.
940
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
Let us obtain an expression of the objective function H0;2 m ðmÞ in terms of canonical moments pk by the following three s,t s,t steps. First, we can obtain an expression of the objective function H0;2 m in terms of zk or of xk by the following theorem. s,t
s,t
Theorem 3.3. The Hankel determinants Hs,t m have the following expressions in terms of zk or of xk : s,t m Hs,t m ¼ ðc0 Þ
m1 Y
s,t
s,t
ðz2k1 z2k Þmk ¼ ðc0s,t Þm
k¼1
m1 Y
s,t
s,t
ðx2k1 x2k Þmk :
ð19Þ
k¼1
The proof is given in Appendix. From Theorem 3.3, we find out that the objective function is described as m1 Y
0;2 m H0;2 m ¼ ðc0 Þ
0;2
0;2
ðx2k1 x2k Þmk :
k¼1 0;0
Next step, let us obtain an expression of the objective function in terms of zk . Note that the canonical moments pk are 0;0
s,t
s,t
similar to zk . To this aim, we can use the nd-Toda equation (18). The relationship among the variables zk and xk is often illustrated as in Figs. 1 and 2. We can find out that Eqs. (18a) and (18b) describe mutual relations among four vertices in 0;2
each rhombus appeared in Figs. 1 and 2, respectively. Figs. 1 and 2 imply that the variables xk of
z0;0 k
via
x0;1 k ,
x0;0 k
and
xk0,1 .
Here, we consider that the variables
x0,1 k
can be expressed in terms
defined by (17) with some constant c1 . Naturally,
the value of objective function H0;2 m does not depend on c1 . s,t s,t The remaining work in this step is to express c0;2 0 in terms of zk and of xk . We can use the
zst 1 ¼
cs0þ 1t , cst 0
xs,t 1 ¼
þ1 cs,t 0
ð20Þ
cs,t 0 s,t
s,t
that can be obtained by (16) and (17). From Eq. (20), we can express c0;2 0 in terms of zk and of xk , namely, 0;0 0;1
c0;2 0 ¼ x1 x1 : In the last step, let us obtain an expression of the objective function H0;2 m in terms of the canonical moments pk. This is 0;0 the easiest step. We can represent zk ¼ zk by pk by using (8).
s,t
Fig. 1. The relationship among zk .
s,t
Fig. 2. The relationship among xk .
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
941
By putting it all together, the algorithm for calculating D-optimal designs for PRMm ðlÞ is described as follows. The algorithm for calculating D-optimal designs for PRMm ðlÞ Step Step Step Step
1. 2. 3. 4.
s,t
s,t
By using Theorem 3.3, describe the objective function H0;2 m in terms of zk and xk . 0;0 By using the nd-Toda (18), describe the objective function H0;2 m in terms of zk . 0;2 By using (8), describe the objective function Hm in terms of canonical moments pk. Find canonical moments which maximize the objective function H0;2 m .
The D-optimal designs can be obtained by using Theorem 3.4.1 in the book by Dette and Studden (1997) after canonical moments are given. Here, if N is less than 2m, then the canonical moments: pk ,
k ¼ N þ1,N þ 2, . . . ,2m
ð21Þ
and
zk , k ¼ N þ 1,N þ 2, . . . ,2m
ð22Þ
cannot be defined from (5) and (8), respectively. However, it is no matter, because the objective function H0;2 m does not depend on (21) and (22). There are some notes about above algorithm. An expression of the objective function can be obtained by using a computer algebra system, and we can find canonical moments which maximize the objective function numerically for fixed l. There are the relations: þ1 s,t þ 1 s,t s,t þ1 s,t þ 1 s,t s,t zs,t zs,t m Z0, 2 m þ z2 m þ 1 ¼ x2 m þ x2 m þ 1 þ l, 2 m þ 1 z2 m þ 2 ¼ x2 m þ 1 x2 m þ 2 ,
ð23Þ
that are the variant of the nd-Toda equation (18). These relations (23) save the cost of the calculation in the second step.
4. Examples Let us calculate D-optimal designs for PRM2 ðlÞ. In this case, the objective function is H0;2 2 . Perform first three steps in proposed algorithm, we can obtain H0;2 2 as a function of p1 ,p2 ,p3 ,p4 and l by using a computer algebra system. In this case, we can obtain D-optimal design analytically by using a computer algebra system. The D-optimal design for PRM2 ðlÞ is as follows: pp ffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffi ffi . 21Þ 2 r l r1=2: (1) Case l r 1 or ð1 The canonical moments pk of the D-optimal design mn are p1 ¼ 1=2,
p2 ¼ 1,
and the D-optimal design mn has the two support points such as Probmn ð0Þ ¼ 1=2,
Probmn ð1Þ ¼ 1=2,
where Probm ðxÞ denote the probability mass function. (2) Case 1r l r1=6: The canonical moments pk of the D-optimal design mn are p1 ¼
3þl , 4
p2 ¼
1l , 3þl
p3 ¼ 1,
and the D-optimal design mn has the two support points such as Probmn ðð1 þ lÞ=2Þ ¼ 1=2,
Probmn ð1Þ ¼ 1=2:
pffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi ffi . 21Þ 2: (3) Case 1=6r l r ð1 The canonical moments pk of the D-optimal design mn are p1 ¼
lð4l3 4l þ 1Þ , 2 ð2l1Þð8l 8l þ 1Þ 4
p2 ¼
3
2
4ðl1Þð20l 40l þ30l 10l þ1Þl 2
3
3
ð1 þ 8l12l þ 4l Þð14l þ 4l Þ
,
942
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
1 support
0.8 0.6 0.4 0.2 0 -1
-0.5
0 lambda
0.5
1
Fig. 3. A relationship between D-optimal designs and l ð1:2 r l r 1:2Þ.
1 support
0.8 0.6 0.4 0.2 0 0.15
0.155
0.16
0.165 lambda
0.17
0.175
0.18
Fig. 4. A relationship between D-optimal designs and l ð0:15 r l r 0:18Þ. 2
p3 ¼
2
3
l ð6l5Þð1 þ 8l12l þ 4l Þ 4
3
2
ð2l1Þð20l 40l þ 30l 10l þ 1Þ
,
p4 ¼ 1, and the D-optimal design mn has the three support points such as Probmn ð0Þ ¼
Probmn
4ð6l1Þðl1Þ4 2
ð8l 8l þ1Þð4l3Þð2l1Þ
,
4 3 2 lð4l3Þ 8l 16l þ 16l 8l þ 1 ¼ , 2 2l1 ð8l 8l þ 1Þð4l1Þð4l3Þ
Probmn ð1Þ ¼
4ð6l5Þl
4
2
ð2l1Þð4l1Þð8l 8l þ 1Þ
:
(4) Case l 41=2: The D-optimal design mn is the measure corresponding to the D-optimal design for PRMs ð1lÞ under the transformation y ¼ 1x. Figs. 3 and 4 show the relationship between D-optimal designs for PRMm ðlÞ and l. The horizontal axis indicates and the vertical axis indicates indicate the value of l and support points of the D-optimal designs for PRMm ðlÞ. The thickness of the lines, here, corresponds to the respective weights. 5. Conclusions In this paper, we first indicate a relationship between canonical moments and discrete integrable systems. By using this relationship, we propose a new algorithm for calculating an expression of the Hankel determinant H0;2 in terms of the k canonical moments pk, where H0;2 is the objective function of the optimization problem (15). This algorithm enable us to k calculate D-optimal designs for polynomial regression model through a fixed point. In this paper, we use the nd-Toda equation that is a kind of discrete integrable systems. Here, discrete integrable systems describe rhombus diagrams such as Figs. 1 and 2. Therefore they help us to calculate D-optimal designs algorithmically. There are many types of integrable systems that have useful properties such as determinant solutions. We wish that the relationship between canonical moments and discrete integrable systems found in this paper will contribute to advance statistical experiments further.
H. Sekido / Journal of Statistical Planning and Inference 142 (2012) 935–943
943
Appendix A. Proof of Theorem 3.3 From (16), we obtain s,t zs,t 2 m z2 m1
s þ 1,t s þ 1,t s,t Hs,t Hm1 m þ 1 H m1 H m
¼
Hsmþ 1,t Hs,t m
s þ 1,t Hs,t m H m1
¼
Hs,t mþ1 Hs,t m
,
Hs,t m : Hs,t m1
ð24Þ
Eq. (24) leads Hs,t mþ1 Hs,t m
s,t
s,t
¼ z2 m z2 m1
m Y Hs,t m ¼ cs,t ðz2k1 z2k Þ: 0 s,t Hm1 k¼1
This implies that the first equation of (19): s,t m Hs,t m ¼ ðc0 Þ
m1 Y
s,t
s,t
ðz2k1 z2k Þmk :
k¼1
From (16) and (17), s,t xs,t 2 m x2 m1 ¼
s,t Hs,t m þ 1 H m1 s,t Hs,t m Hm
s,t
s,t
¼ z2 m z2 m1
that implies m ðcs,t 0 Þ
m1 Y k¼1
s,t
s,t
m ðz2k1 z2k Þmk ¼ ðcs,t 0 Þ
m1 Y
s,t
s,t
ðx2k1 x2k Þmk :
k¼1
References Arnold, V.I., 1978. Mathematical Methods of Classical Mechanics. Springer-Verlag, New York, Heidelberg. Biedermann, S., Dette, H., Hoffmann, P., 2009. Constrained optimal discrimination designs for Fourier regression models. Annals of the Institute of Statistical Mathematics 61, 143–157. Chihara, T.S., 1978. An Introduction to Orthogonal Polynomials. Gordon and Breach, New York. Dette, H., 1992. Optimal designs for a class of polynomial regression of odd or even degree. Annals of Statistics 20, 238–259. Dette, H., Franke, T., 2000. Constrained D- and D1-optimal designs for polynomial regression. Annals of Statistics 28 (6), 1702–1727. Dette, H., Franke, T., 2001. Robust designs for polynomial regression by maximizing a minimum of D- and D1-efficiencies. Annals of Statistics 29 (4), 1024–1049. Dette, H., Haines, L.M., Imhof, L.A., 2007. Maximin and Bayesian optimal designs for regression models. Statistica Sinica 17, 463–480. Dette, H., Studden, W., 1997. The Theory of Canonical Moments with Applications in Statistics, Probability and Analysis. Wiley, New York. Dette, H., Trampisch, M., 2010. A general approach to D-optimal designs for weighted univariate polynomial regression models. Journal of Korean Statistical Society 39 (1), 1–26. Fang, Z., 2002. D-optimal designs for polynomial regression models through origin. Statistics & Probability Letters 57, 343–351. Fang, Z., 2006. Some robust designs for polynomial regression models. Canadian Journal of Statistics 34, 623–638. Fang, Z., Wiens, D.P., 2003. Robust regression designs for approximate polynomial models. Journal of Statistical Planning and Inference 117 (2), 305–321. Hirota, R., 1997. Conserved quantities of random-time Toda equation. Journal of the Physical Society of Japan 66, 283–284. Huang, M.N.L., Chang, F.C., Wong, W.K., 1995. D-optimal designs for polynomial regression without intercept. Statistica Sinica 5, 441–458. Kajiwara, K., Mukaihira, A., 2005. Soliton solutions for the non-autonomous discrete-time Toda lattice equation. Journal of Physics A: Mathematical and General 38 (28), 6363–6370. Lau, T.S., 1988. D-optimal designs on the unit q-ball. Journal of Statistical Planning and Inference 19, 299–315. Lau, T.S., Studden, W.J., 1985. Optimal designs for trigonometric and polynomial regression using canonical moments. Annals of Statistics 13, 385–394. Lim, Y.B., Studden, W.J., 1988. Efficient Ds-optimal designs for multivariate polynomial regression on the q-cube. Annals of Statistics 16, 1225–1240. Nakamura, Y., 1988. A tau-function of the finite nonperiodic Toda lattice. Physical Letters A 195 (5–6), 346–350. Pukelsheim, F., 1993. Optimal Design of Experiments. Wiley, New York. Skibinsky, M., 1967. The range of the (nþ 1)th moment for distributions on [0, 1]. Journal of Applied Probability 4 (3), 543–552. Skibinsky, M., 1968. Extreme nth moments for distributions on ½0; 1 and the inverse of a moment space map. Journal of Applied Probability 5 (3), 693–701. Skibinsky, M., 1969. Some striking properties of binomial and beta moments. Annals of Mathematical Statistics 40 (5), 1753–1764. Spiridonov, V., Zhedanov, A., 1995. Discrete Darboux transformations, the discrete-time Toda lattice, and the Askey–Wilson polynomials. Methods and Applications of Analysis 2, 369–398. Spiridonov, V., Zhedanov, A., 1997. Discrete-time Volterra chain and classical orthogonal polynomials. Journal of Physics A: Mathematical and General 30, 8727–8737. Studden, W.J., 1980. Ds-optimal designs for polynomial regression using continued fractions. Annals of Statistics 8 (5), 1132–1141. Studden, W.J., 1982a. Optimal designs for weighted polynomial regression using canonical moments. Statistical Decision Theory and Related Topics III 2, 335–350. Studden, W.J., 1982b. Some robust-type D-optimal design in polynomial regression. Journal of American Statistical Association 77, 916–921. Studden, W.J., 1989. Note on some Fp optimal design for polynomial regression. Annals of Statistics 17, 618–623. Symes, W.W., 1982. The QR algorithm and scattering for the finite nonperiodic Toda lattice. Physica D 4 (2), 275–280. Tsujimoto, S., Nakamura, Y., Iwasaki, M., 2001. The discrete Lotka–Volterra system computes singular values. Inverse Problems 17 (1), 53–58. Zen, M.M., Tsai, M.H., 2004. Criterion-robust optimal designs for model discrimination and parameter estimation in Fourier regression models. Journal of Statistical Planning and Inference 124 (2), 475–487.