AN APPROACH TO OPTIMIZATION IN RESPONSE SURFACE METHODOLOGY.

By JOSÉ A. DÍAZ-GARCÍA ^{∗}

Department of Statistics and Computation Universidad Autónoma Agraria Antonio Narro 25350 Buenavista, Saltillo, Coahuila, MÉXICO.

and

ROGELIO RAMOS-QUIROGA

Department of Probability and Statistics Centro de Investigación en Matemáticas A. C. 36240 Guanajuato, MÉXICO

Key Word and Phrases: Stochastic optimization, optimum conditions, response surfaces, doubly noncentral tt-distribution.

ABSTRACT

We consider the optimization of a response surface given by a polynomial model. In this article the problem is stated as a stochastic optimization problem. For the first-order model, the E-Model technique is used to find the maximum and stationary points; the exact distributions are derived under the linearity assumption.

1. INTRODUCTION

Consider the model Y = z′(x)β + ϵ, where ϵ ∼ N(0, σ^{2}) and z(x) contains polynomial terms in the covariables. The usual approach in the last steps in the methodology of response surfaces is to obtain the least squares fit, Ŷ (x), and then proceed to the optimization phase by computing the critical point, x̂_{0}, such that it maximizes (minimizes) the response, possibly subject to restrictions on the factor space. We can state this as

max_{x∈Ξ} Ŷ (x)  (1)

∗This article was written by the author during his stay as Visiting Professor in the Dept. of Mathematics at the ITESM, Monterrey, México and at CIMAT, Guanajuato, México.

where Ξ is the experimental region, usually defined as an hypercube or an hypersphere. To solve this problem several approaches have been proposed in the literature, see Khuri and Cornell (1987, Chapter 5), all these relay on ideas of mathematical programming, i.e. deterministic optimization. Specifically, methods of differential calculus are applied to optimize an estimated second-order response surface, see Khuri and Cornell (1987, pp. 174-175) and Lagrange methods are applied for the corresponding problem restricted to a spherical region, see Khuri and Cornell (1987, pp. 158-160). It is important to note that the problem stated in (1) depends on parameters that have been estimated and therefore the objective function is not deterministic. The last four decades have seen the development of the area of mathematical programming called Stochastic Programming; this area, which can be considered to be a part of statistics, studies mathematical programming problems where some of the parameters or the decision variables themselves are random variables. Garstka and Wets (1974), classify stochastic optimization problems into two clases: the penalty models called stochastic programs with recourse and the reliability models called stochastic programs with change constraints. It can be shown that, under certain conditions, these two problems are equivalent in the sense that each one is a special case of the other. In general, a stochastic optimization problem can be stated as

max_{x∈D} f(x, α)

where α is a (possibly vectorial) random variable and D is determined by a set of restrictions which can be of stochastic nature. Prékopa (1978), restated several optimization problems in statistics as problems in stochastic optimization. For the optimization problem in the area of response surface methodology, note that it really is a problem of stochastic optimization. Following the notation in Khuri and Cornell (1987) and assuming a second-order polynomial model we restate (1) as

max_{x∈Ξ} Ŷ (x, β̂) = z′(x)β̂, β̂ ∼ N_{k}(β, σ^{2}(X′X)^{−1})  (2)

this is readily recognized as a stochastic optimization problem. In this paper the E-Model technique of stochastic optimization is used to find the maximum and stationary point for a first-order response surface model (see Theorem 1). In sections 3 and 4 we obtain the exact distributions for the maximum and stationary point respectively; the expressions for these are obtained in section 2. Finally, we present in the appendix the density for the doubly noncentral tt distribution.

2. STATEMENT OF THE PROBLEM.

Essentially, the stochastic optimization methods transform a stochastic problem into its “equivalent” deterministic one, Charnes and Cooper (1963). Under the E-Model technique we can restate problem (2) as the following deterministic counterpart (Charnes and Cooper (1963))

max_{x∈Ξ} E(Ŷ (x, β̂)) = z′(x)β,  (3)

in what follows, we will assume that Ξ = {x : ||x||^{2} ≤ c^{2}}. To illustrate the E-Model technique, assume that E(Ŷ (x, β̂)) = β′x i.e. a linear model (an hyperplane). In this case the maximum is reached on the border of Ξ, i.e. ||x||^{2} = c^{2}; so, applying the Lagrange multipliers technique (note that the problem is already deterministic), we have

L(x, λ) = β′x + λ(||x||^{2} − c^{2})

taking derivatives and equating them to zero we have

\begin{aligned}\frac{∂L(x, λ)}{∂x} & = & β + 2λx & = & 0 \\ \frac{∂L(x, λ)}{∂λ} & = & ||x||^{2} − c^{2} & = & 0 \end{aligned}

if we denote the solution to this system by (x_{0}, λ_{0}) then x_{0} = −β/2λ and

||x_{0}||^{2} = (−\frac{β}{2λ})′ (−\frac{β}{2λ}) = \frac{||β||^{2}}{4λ^{2}} = c^{2}

hence

λ = ±\sqrt{\frac{||β||^{2}}{4c^{2}}}

since λ ≤ 0 and both 4c^{2} and ||β||^{2} are > 0, we have that the critial point is given by

(x_{0}, λ_{0}) = (\frac{cβ}{||β||}, −\frac{||β||}{2c}) .

if we denote the maximum of E(Ŷ (x, β̂)) by E_{max}(Ŷ (x_{0}, β̂)), we have that

E_{max}(Ŷ (x_{0}, β̂)) = β′ (\frac{cβ}{||β||}) = c||β||.

Note that the maximum E_{max}(Ŷ (x_{0}, β̂)) as well as the critical point x_{0} are parametric functions, i.e., they depend on β; therefore the maximum likelihood estimators for them will be based on the maximum likelihood estimator β̂. The estimators for E_{max}(Ŷ (x_{0}, β̂)) and x_{0} are respectively

Ŷ_{max}(x̂_{0}) = c||β̂|| y x̂_{0} = \frac{cβ̂}{||β̂||}  (4)

where, for simplicity, we denoted Ê_{max}(Ŷ (x_{0}, β̂)) by Ŷ_{max}(x̂_{0}). Therefore, we have the following result.

Theorem 1. For a response surface model, with Ξ = {x : ||x||^{2} ≤ c^{2}}, we have that

max_{x∈Ξ} E(Ŷ (x, β̂))

is given by Y_{max}(x_{0}) = c||β|| with x_{0} = \frac{cβ}{||β||} and their corrresponding estimators by (4).

Remark 1. Observe that, without loss of generality, we can choose c = 1, this can be done merely by standarizing the variables x′s dividing by c. We will assume this standarization in the rest of the paper.

3. DISTRIBUTION OF THE MAXIMUM, Ŷ_{max}(r̂_{0}).

It is not straightforward to establish the distributions of the maximum and critical point when we see them solely as functions of β̂. We consider a reparameterization of the original model that facilitates this task. Let us consider the linear model

Y = Xβ + ϵ, ϵ ∼ N_{n}(0, σ^{2}I_{n})

let Γ be the nonnegative definite square root of (X′X)^{−1}, i.e. Γ^{2} = (X′X)^{−1} and consider the following reparameterization

Y = XΓΓ^{−1}β + ϵ = Rδ + ϵ

where R = XΓ and δ = Γ^{−1}β, for which we have δ̂ = R^{−}Y = (XΓ)^{−}Y = Γ^{−}X^{−} = Γ^{−1}β̂ so that δ̂ ∼ N_{k}(δ, σ^{2}I_{k}).

Now let us consider the linear model for the response surface

Y (x) = x′β + e, e ∼ N (0, σ^{2})

with the corresponding reparameterization

Y (r) = r′δ + e, r = Γx, δ = Γ^{−1}β.

for which the estimated response is given by

Ŷ (r) = r′δ̂,  Ŷ (r) ∼ N (r′δ, σ^{2}||r||^{2})

thus, we can see that

max_{r∈Ξ_{r}} E(Y (r))

where Ξ_{r} = {r : ||r||^{2} = 1}, is given by

E_{max}(Y (r_{0})) = ||δ||, where r_{0} = \frac{δ}{||δ||}

and their respective estimators are

Ŷ_{max}(r̂_{0}) and r̂_{0} = \frac{δ̂}{||δ̂||}.

To get the distribution of the maximum recall that δ̂ ∼ N_{k}(δ, σ^{2}I_{k}), from where we get

u = \frac{||δ̂||^{2}}{σ^{2}} ↦ χ^{2}(u; k, λ), λ = \frac{1}{2σ^{2}} ||δ||^{2}

where ↦ means “with density”, and χ^{2}(u; k, λ) is the noncentral χ^{2} density with k degrees of freedom and noncentrality parameter λ. So that

v = ||δ̂||^{2} ↦ \frac{1}{σ^{2}} χ^{2} (\frac{v}{σ^{2}} ; k, λ)

finally, we have

\sqrt{v} = ||δ̂|| ↦ \frac{2v}{σ^{2}} χ^{2} (\frac{v^{2}}{σ^{2}} ; k, λ)  (5)

In summary, we have the following result.

Theorem 2.Under the same assumptions of Theorem 1, with the reparameterized model, we have that the distribution of the estimated maximum is given by (5).

4. DISTRIBUTION OF THE CRITICAL POINT, r̂_{0} .

It is easy to see that δ̂ and ||δ̂|| are not independent and therefore this makes it more difficult to obtain the distribution of r̂_{0}, however, we can get the marginal distribution of r̂_{i0}, i = 1, …, k, by extending to the noncentral case a procedure presented by Cramér (1974, pp. 240-241). To do this, define θ_{i} as

θ_{i} = \sqrt{k}r̂_{i0} = \frac{\sqrt{k}δ̂_{i}}{||δ̂||}

note that θ_{i}^{2} ≤ k or, equivalently |θ_{i}| ≤ \sqrt{k}. Hence the density of θ_{i} is zero outside the interval (−\sqrt{k}, \sqrt{k}).

Now, let

v_{i} = \sqrt{\frac{k − 1}{k}} \frac{θ_{i}}{\sqrt{(1 − \frac{θ_{i}^{2}}{k} )}}

then

\begin{aligned}v_{i} & = & \frac{\sqrt{k − 1}θ_{i}}{\sqrt{k}\sqrt{\frac{k − θ_{i}^{2}}{k}}} & = & \frac{\sqrt{k − 1}θ_{i}}{\sqrt{k − θ_{i}^{2}}} & = & \frac{\frac{k\sqrt{k − 1}δ̂_{i}}{||δ̂||}}{\sqrt{k − \frac{kδ̂_{i}^{2}}{||δ̂||^{2}}}} \\ & = & \frac{δ̂_{i}}{\sqrt{\frac{1}{k − 1}(||δ̂||^{2} − δ̂_{i}^{2})}} & = & \frac{δ̂_{i}}{\sqrt{\frac{1}{k − 1}∑_{ν≠i}^{k}δ̂_{ν}^{2}}} & & . \end{aligned}

Observe that δ̂_{i} ∼ N (δ_{i}, σ^{2}) independently of

s = ∑_{ν≠i}^{k} δ̂_{ν}^{2} ↦ \frac{1}{σ^{2}} χ^{2} (\frac{s}{σ^{2}} , (k − 1), λ^{∗}) , λ^{∗} = \frac{1}{2σ^{2}} ||δ^{∗}||^{2}

from where δ^{∗′} = (δ_{1} … δ_{i−1}δ_{i+1} … δ_{k}). Hence

\frac{δ̂_{i}}{σ} ∼ N (\frac{δ_{i}}{σ} , 1)

\frac{s}{σ^{2}} ∼ χ_{(k−1),λ^{∗}}^{2}

therefore

\frac{\frac{δ̂_{i}}{σ}}{\sqrt{\frac{∑_{ν≠i}^{k}δ̂_{ν}^{2}}{(k − 1)σ}}} = \frac{δ̂_{i}}{\sqrt{\frac{∑_{ν≠i}^{k}δ̂_{ν}^{2}}{(k − 1)}}} = v_{i} ↦ tt′′(v_{i}; (k − 1), λ_{1}, λ^{∗})

where tt′′(v_{i}; (k − 1), λ_{1}, λ^{∗}) is the doubly noncentral t density function with (k − 1) degrees of freedom and noncentrality parameters λ_{1} = δ_{i}/σ and λ^{∗} (see Appendix).

Note that if θ_{i} varies in the interval (−\sqrt{k}, \sqrt{k}), v_{i} will range on the whole real line (−∞, ∞), thus

v_{i} < \sqrt{\frac{k − 1}{k}} \frac{u}{\sqrt{1 − \frac{u^{2}}{k}}}

is equivalent to have θ_{i} < u, so

P (θ_{i} < u) = P (v_{i} < \sqrt{\frac{k − 1}{k}} \frac{u}{\sqrt{1 − \frac{u^{2}}{k}}} ) = ∫−∞^{\sqrt{\frac{k−1}{k}} \frac{u}{\sqrt{1− \frac{u^{2}}{k}}}} tt′′(v_{i}; (k − 1), λ_{1}, λ^{∗})dv_{i}

taking the derivative with respect to u we obtain the density of θ_{i}, thus

θ_{i} = \sqrt{k}r̂_{0i} = \frac{δ̂_{i}}{\sqrt{\frac{∑_{ν≠i}^{k}δ̂_{ν}^{2}}{(k − 1)}}} ↦ \sqrt{\frac{k − 1}{k}} (1 − \frac{u^{2}}{k} )^{−3/2} tt′′ (\sqrt{\frac{k − 1}{k}} \frac{u}{\sqrt{1 − \frac{u^{2}}{k}}} ; (k − 1), λ_{1}, λ^{∗})

(6)

where |u| ≤ \sqrt{k}. Hence, we have the following result.

Theorem 3.Under the same assumptions of Theorem 1, with the reparameterized model, we have that the distribution of the critical point is given by (6).

APPENDIX. DOUBLY NONCENTRAL tt DISTRIBUTION

In this section we present an expression for the doubly noncentral tt density function.

Lemma. Let Y ∼ N (ν, σ) and suppose X/σ^{2} ∼ χ_{n,λ}^{2} independently of Y. Define the random variable

U = \frac{\frac{Y}{σ}}{\sqrt{\frac{X}{nσ^{2}}}} = \frac{Y}{\sqrt{\frac{X}{n}}}

the distribution of U is the doubly noncentral tt distribution with n degrees of freedom and noncentrality parameters λ_{1} = μ/σ and λ. The density of this distribution is given by

f_{U} (u) = ∑_{w=0}^{∞} P (W = w)tt′(u; n + 2w, λ_{1}).

where tt′(u; n + 2w, λ_{1}) is the noncentral tt density function with n + 2w degrees of freedom and noncentrality parameter λ_{1}, see Rao (1973, section 3a.5(iii)), and P (W = w) is the Poisson probability function with parameter λ, see Rao (1973, p. 88). We denote this as

U ∼ tt′′(v_{i}; (k − 1), λ_{1}, λ)

Proof. Let N (y; μ, σ^{2}) denote the normal density. We know that the density of the noncentral χ^{2} distribution with n degrees of freedom and noncentrality parameter λ is given by

f_{S}(s) = ∑_{w=0}^{∞} \frac{P (W = w)}{σ^{2}}χ^{2}(s/σ^{2}; n + 2w)

where P (W = w) is the Poisson probability function with parameter λ and χ^{2}(s; n + 2w) is the central χ^{2} density function with n + 2w degrees of freedom. Then, the joint density of Y and X is given by

f_{Y,X} (y, x) = ∑_{w=0}^{∞} \frac{P (W = w)}{σ^{2}}χ^{2}(s; n + 2w)N (y; μ, σ^{2}).

If we define U = Y/\sqrt{X/n}, then we note that the change of variable will only involve the two continuos densities, therefore the problem is reduced to find the density of the noncentral tt distribution, hence (see Rao (1979, section 3a.5(iii))) we get the desired result. □

References

[1] Charnes, A., and Cooper, W. W. (1963). “Deterministic equivalents for optimizing and satisficing under chance constraints”, Operation Research, 11, 18-39.

[2] Cramér, H. (1974). Mathematical Methods of Statistics, Princeton University Press, Princeton, N. J.

[3] Garstka, S. J., and Wets, R. J.-B. (1974). “On decision rules in stochastic programming”, Mathematical Programming, 7, 117-143.

[4] Khuri, A. I., and Conlon, J. A. (1981). “Simultaneous Optimization of Multiple Responses Represented by Polynomial Regression Functions”, Technometrics, 23, 363-375.

[5] Khuri, A. I., and Cornell, J. A. (1987). Response Surfaces: Designs and Analyses, Marcel Dekker, Inc. New York.

[6] Prékopa, A. (1978).“The use of stochastic programming for solution of some problems in statistics and probability ”, Technical summary report #1834, University of Wisconsin-Madison, Mathematics Research Center.

[7] Rao, R. (1973). Linear Statistical Inference and Its Applications, John Wiley & Sons, New York.