To apply equivalence theorem to a wide variety of practical problems, it is essential to have a computer algorithm to search for these optimal minimax designs numerically, because it is hard to find these optimal minimax designs when the order of the polynomial models gets higher. Therefore, we introduce an algorithm for generating the optimal minimax designs numerically, and this algorithm is based on the iterative scheme of Wong (1998).
First we transform our design problems into the optimization problems with the objective function, Φi, i = 1, 2, 3. Since we require the mass w of the design to satisfy 0 < w ≤ 1 and the supports x to satisfy x ∈ X , we need to add a penalty function to the objective function Φ to limit constraints violation. For example when X = [−1, 1] and there is only one mass w, the penalty can be chosen as
L = r[(max[0, x − 1])2+ (max[0, −1 − x])2+ (max[0, w − 1])2+ (max[0, −w])2], where r is some appropriately chosen positive number. In Wong (1998), Powell’s method and Golden section method is applied to solve our constrained optimization problem, and usually the penalty functions are chosen to be the squared functions. Here Powell’s method (Powell, 1964) is based on the concept of conjugate directions, where direction S1 and S2 are said to be conjugate if S1THS2 = 0 and H is some approximation to the Hessian matrix. A particularly attractive and significant feature of Powell’s method is that, if the objective function is quadratic in each of the n variables, then the function will be minimized in n or fewer conjugate search directions. Please see Vanderplaats (1984) for further details. Therefore, when we give a set of initial design parameters P0
in Powell’s method, the optimizer will output a set of design parameters P∗, which may be a local optimum for the objective function. The optimization algorithm is an iterative application of the equation:
Pt= Pt−1+ α∗tSt,
where Pt is the set of design variables, St is the search direction, and α∗t is a scalar determining the optimal amount of change in St at the t-th iteration. Here we use the Golden Section method to determine the best α∗t for given St. Finally, since designs found may be locally optimal, Theorem 2.1 is used to verify if the design is optimal within the class of all designs.
The generating algorithm in Wong (1998) is described as follow. Since the algorithm is an iterative method, we use Ptto denote the design at the t-th iteration. The algorithm for finding an MV-optimal design for fT(x) = (1, x, . . . , xd) on X = [−1, 1] is in the following :
Step 0. Set objective function: Φ1 and variables: d, Pt, λ(x), Z, ε.
Step 1. Initialize St to be coordinate unit vectors, t = 1, . . . , 2d + 1.
Step 2. t ← 0.
Step 3. t ← t + 1.
Step 4. Find α∗t to minimize Φ1(Pt−1+ αtSt) + rL using the Golden Section method.
Step 5. Pt← Pt−1+ α∗tSt
Step 6. Check if t = 2k + 1. If yes, proceed to Step 7. Otherwise, proceed to Step 3.
Step 7. Create a conjugate direction St+1 ← Pt− P0
Step 8. Find αt+1∗ to minimize Φ1(Pt+ αt+1St+1) + rL using the Golden Section method.
Step 9. π∗ ← Pt+ α∗t+1St+1.
Step 10. Check if π∗ a local minimum by Powell’s method(εp = 10−5). If yes, proceed to Step 11. Otherwise, let St ← St+1, t = 1, . . . , 2k + 1 and P0 = π∗, then proceed to Step 3.
Step 11. Check if π∗ a global minimum.
Step 12. Find u1, u2, . . . in A1(π∗) to minimize −d(z, π∗) using the Golden Section method.
Step 13. Find µ∗ to minimize c1(x∗, µ, π∗) using Powell’s method (x∗: any support point of π∗).
Step 14. Check if c1(x∗, µ∗, π∗) ≤ ε for all x ∈ [−1, 1]. If yes, stop and conclude that π∗ is MV-optimal for the given tolerance level ε = 10−3. Otherwise, Randomly select another design P0, then proceed to Step 1.
For generating the other two types minimax designs, we need to replace A1 to be A2 or A3; to change the objective function Φ1 to be Φ2 or Φ3, and the corresponding equivalence theorem c1 to be c2 or c3.
Here we simply implement the algorithm to find an MV-optimal design for a simple linear model with a given λ(x) and Z. This is a 3-dimensional problem which requires input variables: x1, x2 and w1. Suppose we choose P0 = (x1, x2, w1) = (1, −1, 0.5).
Then the first cycle of minimization will be in the coordinate directions S1 = (1, 0, 0), S2 = (0, 1, 0) and S3 = (0, 0, 1) to determine α∗i, i = 1, . . . , 2d + 1, sequentially in accordance to
Pt= Pt−1+ α∗tSt, where each αt∗ is selected to minimize
Φ1(Pt−1+ αtSt) + rL
over the real line. Then a conjugate direction S4 = X3
i=1
α∗iSi is formed, and α∗4 is de-termined to yield a new design with parameters P∗ = P4∗ = P3 + α∗4S4. Because of the constraints imposed on the design problem (−1 ≤ x1 < x2 ≤ 1 and 0 < w1 ≤ 1), our objective is to minimize
Φ1+ rL,
where r is a multiplier and L is the penalty function defined by L =
X2 i=1
[(max[0, xi− 1])2 + (max[0, −1 − xi])2] + (max[0, w1− 1])2+ (max[0, −w1])2. Note that L does not penalize Φ when the solution is in the feasible set for the problem considered here, i.e. −1 ≤ x1 < x2 ≤ 1 and 0 < w1 ≤ 1.
3 Optimal minimax designs for simple linear model
In this section, we consider the simple linear model, i.e. fT(x) = (1, x), on X = [−1, 1] and three different types efficiency function λ(x) are chosen. These three efficiency
functions are linear function, x + 5; concave function, exp(−5x2), and convex function, 0.5x2 + 1. Here for the intervals Z, we choose three types. The first one is symmetric interval, [−0.8, 0.8], [−1, 1] and [−2, 2]; the second type is the extrapolation interval, [1, 1.5], and the third type of Z is [−1, 1.5].
3.1 Numerically optimal minimax designs
Now the generating algorithm is applied to find the optimal minimax designs nu-merically. In the generating algorithm, first the starting design is chosen by randomly picking two distinct support points with equal weights 1/2, and we fix r = 5000 in the generating algorithm for all cases.
Given different λ(x) and Z, the corresponding MV-optimal, MP-optimal and MD-optimal designs on X = [−1, 1] are shown in following:
1. λ(x) = x + 5 : The first case is that λ(x) = x + 5. To make sure that the numerically found designs are indeed optimal, we need to check the corresponding ci(x, µ∗, π∗), i = 1, 2, 3, to see if ci(x, µ∗, π∗) ≤ 0, ∀x ∈ X . For example, for MV-optimal design with λ(x) = x + 5 and Z = [−0.8, 0.8], the numerical MV-opitmal design is πM V =
( −1 1 0.6 0.4
)
with µ∗ =
( −0.8 0.8 0.397 0.603
)
, and
maxz∈Z d(z, πM V) = Φ1(±0.8, πM V) = 0.3417.
Now
c1(x, µ∗, πM V) = −0.1247 − 0.0278x + 0.1247x2 + 0.0278x3.
From Figure 1, it is clear that πM V is MV-optimal. For the other cases, the figures of ci(x, µ∗, π∗)’s are displayed in Appendix A.1.1. The numerically optimal designs for these three criteria are shown in Table 1.
-1 -0.5 0.5 1
-0.12 -0.1 -0.08 -0.06 -0.04 -0.02
Figure 1: Plot for c1(x, µ∗, πM V)
Table 1 : λ(x) = x + 5
Z πM V µ∗ πM P µ∗ πM D µ∗
[−0.8, 0.8] 1 0.4 0.8 0.397 1 0.352 0.8 0.301 1 0.480 0.8 0.483
-1 0.6 -0.8 0.603 -1 0.648 -0.8 0.699 -1 0.520 -0.8 0.517
[−1, 1] 1 0.4 1 0.4 1 0.356 1 0.315 1 0.5 1 0.5
-1 0.6 -1 0.6 -1 0.644 -1 0.685 -1 0.5 -1 0.5
[−2, 2] 1 0.4 2 0.375 1 0.348 2 0.25 1 0.667 2 0.885
-1 0.6 -2 0.625 -1 0.652 -2 0.75 -1 0.333 -2 0.115
[1, 1.5] 1 0.803 1.5 1 1 0.803 1.5 1 1 0.803 1.5 1
-1 0.197 -1 0.197 -1 0.197
[−1, 1.5] 1 0.526 1.5 0.561 1 0.48 1.5 0.466 1 0.653 1.5 0.727
-1 0.474 -1 0.439 -1 0.52 -1 0.534 -1 0.347 -1 0.273
Table 2 : λ(x) =exp(−5x2)
Z πM V µ∗ πM P µ∗ πM D µ∗
[−0.8, 0.8] 0.4 0.5 0.8 0.5 0.4 0.5 0.8 0.5 0.316 0.5 0.316 0.5
-0.4 0.5 -0.8 0.5 -0.4 0.5 -0.8 0.5 -0.316 0.5 -0.316 0.5
[−1, 1] 0.415 0.5 1 0.5 0.416 0.5 1 0.5 0.316 0.5 0.316 0.5
-0.415 0.5 -1 0.5 -0.416 0.5 -1 0.5 -0.316 0.5 -0.316 0.5
[−2, 2] 0.436 0.5 2 0.5 0.437 0.5 2 0.502 0.316 0.5 0.316 0.5
-0.436 0.5 -2 0.5 -0.437 0.5 -2 0.498 -0.316 0.5 -0.316 0.5 [1, 1.5] 0.447 0.649 1.5 1 0.447 0.649 1.5 1 0.447 0.724
-0.447 0.351 -0.447 0.351 -0.447 0.276 1 1
[−1, 1.5] 0.447 0.649 1.5 1 0.477 0.649 1.5 1 0.316 0.5 0.316 0.5
-0.447 0.351 -0.477 0.351 -0.316 0.5 -0.316 0.5
2. λ(x) = exp(−5x2) : Wong (1998) has also studied this case. From his results, we assume that the support points of optimal minimax designs are symmetric for saving computing time and cost. That is, in the algorithm, (x1, x2, w1) is replaced with (x1, −x1, w1). To check if numerical designs are optimal, we still plot ci(x, µ∗, π∗), and all figures of ci(x, µ∗, π∗)’s are shown in Appendix A.1.2. The corresponding numerically optimal minimax designs are in Table 2.
3. λ(x) = 0.5x2+ 1 : Here λ(x) = 0.5x2+ 1 is a convex function, and the corresponding numerical results are shown in Table 3. We also plot ci(x, µ∗, π∗) to make sure our numerical results are optimal, and all figures are shown in Appendix A.1.3.
From our results, three numerically optimal minimax designs are the same when Z = [−0.8, 0.8], [−1, 1], [−2, 2], [1, 1, 5] and [−1.5, −1].
According to our numerically optimal designs, we compare the MV-, MP- and
MD-Table 3 : λ(x) = 0.5x2 + 1
Z πM V µ∗ πM P µ∗ πM D µ∗
[−0.8, 0.8] 1 0.5 0.8 0.5 1 0.5 0.8 0.5 1 0.5 0.8 0.5
-1 0.5 -0.8 0.5 -1 0.5 -0.8 0.5 -1 0.5 -0.8 0.5
[−1, 1] 1 0.5 1 0.5 1 0.5 1 0.5 1 0.5 1 0.5
-1 0.5 -1 0.5 -1 0.5 -1 0.5 -1 0.5 -1 0.5
[−2, 2] 1 0.5 2 0.5 1 0.5 2 0.5 1 0.5 2 0.5
-1 0.5 -2 0.5 -1 0.5 -2 0.5 -1 0.5 -2 0.5
[1, 1.5] 1 0.833 1.5 1 1 0.833 1.5 1 1 0.833 1.5 1
-1 0.167 -1 0.167 -1 0.167
[−1.5, −1] 1 0.167 1 0.167 1 0.167
-1 0.833 -1.5 1 -1 0.833 -1.5 1 -1 0.833 -1.5 1
[−1, 1.5] 1 0.625 1.5 0.667 1 0.597 1.5 0.607 1 0.708 1.5 0.777
-1 0.375 -1 0.333 -1 0.403 -1 0.393 -1 0.292 -1 0.223
Table 4 : Φi-relative efficiency for λ(x) = x + 5
Φ1 Φ2 Φ3
Z πM P πM D πM V πM D πM V πM P
[−0.8, 0.8] 0.882 0.870 0.966 0.887 0.837 0.738
[−1, 1] 0.890 0.833 0.957 0.851 0.800 0.712
[−2, 2] 0.888 0.595 0.957 0.631 0.721 0.640
[1, 1.5] 1 1 1 1 1 1
[−1, 1.5] 0.922 0.733 0.940 0.753 0.841 0.775
optimal designs by their relative efficiencies, and these relative efficiencies are shown in Tables 4, 5 and 6. Here we highlight these results in the following :
• λ(x) = x + 5 : From Table 4, we get that three optimal minimax designs are the same with Z = [1, 1.5]. For the other intervals, the Φ1-relative efficiencies of πM P and Φ2 -relative efficiencies of πM V are all larger than 0.88. Therefore, we think that πM V’s and πM P’s have the similar performances for simple linear model. However, for πM D, these corresponding relative efficiencies are low. Especially, when Z = [−2, 2], the Φ1- and Φ2-relative efficiencies of πM D are 0.595 and 0.631, respectively.
• λ(x) = exp(−5x2) : As shown in Table 5, πM V’s are equal to πM P’s for all Z’s. For πM D, the Φ1- and Φ2-relative efficiencies of πM D are at least higher than 0.78. That is, three optimal minimax designs are about equal for MV- and MP-criteria. For the Φ3-relative efficiencies, the corresponding values of πM V and πM P are all higher than 0.75, except the case of Z = [−1, 1.5]. When Z = [−1, 1.5], these two values
Table 5 : Φi-relative efficiency for λ(x) = exp(−5x2)
Φ1 Φ2 Φ3
Z πM P πM D πM V πM D πM V πM P
[−0.8, 0.8] 1 0.912 1 0.887 0.878 0.878
[−1, 1] 1 0.887 1 0.988 0.836 0.834
[−2, 2] 1 0.843 1 0.999 0.772 0.769
[1, 1.5] 1 0.973 1 0.999 0.976 0.976
[−1, 1.5] 1 0.789 1 0.999 0.556 0.556
Table 6 : Φi-relative efficiency for λ(x) = 0.5x2+ 1
Φ1 Φ2 Φ3
Z πM P πM D πM V πM D πM V πM P
[−0.8, 0.8] 1 1 1 1 1 1
[−1, 1] 1 1 1 1 1 1
[−2, 2] 1 1 1 1 1 1
[1, 1.5] 1 1 1 1 1 1
[−1.5, −1] 1 1 1 1 1 1
[−1, 1.5] 0.962 0.779 0.949 0.787 0.908 0.873
are 0.556 for both πM V and πM P.
• λ(x) = 0.5x2+ 1 : According to Table 6, three optimal minimax designs are the same except the case of Z = [−1, 1.5]. For Z = [−1, 1.5], all these corresponding relative efficiencies are all larger than 0.75. Therefore, we think that three optimal minimax designs have the similar performance.
From Tables 4, 5 and 6, three optimal minimax designs are equivalent when λ(x) = 0.5x2+ 1 and Z = [−0.8, 0.8], [−1, 1], [−2, 2], [1, 1.5] and [−1.5, −1]. Here we conjecture that these three types of optimal minimax designs are the same when λ(x) is a convex function, and Z is a symmetric interval or an extrapolation interval. This is denoted as conjecture 1.
Conjecture 1. Let design space X = [−1, 1] and fT(x) = (1, x). Suppose efficiency function is symmetric convex function over X , and Z is a symmetric interval or an extrapolation interval. Then these three types of optimal minimax designs are equivalent.
Thus, two more cases for symmetric convex efficiency functions are studied here. In these two case, efficiency functions are set to be x2 + 1 and 10x2 + 1 with Z = [−1, 1] and
Table 7 : λ(x) = cx2+ 1
c = 1 πM V µ∗ πM P µ∗ πM D µ∗
Z = [−1, 1] 1 0.5 1 0.5 1 0.5 1 0.007 1 0.5 1 0.5 0 0.986
-1 0.5 -1 0.5 -1 0.5 -1 0.007 -1 0.5 -1 0.5
Z = [1, 1.5] 1 0.833 1.5 1 1 0.833 1.5 1 1 0.833 1.5 1
-1 0.167 -1 0.167 -1 0.167
c = 10 πM V µ∗ πM P µ∗ πM D µ∗
Z = [−1, 1] 1 0.5 1 0.5 1 0.5 0 1 1 0.5 1 0.5
-1 0.5 -1 0.5 -1 0.5 -1 0.5 -1 0.5
Z = [1, 1.5] 1 0.833 1.5 1 1 0.833 1.5 1 1 0.833 1.5 1
-1 0.167 -1 0.167 -1 0.167
[1, 1.5]. The numerical results are shown in Table 7. From Table 7, numerically optimal designs of three minimax criteria are still equivalent for Z = [−1, 1] and [1, 1.5].