Arts and Sciences, Irving K. Barber School of (Okanagan)
Computer Science, Mathematics, Physics and Statistics, Department of (Okanagan)
DSpace
UBCO
Foster, Jodie Lynn
Master of Science - MSc
University of British Columbia
Mathematical models implemented as computer code are gaining widespread use across the sciences and engineering. In some cases these model are even replacing physical experiments. The computational complexity of the models typically limits the number of runs that can be performed. In such cases, the Gaussian process is used to emulate the true computer model so that experiments are performed on the emulator. The applicability of the emulator is closely related to the quality of the fitted Gaussian process model. A key step in fitting the Gaussian process is estimating the unknown correlation parameters using a numerical optimizer. It is well known that maximum likelihood estimation is invariant to the parameterization of the model. When the mapping from one parametrization to another is injective this leads to an equivalence of the likelihood and in a Bayesian context an equivalence of the posterior for carefully chosen prior distributions. In the context of a Gaussian process model, we show that the parameterization of the model is in fact critical to achieving meaningful results and ensuring convergence of the optimizer. This thesis is aimed at providing practical advice on the best parameterization for the Gaussian process model. The approach presented here implements a simulation study on a wide range of test problems and a large number of parameterizations. We illustrate that the parameterization can have a huge effect on the fitted model and show that many of the commonly used parameterizations are in fact sub-optimal for fitting the Gaussian process.
A Numerical Study of the Effects ofthe Parameterization of the GaussianProcessbyJodie Lynn FosterB.Sc., The University of British Columbia, 2013A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe College of Graduate Studies(Computer Science)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)August 2017© Jodie Lynn Foster, 2017The undersigned certify that they have read, and recommend to theCollege of Graduate Studies for acceptance, a thesis entitled: A NumericalStudy of the Effects of the Parameterization of the GaussianProcess submitted by Jodie Lynn Foster in partial fulfilment of therequirements of the degree of Master of ScienceJason Loeppky, Irving K. Barber School of Arts and SciencesSupervisor, ProfessorJohn Braun, Irving K. Barber School of Arts and SciencesSupervisory Committee Member, ProfessorParamjit Gill, Irving K. Barber School of Arts and SciencesSupervisory Committee Member, ProfessorJulian Cheng, School of EngineeringUniversity Examiner, ProfessorJuly 31, 2017(Date Submitted to Grad Studies)iiAbstractMathematical models implemented as computer code are gaining widespreaduse across the sciences and engineering. In some cases these model areeven replacing physical experiments. The computational complexity of themodels typically limits the number of runs that can be performed. In suchcases, the Gaussian process is used to emulate the true computer modelso that experiments are performed on the emulator. The applicability ofthe emulator is closely related to the quality of the fitted Gaussian processmodel. A key step in fitting the Gaussian process is estimating the unknowncorrelation parameters using a numerical optimizer. It is well known thatmaximum likelihood estimation is invariant to the parameterization of themodel. When the mapping from one parametrization to another is injectivethis leads to an equivalence of the likelihood and in a Bayesian context anequivalence of the posterior for carefully chosen prior distributions. In thecontext of a Gaussian process model, we show that the parameterizationof the model is in fact critical to achieving meaningful results and ensuringconvergence of the optimizer. This thesis is aimed at providing practicaladvice on the best parameterization for the Gaussian process model. Theapproach presented here implements a simulation study on a wide range oftest problems and a large number of parameterizations. We illustrate thatthe parameterization can have a huge effect on the fitted model and showthat many of the commonly used parameterizations are in fact sub-optimalfor fitting the Gaussian process.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . ixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 1Chapter 2: The Gaussian Process . . . . . . . . . . . . . . . . . 42.1 Computer Models . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 The Gaussian Process . . . . . . . . . . . . . . . . . . . . . . 42.3 Gaussian Process Emulation . . . . . . . . . . . . . . . . . . . 52.4 Estimating Theta . . . . . . . . . . . . . . . . . . . . . . . . . 72.5 Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.5.1 Correlation Function . . . . . . . . . . . . . . . . . . . 82.5.2 Prediction . . . . . . . . . . . . . . . . . . . . . . . . . 92.6 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.6.1 Gradient . . . . . . . . . . . . . . . . . . . . . . . . . . 11Chapter 3: Implementation . . . . . . . . . . . . . . . . . . . . . 133.1 Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133.2 Correlation Functions . . . . . . . . . . . . . . . . . . . . . . 153.2.1 Invariance . . . . . . . . . . . . . . . . . . . . . . . . . 183.3 Optimizing in R . . . . . . . . . . . . . . . . . . . . . . . . . 193.4 Assessing Fit . . . . . . . . . . . . . . . . . . . . . . . . . . . 20Chapter 4: Borehole Example . . . . . . . . . . . . . . . . . . . 22ivTABLE OF CONTENTSChapter 5: Numerical Testing . . . . . . . . . . . . . . . . . . . 295.1 Full Simulation Study . . . . . . . . . . . . . . . . . . . . . . 295.2 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . 295.3 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 305.4 Performance Analysis . . . . . . . . . . . . . . . . . . . . . . 325.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 34Chapter 6: Conclusions . . . . . . . . . . . . . . . . . . . . . . . 42Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47Appendix A: Test Functions . . . . . . . . . . . . . . . . . . . 47vList of TablesTable 3.1 Summary table specifying: parameterization, optimiza-tion method, and type of gradient used in the Gaussianprocess. . . . . . . . . . . . . . . . . . . . . . . . . . . 18Table 4.1 Log-likelihood values for the 10 optimizations for onefitting. . . . . . . . . . . . . . . . . . . . . . . . . . . 24viList of FiguresFigure 1.1 Root mean squared errors for 5 separate Gaussianprocess fittings for 11 different methods. . . . . . . . 2Figure 3.1 Three different samples for the initial input points forn = 5 input points with d = 2. . . . . . . . . . . . . . 14Figure 4.1 Log-likelihood values for the 5 different fittings of theBorehole test function with 11 parameter-optimizercombinations. . . . . . . . . . . . . . . . . . . . . . . 25Figure 4.2 Standardized root mean squared errors for the 5 dif-ferent fittings of the Borehole test function with 11parameter-optimizer combinations. . . . . . . . . . . 26Figure 4.3 A closer examination of the log-likelihood values ofthe two φ2 parameterizations using upper bounds of5 (left) and 100 (right). . . . . . . . . . . . . . . . . 28Figure 5.1 A boxplot of the log-likelihood values across all pa-rameterizations for the full test. . . . . . . . . . . . . 31Figure 5.2 Box plots of the logarithm (base ten) of both of theRMSE and AME values . . . . . . . . . . . . . . . . 33Figure 5.3 Empirical cumulative distribution functions for thelogarithm (base ten) of the RMSE and AME values. 35Figure 5.4 Empirical cumulative distribution functions for theperformance of each parameterization. . . . . . . . . 36Figure 5.5 Empirical cumulative distribution functions for thetop three parameterizations with respect to perfor-mance. . . . . . . . . . . . . . . . . . . . . . . . . . . 37Figure 5.6 Empirical cumulative distribution functions of RMSEand AME for the top four parameterizations with re-spect to errors. . . . . . . . . . . . . . . . . . . . . . 39viiLIST OF FIGURESFigure 5.7 Empirical cumulative distribution functions of RMSEand AME for the top thee parameterizations (withrespect to performance) where performance is greaterthan -8. . . . . . . . . . . . . . . . . . . . . . . . . . . 41viiiAcknowledgementsFirst and foremost I would like to thank to my supervisor, Dr. JasonLoeppky. Your inspiration, guidance, encouragement, and support havebeen invaluable for my success during this project.I would also like to thank my supervisory committee: Dr. Paramjit Gilland Dr. John Braun. I sincerely appreciate your support, and thank youfor dedicating your time to my work.My gratitude goes to Dr. Patricia Lasserre for her constant support andsuggestions, as well as to Dr. Shirin Golchi for her invaluable insight andhelp in the work presented here.I wish to say thank you to the many students and friends I have hadthe pleasure of knowing or working with during my studies. My lab matesespecially, provided me with encouragement, kind words, and coffee breaksto ensure that I would be able to power through debugging. Words cannotexpress how thankful I am to have had you by my side during this experience.Finally, a special thanks to Spencer Hunt, Trent Jensen, and RachelStone who were the first to encourage me to peruse graduate studies; withoutyou, this thesis would never have been started. For this, I would like to say“Thank you.”ixDedicationFor my friends and especially my family; thank you for your uncondi-tional love, encouragement, and support though my, at times tumultuous,journey.xChapter 1IntroductionMathematical models are used across the statistical, mathematical, andengineering sciences. Mathematical models are implemented as computercode when physical experiments are costly or impossible to implement, re-placing the physical experiments. The science being modeled by the mathe-matical code is often complex, requires a large number of variables and theycan be slow to run. Statistical models, such as a Gaussian process can beused to approximate the code [5, 20–22] when the code is expensive to run.Successfully implementing a Gaussian process surrogate model involvesthree major steps: choosing a design, specifying the form of the Gaussianprocess and then estimating the parameters in the Gaussian process. Tochoose a design that is appropriate for fitting a Gaussian process involvesselecting an appropriate sample size as well as the design itself.The exact choice of the form of the Gaussian process depends on thechoice of the regression function as well as the choice of the correlationfunction. This thesis focuses on different selections for the correlation func-tions and the resultant effect on the fittings of the Gaussian process. Figure1.1 shows error values from five separate runs of fitting the Gaussian processmodel (denoted as five separate lines) for eleven different Gaussian processcorrelation functions. In the figure, the correlation functions used for fittingthe Gaussian process produced different final results (shown by the differenterror values on the y-axis. These different results suggest a problem of in-variance; the different correlation functions should produce the same results,however, as shown in Figure 1.1, this is not the case. [3] consider efficientlynumerically optimizing the Gaussian process. They implement several al-gorithms for choosing the starting values and optimize the GP using oneparametrization and illustrate this on a few examples. In this thesis wegreatly extend this work to consider several parameterizations and a muchlarger class of test functions.This thesis focuses on determining which parameterization of the cor-relation function performs the best over a set of test functions and canbe denoted as the ‘go-to’ parameterization to fit a Gaussian process. Theremainder of the thesis is organized as follows, in Chapter 2 we formally1Chapter 1. Introductionllllllllllllllllllllllllllllllllllllllllllllllllll−4.0−3.5−3.0−2.5−2.0Bt Bt_ng Lam Ph1 Ph2 Rh Rh_ng Th Th_un Th1 ThBStandardized RMSE Figure 1.1: Root mean squared errors for 5 separate Gaussian process fittingsfor 11 different methods.2Chapter 1. Introductionintroduce the Gaussian process and outline the method of fitting a Gaus-sian process model. In Chapter 3 we delve into the three steps outlinedabove, choosing a design, discussing various parameterization options, andestimating the parameter values, while hinting at the problem of invariance.Chapter 4 introduces an example test function (the Borehole function) andthe testing strategies used for the full simulation study. Chapter 5 performsthe full test simulations and discusses the analysis and results. Finally, weconclude in Chapter 6 with recommendations.3Chapter 2The Gaussian ProcessComputer ModelsComputer models are mathematical representations of physical experi-mentations that are implemented using computer code. These techniques areused commonly throughout the physical sciences when physical experimen-tation is expensive or time-consuming. Observations are made by runninga computer model at a set of input points. This is typically referred toas a computer experiment. The computer code is deterministic - producingidentical output for different runs of a given input - though the output is un-known before running the model. Consider a black-box function η(x) withinput x = (x1, . . . , xd), where (for simplicity) x ∈ [0, 1]d. The black-boxfunction, η(x), is known in the way that we know the code takes an inputx and produces an output y(x), however, the complexity of the code essen-tially makes it an unknown function that requires modelling. Thus, we treatthe function η(x) as unknown. From a Bayesian standpoint, the uncertaintyof the output can be expressed using a stochastic process which is knownas a statistical emulator. We can use a Gaussian process (GP) to representthe unknown function [15]. After building a GP emulator using runs of thecomputer code, analyses can be made without making more computer coderuns. One such analysis is to predict output for a new set of locations notrun.The Gaussian ProcessAs stated above, we can use a Gaussian Process to represent the un-known function η(x). The Gaussian process can be formally defined asa collection of random variables, any finite number of which have a jointGaussian Distribution [17]. Given two vectors of data x and x∗, the meanfunction m(x) and the covariance function V (x,x∗) are defined asm(x) = E[f(x)], (2.1)42.3. Gaussian Process EmulationandV (x,x∗) = E[(f(x)−m(x))(f(x∗)−m(x∗))] (2.2)respectively.Gaussian Process EmulationIn this section, we will discuss the idea of the Gaussian process emula-tor. Given the unknown function η(x) with input x = (x1/ldots, xd) wherex ∈ [0, 1]d and corresponding output y = (y(x1), . . . , y(xd))T . The uncer-tainty of the function η(x) is described as a Gaussian process with meanfunction m(x) and covariance function V (x,x∗). Formally, if η(x) has aGaussian process distribution, then for every n = 1, 2, . . ., the joint distri-bution of η(x(1)), . . . , η(x(n)) is multivariate normal for all x(1), . . . ,x(n).The choice of m(x) can be any function of x, but the choice of V (x,x∗)must produce a covariance matrix with elements {x(i),x(j))} that satisfiesthe property of being nonnegative definite. Prior information of η(x) is rep-resented by a Gaussian process with mean and covariance functions, m(x)and V (y(x), y(x∗)) respectively. The distribution of η(x) isη(x)|β, σ2, θ ∼ GP (m(x), V (y(x), y(x∗))), (2.3)where the mean function m(x) is given bym(x) = F T (x)β, (2.4)where F T (x) = (F 1(x), . . . ,F p(x))T is a set of known basis functions rep-resenting the mean and β is a set of unknown regression coefficients. Con-sidering that η(x(i)) is unknown, it is feasible to assume that F (x) = 1 with1 being a vector of ones as it allows the covariance function to model all ofthe signal [4]. The covariance function V (y(x), y(x∗)) is given byV (y(x), y(x∗)) = σ2R (x,x∗|θ) , (2.5)where σ2 is an unknown scale parameter and R is a known correlation func-tion with a vector of unknown correlation parameters θ. The choice ofR(x,x∗|θ) should ensure that the covariance matrix is nonnegative definitefor any set of input points. The choice for the exact form of R will be dis-cussed in Section 2.5.1 and Section 3.2. At this point, we will consider thegeneral function R. Suppose that we have a set of n inputs x(1), . . . ,x(n),and the corresponding outputs y = (η(x(1)), . . . , η(x(n)))T . These are de-noted as the training data. According to Equation 2.3, the distribution of52.3. Gaussian Process Emulationthe outputs y is multivariate normal,y|β, σ2, θ ∼ N(F Tβ, σ2R), (2.6)where R is the correlation matrix comprised of elementsRij = R(x(i),x(j), |θ). (2.7)Using multivariate normal distribution conditioning techniques [7], it can beshown thatη(x)|β, σ2, θ,y ∼ GP(m∗(x), V ∗(y(x), y(x∗))), (2.8)wherem∗(x) = F Tβ +(rT (x)R−1(y − Fβ)) (2.9)andV ∗(x,x∗) = σ2(R(x,x∗|θ) + rT (x)R−1r(x∗)) (2.10)where r(x) =(R(x,x(1)), . . . , R(x,x(n)))T. Using a weak prior for (β, σ2),p(β, σ2) ∝ σ−2, and combining with Equation 2.6 using Bayes’ Theorem,the posterior for (β, σ2) is a normal inverse-gamma distribution,β|y, σ2,θ ∼ N(βˆ, σ2(F TR−1F )−1), (2.11)whereβˆ = (F TR−1F )−1F TR−1y, (2.12)andσ2|y,θ ∝ InvGam(n− q2,(n− q − 2)σˆ22), (2.13)whereσˆ2 =yT (R−1 −R−1F (F TR−1F )−1F TR−1)yn− q − 2 . (2.14)Again, following [7], we can then integrate β our of the product of Equation2.8 and Equation 2.11 to achieveη(x)|y, σ2 ∝ GP (m′(x), V ′(x,x∗)), (2.15)wherem′(x) = F T (x)βˆ + rT (x)R−1(y − F βˆ), (2.16)62.4. Estimating ThetaandV ′(x,x) =σ2(R(x,x∗|θ)− rT (x)R−1r(x∗)+ (F (x)− rT (x)R−1F )(F TR−1F )−1(F ∗T − rT (x∗)R−1F )T .(2.17)We can now integrate σ2 out of the product of Equation 2.13 and Equa-tion 2.15 to obtainη(x)|y,θ ∝ StudentProcess(n− q,m′(x), V ′′(x,x∗)), (2.18)whereV ′′(x,x∗) =σˆ2σ2V ′(x,x∗). (2.19)Estimating ThetaFollowing [5], we use an empirical Bayesian approach to find values forθ by maximizing the likelihood. The likelihood for the GP isL(β, σ2,θ|y) =(2pi)−n/2(σ2)−n/2|R|−1/2· exp{− 12σ2(y − Fβ)TR−1/2(y − Fβ)}.Taking the natural logarithm yields the log-Likelihood function [28], whichis,l(β, σ2,θ|y) =− n2log(2pi)− n2log(σ2)− 12log(|R|)− 12σ2(y − Fβ)TR−1/2(y − Fβ).Maximizing with respect to each of the parameters β and σ2 results in theestimates ofβˆ = (F TR−1F )−1F TR−1y (2.20)andσˆ2 =1n(y − F βˆ)TR−1(y − F βˆ). (2.21)We can use the above two substitutions for β and σ2 to provide the profilelikelihood function [28] that relies only on θ, which is contained within Rl(θ|y) = −n2log(2pi)− n2log(σˆ2)− 12log(|R|)− n2. (2.22)Finding the maximum of the likelihood function with respect to the param-eter, θ, requires numerically optimizing Equation 2.22.72.5. ModellingModellingCorrelation FunctionThe successfulness of the Gaussian process is dependent upon the choiceof the correlation function R. There are some general properties that wewould like to see characterized by the function. They are as follows:− Replicate observations should be perfectly correlated, so x(i) = x(j)then R(x(i),x(j)) = 1.− Points that are close to one another should have function values thatare highly correlated. That is, if ||x(i) − x(j)||2 = ∑dk=1(x(i)k − x(j)k )2is small, then R(x(i),x(j)) should be close to one.− Points that are far away from one another should have function valuesthat are basically unrelated (not correlated). That is, if ||x(i)−x(j)||2is large then R(x(i),x(j)) should be close to zero.There are many different correlation functions available for use. The mostcommon are the Mate´rn, power exponential, and the Gaussian correlationfunctions. The Gaussian correlation function is widely used as the correla-tion function of choice [18, 25, 27] and we choose to use it due to its char-acteristic of being indefinitely differentiable, and therefore is a very smoothfunction. The Gaussian correlation function for a d-dimensional functionhas the general form ofR(x(i),x(j)|θ)=d∏k=1exp{−θk(x(i)k − x(j)k)2}, (2.23)where θk ∈ [0,∞) for all k = 1, . . . , d. The values of the parameter θ controlthe activity of the function being modelled. A value of θk = 0 removesthe kth dimension from the correlation completely. Specifically, a value ofθk = 0 in the kth dimension corresponds to an overall value of 1 in theproduct for that dimension. This value of 1 results in no change, or effect,on the final value of R due to the kth dimension. This means that the kthdimension does not explain any of the variability of the function output,y. Oppositely, large values of θk correspond to a more active or fluctuatingfunction in the kth dimension. Some implementations of a GP use varyingforms of the parameter θ [19] and [21]. These various parameterizations willbe explored further in Section 3.2.82.5. ModellingPredictionConsider the case of predicting the output, m′(x) for a single pointx∗, where x∗ = x(i) for some i = 1, . . . , n. In this case, F ∗ = fT (x(i)),and rT =(R(xi,x(1), . . . , R(xi,x(n))), which is the ith row of the originalvariance-covariance matrix R. That is,r = R (0, . . . , 0, 1, 0, . . . , 0) ,where the 1 is in the ith position. We will denote this vector as ui as itis the ith unit vector, or ith row of the identity matrix I. r can then bedenotd as:r = Rui.Multiplying both sides of the above equation by the inverse of R we obtain,R−1r = ui.Take the transpose of each side to obtainrTR−1 = uTi .We can then substitute the above result (as well as the simplification ofF ∗ = fT (x(i)) into Equation 2.16:m∗ = F ∗β + (rTR−1(y − Fβ))= fT (x(i))β + (uTi (y − Fβ))= fT (x(i))β + y(i) − fT (x(i))β= y(i)obtaining the original output value y(i) as the prediction of the point x(i).Now let us consider the case of predicting m′ for the full set of n trainingpoints; x∗ = x. In this case, our original dataset x = x∗ which meansthat r = R. Recalling that R is the correlation matrix and is thereforesymmetric, we have RT = R. Finally, seeing as there are the same numberof points in x and x∗ then F ∗ = F . We can then calculate m′:m∗ = F ∗β + (rTR−1(y − Fβ))= Fβ + (RTR−1(y − Fβ))= Fβ + (RR−1(y − Fβ))= Fβ + (I(y − Fβ))= Fβ + y − Fβ= y92.6. OptimizationOptimizationFinding the maximum likelihood estimates requires maximizing Equa-tion 2.22 with respect to the unknown vector θ. There are several optimiza-tion methods that could be used to find the maximum likelihood estimated.These include genetic algorithms [14], mesh adaptive seeded algorithms [2],simplex algorithms (such as the Nelder-Mead simplex algorithm [6]), andgradient based methods [24]. Genetic algorithms will converge to the globalmaximum if enough iterations are performed, which means they can be veryslow. Both the mesh adaptive seeded algorithms and simplex algorithm donot require the gradient in order to perform the optimization. Typically,gradient-based algorithms are faster than non-gradient based algorithms,providing the gradient is easy to compute. In the case of the GP, the com-putation of the gradient is straightforward and the computational cost is lessthan computing the likelihood meaning that gradient based methods shouldbe ideal in this situation. In order to deal with the possibly bounded searchspace on θ we suggest using the bounded Limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS-B) method or its unbounded counterpart, theBroyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. The purpose of theL-BFGS-B algorithm is to minimize a d dimensional nonlinear function, f(x)wherel ≤ x ≤ uwhere l and u are vectors containing the lower and upper bounds on thevariable x respectively. Some, or all of the variables may not have boundsresulting in unconstrained optimization. At each iteration of the algorithman approximation to the Hessian matrix is updated. This matrix is usedto define a quadratic model of the function. Using a two-stage approach, asearch direction is then calculated; first the gradient projection method isused to identify active variables and then the quadratic model is approxi-mately minimized with respect to the remaining free variables. This searchdirection is defined as the vector leading from the currently position to theapproximate minimizing location. Finally, a line search is performed alongthe search direction [30]. The Broyden-Fletcher-Goldfarb-Shanno method isimplemented using the gradient, but will require a method to choose startingvalues for the parameter with which we are optimizing in respect to. It ispossible to implement BFGS without a gradient as the optim function willcompute a finite-difference approximation. [16] shows that fitting the Gaus-sian process without using the gradient works well, but we could make usethe knowledge of the gradient function which we compute below in Equation2.26 for our optimization.102.6. OptimizationGradientThe exact form of the analytical gradient will depend on the correlationfunction used, but the general form of the gradient can be calculated withoutit. To calculate the gradient we use the profile log-likelihood function fromEquation 2.22l(θ|βˆ, σˆ2,y) = −n2log(2pi)− n2log(σˆ2)− 12log(|R|)− n2(2.25)The gradient is then the derivative of this function with respect to theparameter of interest. In this derivation we will use the parameter θ` forillustrative purposes. The `th entry of the gradient vector, G, taken withrespect to θ` is calculated asG`,θ` =∂∂θ`[l(θ|βˆ, σˆ2,y)]=∂∂θ`[−n2log(2pi)− n2log(σˆ2)− 12log(|R|)− n2]= −n21σˆ2∂∂θ`[σˆ2]− 121|R|∂∂θ`[|R|]= −n21σˆ2∂∂θ`[1n(y − F βˆ)TR−1(y − F βˆ)]− 121|R|∂∂θ`[|R|]G`,θ` =12σˆ2(y − F βˆ)TR−1 ∂∂θ`[R]R−1(y − F βˆ)− 121|R| |R|tr(R−1∂∂θ`[R]).(2.26)The gradient is now solely in terms of ∂∂θ` [R]. For the case of the Gaussiancorrelation function introduced in Section 2.5.1, we can calculate ∂∂θ`R(x(i),x(j)|θ)by taking the derivative of Equation 2.23 with respect to θ`.∂∂θ`R(x(i),x(j)|θ)= exp{d∑k=1−θk(x(i)k − x(j)k)2}(−(x(i)` − x(j)` )2) .(2.27)Substituting the above into Equation 2.26 and repeating for each of thedimensions, ` = 1, . . . , d, we can achieve the full analytical gradient for use112.6. Optimizationin optimization. In the next chapter we discuss the full implementation ofthe GP model which involves selecting an initial design, specifying the exactparameterization of the correlation matrix, maximizing the likelihood andpredicting a set of new observations.12Chapter 3ImplementationIn practice, fitting a Gaussian process requires selecting the initial design,optimizing the log-likelihood function in Equation 2.22 with respect to theparameter of choice by estimating the parameter values, and then predictingvalues for a new set of input points.DesignThe problem of selecting the initial input values, x, is a common topic incomputer modelling. The approach we will take is to treat the input valuesthemselves as random variables. This method models the uncertainty of theinput values themselves.We will refer to the initial input data as the design of the experiment.Due to our predicting functions defined in Equations 2.16 and 2.17, it isimportant to create a design that is space-filling.The prediction equations are interpolators, as discussed in Section 2.5.2,which results in the prediction error at a given input site being a functionof its location (relative to the design points x). If we were to select a non-space-filling design, it may yield predictors that perform poorly in the areasthat are sparsely observed.There are many different approaches that we could use in the designprocess; we could select points using a line, a simple random sample, or alatin hypercube sample, just to name a few. The random latin hypercubesampling method has been used by many as the approach in designing ex-periments, [13, 20, 23, 28] and it is the method that we will use to createthe initial design for our study. First introduced by [12], the latin hyper-cube sampling method allows all portions of the sample space for a singledimension, S, to be represented. The method starts by dividing the sam-ple space into n equal-width intervals, Si ∈ S, where n is the sample size,and i = 1 . . . n. For simplicity, we will continue with a sample space ofS = [0, 1]. To form the latin hypercube sample, we will select points fromv, where v ={Ri+εin−1 |Ri = i− 1, i = 1 . . . n}. Using the points in v, we willgive them a random ordering forming a vector. We continue this method133.1. Designfor each of the d dimensions we can then match the vectors as columns ina matrix to form the final design. This results in an n-by-d matrix whereeach column is some permutation of the elements contained in v.There are different methods to set the values of εi. We could set aspecific value such as εi = 0, or εi = 1/2, or we could use a distributionsuch as εi ∼ UNIF (0, 1). Selecting a value of εi = 0 forces the points tothe boundary of the space, whereas ε = 1/2 places them at the midpoint ofeach space.lllll0.000.250.500.751.000.00 0.25 0.50 0.75 1.00(a)lllll0.000.250.500.751.000.00 0.25 0.50 0.75 1.00(b)lllll0.000.250.500.751.000.00 0.25 0.50 0.75 1.00(c)Figure 3.1: Three different samples for the initial input points for n = 5input points with d = 2.Three different realizations of a random latin hypercube sample, with143.2. Correlation Functionsd = 2 and n = 5, are contained in Figure 3.1. All of Figure 3.1 (a), (b),and (c) are valid latin hypercube samples with each subsection Si (eachvisual ‘row’ or ‘column’ division of the graph) containing a point. Figure3.1 (a) and (b) are two different realizations of the latin hypercube samplingmethod using εi = 1/2 and (c) uses εi = 0.In terms of being a space filling design, the first graph, (a), would beconsidered a poor design as it leaves two large areas of the sample space notsampled and the same argument can be made for (b) and the area along thediagonal. As well, the use of εi = 1/2 in (a) and (b) does not allow for theboundary points to be sampled which can allow for boundary effects to notbe caught when fitting the GP. Graph (c) uses εi = 0 and forces points tobe along the outer boundaries (0 and 1) which will allow boundary effectsto be captured in the model created by the Gaussian process. For the thisthesis we place points on the boundary by choosing ε = 0.The maximin latin hypercube sample [9], maximizes the distance be-tween the points in the design, can be used to ensure that the resulting latinhypercube design is space filling. However, in order to make the problem ofsolving the Gaussian process for an experiment more difficult, we limit ourdesign to a random latin hypercube sample. This will ensure that resultsgeneralize to designs that are space filling.Correlation FunctionsThe major focus of this thesis is efficiently optimizing the likelihood,which we will show is highly dependent on the choice of the parameterizationof the likelihood. Revisiting the correlation function R introduced in Section2.5.1 we will explore some other forms (parameterizations) found in theliterature. We will denote these different methods listed below as differentparameterizations of the correlation function. We wish to investigate thesedifferent parameterizations to evaluate if a best all around parameterizationcan be found to fit a Gaussian process.Recall the correlation function (parameterization) in Equation 2.23 pre-sented in Section 2.5.1. We will denote this parameterization as theta, orθ.The next parameterization will be called theta2, or θ2, and is a variationof the theta parameterization. The correlation function isR(x(i),x(j)|θ)=d∏k=1exp{−θ2k(x(i)k − x(j)k)2}. (3.1)153.2. Correlation FunctionsThe domain of θ2 is the same as for the θ case, θk ∈ [0,∞),∀k.The θ2 casespreads out the small values close to zero. The main difference between thetwo comes to light in the derivative ∂∂θ` [R]:∂∂θ`R(x(i),x(j)|θ)= exp{d∑k=1−θ2k(x(i)k − x(j)k)2}(−2(x(i)` − x(j)` )2 θ`) .(3.2)The next parameterization is a form of the Gaussian correlation functionthat is used in the DiceKriging package in the program R [19]. We willdenote it as phi2, or φ2. It varies from the previous parameterizationsby introducing the parameter in the denominator of the exponent. Thecorrelation function for phi2 isR(x(i),x(j)|φ)=d∏k=1exp−(x(i)k − x(j)k√2φk)2 . (3.3)Due to the parameter being in the denominator, the domain is now φk ∈(0,∞), ∀k. The partial derivative for the gradient is calculated to be∂∂φ`R(x(i),x(j)|φ)= exp{d∑k=1−12φ2k(x(i)k − x(j)k)2}(φ−3`(x(i)` − x(j)`)2).(3.4)The next parameter beta (β) was introduced by [16] and is the parame-terization used in the study done by [3].It takes a different approach to theGaussian correlation form raising the parameter as an exponent:R(x(i),x(j)|β)=d∏k=1exp{−10βk(x(i)k − x(j)k)2}. (3.5)The domain for β is βk ∈ (−∞,∞),∀k, and the derivative can be calculatedas∂∂β`R(x(i),x(j)|β)= exp{d∑k=1−10βk}(− ln(10)(x(i)` − x(j)`)210β`).(3.6)The next parameterization was suggested by [8] and [10] and it placesthe parameter in the base of the equation rather than the exponent. Thedomain of this parameter, labeled rho (ρ), is restricted to ρk ∈ (0, 1], ∀k and163.2. Correlation Functionsthe correlation function isR(x(i),x(j)|ρ)=d∏k=1ρ4(x(i)k −x(j)k)2k . (3.7)The derivative can be shown to be∂∂ρ`R(x(i),x(j)|ρ)=∏k 6=`ρ4(x(i)k −x(j)k)2k·(4(x(i)` − x(j)`)2ρ4(x(i)` −x(j)`)2−1`),(3.8)which can be simplified to∂∂ρ`R(x(i),x(j)|ρ)=(d∏k=1ρ4(x(i)k −x(j)k)2k)·4(x(i)` − x(j)`)2ρ` .(3.9)The final parameterization will be labeled as lambda (λ). It is a rear-rangement of the ρ parameterization with relationλk = log(ρk)− log(1− ρk),orρk =exp{λk}1 + exp{λk} .Using the above relation, the correlation function isR(x(i),x(j)|λ)=d∏k=1(exp{λk}1 + exp{λk})4(x(i)k −x(j)k )2, (3.10)with corresponding derivativeddλ`R(x(i),x(j)|λ)= d∏k=1(exp{λk}1 + exp{λk})4(x(i)k −x(j)k )2· 4(x(i)` − x(j)`)2 1− 2 exp{λ`}1 + exp{λ`} .(3.11)173.2. Correlation FunctionsThe five parameterizations presented above as well as the one first intro-duced in Section 2.5.1 will be used in fitting the model for the Gaussian pro-cess. As specified in Section 2.6, two methods of optimization (L-BFGS-Band BFGS) along with two versions of the gradient (analytical and numer-ical) will be used. Table 3.1 is a summary table specifying the parameter-ization for the correlation function, label used to differentiate between fits,the optimization method used, the bounds used (in the case of the boundedoptimizer), and the gradient used in fitting the model through the Gaussianprocess.Table 3.1: Summary table specifying: parameterization, optimizationmethod, and type of gradient used in the Gaussian process.Parameter Label Optimizer Bounds Gradientβ Bt BFGS analyticalβ Bt ng BFGS numericalλ Lam BFGS analyticalφ2 Ph1 L-BFGS-B (0, 5] analyticalφ2 Ph2 L-BFGS-B (0, 100] analyticalρ Rh L-BFGS-B (0, 1] analyticalρ Rh ng L-BFGS-B (0, 1] numericalθ Th1 L-BFGS-B [0,∞) analyticalθ2 Th L-BFGS-B [0,∞) analyticalθ2 Th un BFGS analyticalInvarianceLikelihood theory shows that maximum likelihood estimates are invari-ant for one-to-one transformations of the parameters [29]. That is to say, ifθˆ is the maximum likelihood estimator (MLE) of θ, for any function τ(θ),the MLE of τ(θ) is τ(θˆ). To show this, consider η = τ(θ) and the definedinduced likelihood function, L∗, asL∗(η|x) = sup{θ:τ(θ)=η}L(θ|x).Let ηˆ denote the value that maximizes L∗(η|x). We want to show thatL∗(ηˆ|x) = L∗[τ(θˆ)|x]. As previously stated, the maxima of L and L∗ coin-183.3. Optimizing in Rcide, and we haveL∗(ηˆ|x) = supηsup{θ:τ(θ)=η}L(θ|x)= supθL(θ|x)= L(θˆ|x)where the second equality follows because the iterated maximization is equalto the unconditional maximization over θ which is attained at θˆ, by defini-tion. Further we can showL(θˆ|x) = sup {θ : τ(θ) = τ(θˆ)}L(θ|x)= L∗[τ(θˆ)|x].Using the two above equalities we can show that L∗(ηˆ|x) = L∗(τ(θˆ)|x), andthat τ(θˆ) is the maximum likelihood estimator of τ(θ).A problem arises, with respect to invariance, within the φ2 parameteri-zation. Consider the case in the θ parameterization when θk = 0; a value ofθk = 0 corresponds to a value of φk =∞ which is numerically unattainableand results in a lack of invariance. Numerically optimizing φ requires a limitto be placed for the upper bound which renders the situation non-invariant.Forcing φk away from ∞ brings the parameter from non-active to activestatus and can result in overfitting. Different limits on the upper bound canresult not only in different maximum likelihood estimates for φˆ, but alsodiffering values for the corresponding log-likelihood (i.e. differing results).The remainder of this thesis is focused on effectively numerically optimizingthe likelihood function L(θ) and in turn investigating the effects due to theboundary.Optimizing in RUsing the parameterization, optimization, and gradient pairings pro-vided in Table 3.1, we use the statistical programming language R to per-form the optimization. The initial input points for x, or the design of theexperiment, are created using the random latin hypercube method as de-scribed in Section 3.1. As mentioned earlier, the latin hypercube sampleallows for a space filling sample to be chosen. We require the design to bespace filling so that any effects or trends due to a point being a boundarypoint will be captured in the model fit by the Gaussian process.To numerically optimize the likelihood function itself, the R functionoptim() is used. The optim() function will minimize a provided function193.4. Assessing Fitrather than maximize it, so the log-likelihood and gradient functions aremodified to reflect this. For the remainder of this thesis, we will talk ofminimizing the log-likelihood function where a minimal value translates tothe optimal value.The optimization function requires a set of starting values for the pa-rameters. We use a warm start method in order to provide these values.A total of 15 · d points are generated for each of θ2 and φ2 using the ran-dom latin hypercube method used to generate the input points x. For θ2,the points were generated to be within [0,min(0.4 · d, 1.6)], and for φ2 thepoints were within [10× 10−10,min(0.4 · d, 1.9)]. We then calculate the log-likelihood value for each of the 30 · d points. The twenty points that yieldthe smallest log-likelihood values, regardless of initial parameterization (θ2or φ2), are then used as the twenty starting values for the optim() functionfor all parameterization fittings.Running the optim() function twenty times will result in the output ofthe optimal log-likelihood values and their corresponding optimal parame-ters. The best (minimum) log-likelihood value out of these twenty is denotedas the winner of this run and it, along with its corresponding parameter val-ues, are returned as the final result of this run.Assessing FitSince this thesis is devoted to assessing the ‘best’ parametrization of thelikelihood it is important to determine how we define the term best. Theone obvious choice of best is simply based on the finding the largest valueof the likelihood. However, one could also considered measures based onthe predictive performance of the fitted model. Since all of the example arebased on toy problems which mimic real problems we can obtain a large setof points for validation. Given a set of m new locations not used for fittingwe can run the code and obtain outputs y(x(j)). The two obvious methodsfor validation are the Root mean squared error and the maximum error. Inorder to compare different function we consider the standardized values ofthe measures.The standardized root mean squared error is given asRMSE =√∑mi=1(y(x(j))− yˆ(x(j))2∑mi=1(y(x(j))− y¯t)2,where y¯t is the mean of the training data. The standardized max error is203.4. Assessing Fitdefined asAME =maxi=1,...,m∣∣y(x(j))− yˆ(x(j))∣∣maxi=1,...,m∣∣y(x(j))− y¯t∣∣ .If each of the methods is equivalent then each of these measures should allbe the same. If the methods are not different then smaller values of thelikelihood are worse and larger values of the AME and RMSE are worse.We will use each of these three measures together to assess the ‘best’ pa-rameterization.21Chapter 4Borehole ExampleConsider the Borehole function which calculates the flow rate through aborehole using eight input variables (d = 8). The function isf(x) =2pix3(x4 − x6)ln(x2/x1)(1 + 2x7x3ln(x2/x1)x21x8)+ x3x5.The variables and their appropriate ranges are given below:− radius of the borehole (m) – x1 ∈ [0.05, 0.15]− radius of influence (m) – x2 ∈ [100, 50 000]− transmissivity of upper aquifer (m2/yr) – x3 ∈ [63 070, 115 600]− potentiometric head of upper aquifer (m) – x4 ∈ [990, 1 110]− transmissivity of lower aquifer (m2/yr) – x5 ∈ [63.1, 116]− potentiometric head of lower aquifer (m) – x6 ∈ [700, 820]− length of borehole (m) – x7 ∈ [1 120, 1 680]− hydraulic conductivity of borehole (m/yr) – x8 ∈ [9 855, 12 045].This function is considered easy to fit using a Gaussian process. Usingthe procedure outlined in Chapter 3, we will fit the model for the boreholefunction using the Gaussian process. We start by selecting the design ofthe initial points by using the random latin-hypercube sampling design ofsize 7 · d = 56 points. Now that we have the initial input points, we can fitthe model using the Gaussian process. We will use the Gaussian correlationfunction parameterizations as outlined in Sections 2.5.1 and 3.2. A totalof eleven different models are to be created with the parameter-optimizer-gradient pairs outlined in Table 3.1. These different methods all used thesame starting design, of the 7 · d initial points.The starting values of the parameters are chosen using a warm startmethod. Initially, 15 · d starting values for the θ2 parameterization are22Chapter 4. Borehole Examplegenerated using a random latin hypercube sample in [0,min(0.4 · d, 1.6)].Another 15 · d starting points are generated for the φ2 parameterization onthe interval [10−10,min(0.4 ·d, 1.9)]. The log-likelihood values are calculatedfor both sets of 15 · d points. These 30 · d log-likelihood values are sortedand the twenty smallest values are chosen to be the starting values of thewarm start. The parameter values corresponding to the winning fits arethen translated to be in terms of both θ2 and φ2 using the relations:θ =1√2φ2,andφ =1√2θ2.Now we can use the following relations to translate the initial startingpoints to be in terms of all other parameterizationsθ =√θ2,β = log10(θ2),ρ = exp{−θ24},λ =−θ24− log(1− exp{−θ24}).Each parameterization will then start on even ground by starting withequivalent parameter values. The warm start method allows us to use thecomputationally cheaper method (compared to optimization) of calculatingthe log-likelihood values to find parameter locations of lower log-likelihoodvalues from which to start our optimization. The warm start saves in termsof the number of function calls during the optimization as well as compu-tational time. The twenty sets of parameter points are used as the startingparameter values for a run though the optimizer. The output will consist oftwenty sets of optimal parameter points and the corresponding optimal log-likelihood values. The final result of the run, deemed the optimal solutionfor this fitting, is the minimum log-likelihood value of the twenty values andits corresponding parameter values.Table 4.1 contains the optimal log-likelihood values for a single run offitting the model. Multiple methods (β, β(no gradient), θ2, θ2(unbounded),and θ2(β starting points)), were able to achieve the same result (to threedecimal places). Others (specifically λ, φ2(upper bound = 5), ρ, and θ)23Chapter 4. Borehole ExampleTable 4.1: Log-likelihood values for the 10 optimizations for one fitting.Fit Method Log-likelihood RMSE Max error1 Bt 125.016 0.393 6.5142 Bt ng 125.016 0.393 6.5143 Lam 164.347 2.866 15.8584 Ph1 161.518 3.458 20.9675 Ph2 130.226 0.500 7.7876 Rh 214.874 36.496 67.1467 Rh ng 125.984 0.486 6.6698 Th 125.016 0.393 6.5149 Th un 125.016 0.393 6.51410 Th1 136.746 1.019 10.14411 ThB 125.016 0.393 6.514performed poorly, achieving log-likelihood values greater much greater thanoptimal value achieved, 125.016. In this small example we will repeated theabove process four more times to get a total of five different results.We repeated the process by starting with generating a new set of inputpoints, x, for each of the replicates as well as a new set of starting pointsfor the parameters before running them through the optimization processand getting the outputs. Figure 4.1 shows the log-likelihood value resultsfor each of the five simulations. The eleven parameterizations are along thex-axis with labels corresponding to Method column in Table 4.1. Each ofthe simulations is represented by a line. The first simulation ran, the onewhere the values are shown in Table 4.1, is shown as the red line along withthe four replicates. Figure 4.2 shows the standardized root mean squarederrors from fitting the five simulations (the individual root mean squarederrors divided by the root mean squared error from the trivial predictor,y. Looking at both Figure 4.1 and Figure 4.2 we can start to see a pat-tern across the parameterizations. The parameterizations that performedwell on the first simulation (smaller log-likelihood and root mean squarederror values) as seen by the red line, also performed well on the other fourreplicates. Similarly, the parameterizations that performed poorly with highlog-likelihood and root mean square error values on the first simulation alsoperformed poorly on the following replicates.The drastic differences in the final log-likelihood values achieved by thedifferent parameterizations hint at the problem of invariance outlined inSection 3.2.1. As stated, likelihood theory shows that maximum likelihood24Chapter 4. Borehole Examplelllllllllllllllllllllllllllllllllllllllllllllllllllllll125150175200225Bt Bt_ng Lam Ph1 Ph2 Rh Rh_ng Th Th_un Th1 ThBLog−likelihoodFigure 4.1: Log-likelihood values for the 5 different fittings of the Boreholetest function with 11 parameter-optimizer combinations.25Chapter 4. Borehole Examplellllllllllllllllllllllllllllllllllllllllllllllllll−4.0−3.5−3.0−2.5−2.0Bt Bt_ng Lam Ph1 Ph2 Rh Rh_ng Th Th_un Th1 ThBStandardized RMSE Figure 4.2: Standardized root mean squared errors for the 5 different fittingsof the Borehole test function with 11 parameter-optimizer combinations.26Chapter 4. Borehole Exampleestimates are invariant for one-to-one transformations. This would meanthat for the transformations of the theta parameter, we expect to see thesame optimal log-likelihood values produced from optimizing the differentparameterizations. As can be seen in Figure 4.1, this is not the case. Forthe two fits using the phi parameterization there is a difference due to thebounds imposed. In the original parameterization, theta, the bounds wereθ ∈ [0,∞). Looking back at the translation between the φ and θ2 parame-terizations we haveθ =1√2φ2.The lower bound value of θ = 0 corresponds to a value of ∞ for φ, and theupper bound value of ∞ for θ to a value of 0 for φ. A value of φ = 0 is notpossible due to the parameter being in the denominator of the exponent andtherefore a restriction is imposed upon the domain of φ resulting in a domainof (0,∞). This non-equivalence in the domains results in the violation ofthe assumptions of invariance.Taking a closer look at the two φ parameterizations, the only differencein their implementation is the value of the upper bounds used. Figure 4.3shows the log-likelihood values for the five different fits for just the two φparameterizations. The fit using upper bound equal to 5 is denoted as avalue of 1 on the x-axis and upper bound equal to 100 is denoted as 2.We can see that there is clearly a difference between the two log-likelihoodresults. Looking closer at the values of the upper bounds, consider the caseof φ = 5. In the θ2 parameterization, this is equivalent to θ = 1√2·25 asa lower bound. The case of θk 6= 0, ∀k = 1, . . . , d forces activity in eachof the d dimensions. Even using the larger value of φ = 100 results in alower bound to nearly 0. This allows for inactive dimensions to becomeless active than the upper bound value of 5, but there would still be forcedactivity. Comparing these two cases to the default upper bound used by theR package DiceKriging, φ = 2, and we can use the same relation to showthat the DiceKriging scenario leads to a higher lower bound in the θ2 caseand therefore results in more activity than the two bounds that we use here.27Chapter 4. Borehole Examplellllllllll130140150160170Ph1 Ph2Log−likelihoodFigure 4.3: A closer examination of the log-likelihood values of the two φ2parameterizations using upper bounds of 5 (left) and 100 (right).28Chapter 5Numerical TestingFull Simulation StudyThe goal of this thesis is to determine which parameterization worksthe best across a set of various test problems. It is our wish to provideone parameterization to be deemed the ‘go to’ parameterization for anyonewishing to fit a Gaussian process. In order to determine which methodperforms the best we will need to test the different parameterizations on awide range of test problems. We use fourteen different functions of varyingdimensions (d = 2 to d = 10) for the full simulation. The full list of testfunctions can be found in Appendix ??.A total of eight hundred and forty tests are performed. A single test con-sists of fitting the model for each of the parameterizations using the methodoutlined in Chapter 3. Three different sample sizes are used for the creationof the initial design, x: 7·d, 10·d, and 15·d. [11] provide evidence supportingthe initial design choice of 10 · d as being an effective design. Overall, thereare a total of forty-two different sample size and function combinations. Werepeat each of these forty-two tests using the replication procedure used toget the four additional replicates done in Chapter 4. We repeat to get atotal of twenty replicates of each sample size-function combination. For atotal of eight hundred and forty tests. We treat each of the simulationsand replicates as separate problems with each needing analysis to determinewhich parameterization yields the best results for that problem.Numerical ResultsAs we to answer the question of which optimizer is the best using ourlarge set of test problems, the hope is that the best solver on our set oftest problems will also be the best solver on a set of currently untestedproblems. Performing the full simulation outlined in Section 5.1 results insolutions for each of the eight hundred and forty tests. The solutions con-sist of the optimal log-likelihood values for each of the parameterizationsas well as their corresponding parameter values. The obvious choice is to295.3. Error Analysisbenchmark using the log-likelihood value itself to determine which solutionis the optimal solution. The minimum log-likelihood value out of the elevendifferent parameterization options is the optimal solution for that run. How-ever, using the log-likelihood values for all of the eight hundred and fortytests to determine which parameterization performs the best is not feasibleas there is no way to make direct comparisons between the parameteriza-tions. Figure 5.1 shows the distributions of the log-likelihood values foreach of the parameterizations. The solutions come from the runs of fitting aGaussian process for different functions. These different functions’ optimallog-likelihood values may be drastically different from one another renderinga between-parameterization comparison uninformative. The graph does notallow us to link the log-likelihood values to the functions and test cases thatthey belong to, so it is not possible to determine if the log-likelihood valueachieved was good or not.Error AnalysisThere are many different tools that are available for analysis of statisticalmodels. Briefly mentioned in Chapter 4 two of the techniques that we willuse are the root mean squared error, and the absolute maximum error.Recall the RMSE for a set of m prediction pointsRMSE =m∑i=1√(yi − yˆi)2,where yi is the true value of the test function at input xi and yˆi is thepredicted value calculated using the model produced from the Gaussianprocess calculated using xi. Similarly the absolute maximum error AME iscalculated for the prediction points asAME = maxi=1,...,m|yi − yˆi|,with the same definitions of yi and yˆi as above. We standardize the errorsby dividing by the errors obtained by the trivial predictor of yˆi = y¯, where y¯is the mean of the responses originally used in fitting the Gaussian process.Using these standardized measure of errors, any value greater than 1 isconsidered to be extremely ill fitting as they perform worse than the trivialpredictor. As well, any predictors with a value greater than 0.5 are deemedpoor fitting. Finally, we take the logarithm (base 10) to spread the dataoutwards because the original values are small and close to zero. Now, an305.3. Error Analysisllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll−1500−1000−5000500Bt Bt_ng Lm Ph1 Ph2 Rh Rh_ng Th Th_un Th1 ThBfitlog−likelihoodFigure 5.1: A boxplot of the log-likelihood values across all parameteriza-tions for the full test.315.4. Performance Analysisextremely ill fitting model is one that produces an error value above 0, andan poor fit is any fit with an error above log10(0.5) = −0.301 .Both the RMSE and AME create one single measure of error using all ofthe predicted values from the fitting as well as the true function values. Forboth measures, a smaller value represents a better solution; the smaller thevalue the smaller the difference between the predicted value and true func-tion value and therefore the better the fit. Figure 5.2 shows two side-by-sidebox plots for the RMSE and AME values for the full set of simulation testcases. Looking at the RMSE box plot, the majority of the values are belowzero, and there is no distinct difference between the various fits. The AMEbox plot is similar with a smaller spread, and some values greater than zero,but no real difference between the distributions of the parameterizations.Using the box plots, it is difficult to determine a difference between the dis-tributions of the different parameterizations and therefore determine whichmethod yields the best results.Performance AnalysisThe log-likelihood is the most important tool to be used in evaluating theperformance of a model, but the error terms are also an important approachin looking for discrepancies. One of the main purposes of fitting a Gaussianprocess is to gain a model from which we can predict for new points. Thediscrepancies occur when we have a fit-produced model that has the bestlog-likelihood value, but it might produce predictions that yield large errorsand therefore be a bad fit prediction-wise.[1] argue for the use of a simplegraphical method for analyzing statistical methods. The graph they sug-gest, the Empirical-Cumulative-Distribution-Function (ECDF), yields someadvantages in assessing the performance of the different methods presented.The main advantage of the ECDF (when comparing to other commonlyused methods such as the box plot) is that it displays every data point whilealso allowing for easy interpretation. As well, the ECDF allows for multiplemethods to be displayed simultaneously and for direct comparisons betweenthe methods to be made easily. Figure 5.3 shows the ECDF’s of both theRMSE’s and the AME’s for each of the parameterizations.Compared to the boxplot representation in Figure 5.2, in the ECDFwe can see distinct differences in the distributions of the parameterizations.There is a set of parameterizations that are clearly set away from the othersconcentrated on the left side of the graph. The ECDF highlight the func-tions that appear on the left side of the graphs, and which have the largestconcentration of small error values and therefore the largest concentration325.4. Performance Analysisllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll−10−50Bt Bt_ng Lm Ph1 Ph2 Rh Rh_ng Th Th_un Th1 ThBfitlog10(RMSE)lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll−6−4−20Bt Bt_ng Lm Ph1 Ph2 Rh Rh_ng Th Th_un Th1 ThBfitlog10(AME)Figure 5.2: Box plots of the logarithm (base ten) of both of the RMSE andAME values .335.5. Results and Discussionof models (or fits) that perform well. The same evaluation can be performedon the AME graph.As well as using the different error values, we can create a metric fromthe log-likelihood values, Lp, to assess the performance of the fits of theGaussian process. We define the performance using the following definitionof Lp:Lp =(`p − `∗)|`∗| ,where `p is the log-likelihood solution found by the solver, and `∗ is theknown, optimal, log-likelihood value. The problem with the statistic Lp isthat the true solution value, `∗ , is not known and instead we are requiredto use the best (minimum) log-likelihood value found by all of the possiblesolvers for that run of the problem. We will denote this new measure asthe performance of the parameterization. This performance metric allowsus to use information from the log-likelihood values without the problem ofdifferent optimal values for the different test functions. Small values of Lpindicate that the parameterization achieved the optimal log-likelihood valuefor that run or a value close to the optimal value found.Allowing visualization of each of the eleven parameterizations on onegraph, the ECDF allows for a quick and reliable decision to be made regard-ing which method(s) may be outperforming others on the large set of testproblems. Figure 5.4 shows the ECDF for the performance of the elevenparameterizations across all test problems. As with the error ECDF’s,the methods performing well are those that are concentrated on the leftside of the graph. Figure 5.4 shows that there are three methods thathave an ECDF curve outperforming the others. These methods are β,θ2(unbounded), and θ2β. Figure 5.5 ignores the remaining parameteriza-tions and focuses on just the top three just listed. Looking at the line forθ2(unbounded) in Figure 5.5, we can see that it starts at a value of ap-proximately 0.38 on the y-axis. This means that in (roughly) 38% of thetest cases, the θ2(unbounded) parameterization achieved the optimal (theminimum) log-likelihood value. Acknowledging that multiple parameteriza-tions can achieve the optimal log-likelihood value, we can also note that θ2βachieved the optimal value approximately 35% of the test cases, and β, 25%.Results and DiscussionWe have explored three different ways to analyze the large set of data pro-duced by the full test simulation: the root mean squared error, absolute max-imum error, and performance metric (calculated using the log-likelihood).345.5. Results and Discussion0.000.250.500.751.00−10 −5 0log10(RMSE)0.000.250.500.751.00−6 −4 −2 0log10(AME)fitBtBt_ngLmPh1Ph2RhRh_ngThTh_unTh1ThBFigure 5.3: Empirical cumulative distribution functions for the logarithm(base ten) of the RMSE and AME values.355.5. Results and Discussion0.000.250.500.751.00−15 −10 −5 0log10(performance)fitBtBt_ngLmPh1Ph2RhRh_ngThTh_unTh1ThBFigure 5.4: Empirical cumulative distribution functions for the performanceof each parameterization.365.5. Results and Discussion0.000.250.500.751.00−15 −10 −5 0log10(performance)fitBtTh_unThBFigure 5.5: Empirical cumulative distribution functions for the top threeparameterizations with respect to performance.375.5. Results and DiscussionLooking back at the empirical cumulative distribution function graphs of theRMSE and AME values in Figure 5.3, there are four parameterizations thatfall to the left of the other seven and therefore have the highest proportionof small error values. These four are β, β(no gradient), θ2(unbounded), andθ2β. The parameterizations that perform the worst, by falling on the rightside of the graph, are ρ, φ(upper bound 5), and λ. These poor performingparameterizations are the parameterizations with the highest concentrationsof large RMSE and AME values when compared with the rest.In Markov chain Monte Carlo settings, [8] and [10] have had much successwith the rho parameterization. However this is not true when numericallyoptimizing the likelihood. Rho performs quite poorly with the smallest (andbest) error values being greater than a large portion of the error values in thefour best parameterizations. It is interesting to revisit the φ parameteriza-tions and note that the only difference between the two φ parameterizationsis the upper bound values of 5 and 100, but the curves in the error ECDF’s,and thus the solutions achieved, are not the same.Figure 5.6 shows the ECDF’s of just the top four parameterizations: β,β(no gradient), θ2(unbounded), and θ2β. With the exception of the firstinitial increase from 0 to 0.125, the graphs are identical to one another.There is no apparent difference between the distributions of the errors forthe four parameterizations.With no obvious difference between the top four parameterizations interms of errors, we need to delve deeper into the performance results basedon the log-likelihood values in order to pick the top parameterization. Look-ing back at the best parameterizations (β, β(no gradient), θ2(unbounded),and θ2β), we can see that out of the top four, three are the same as the topparameterizations found by analyzing Figure 5.5. The θ2(unbounded) andθ2β parameterizations have very similar curves and perform better than theβ parameterization until roughly the 0.87 and 0.93 marks on the y-axis re-spectively. The problematic area is the points where the β parameterizationperforms better than the theta(unbounded) and θ2β parameterizations (per-formance > -8). With the β curve being to the left of both of the other two,this means that the θ2(unbounded) and θ2β parameterizations have a largerconcentration of points at the higher end of performance values. Recall thatperformance is defined asLp =(`p − `∗)|`∗| ,and the larger the value, the further away the solution is from the bestsolution found. That is to say that the θ2(unbounded) and θ2β parameteri-385.5. Results and Discussion0.000.250.500.751.00−10 −5 0log10(RMSE)0.000.250.500.751.00−6 −4 −2 0log10(AME)fitBtBt_ngTh_unThBFigure 5.6: Empirical cumulative distribution functions of RMSE and AMEfor the top four parameterizations with respect to errors.395.5. Results and Discussionzations have a higher concentration of problems where they found solutionsthat were further away from the best solution found. In order to explorethese values more we will look at only the observations from the three bestparameterizations (β, θ2(unbounded), and θ2β) where performance is greaterthan -8. Figure 5.7 shows the empirical cumulative distribution functionsfor the RMSE and AME values of the three parameterizations for the sub-set of data that fit the criterion listed above. Looking at the ECDF’s forthe errors, we can see that the β curve falls to the right (with higher val-ues) of the other two. The θ2(unbounded) and θ2β parameterizations havesimilar curves with the exception of the area between 0.5 and 0.75 on they-axis. In that section θ2(unbounded) has the upper hand. Overall, the twoparameterizations have a similar performance, error wise.The β parameterization is out of contention for the title of best parame-terization because in this section of problematic performance values, it mayhave a better performance value than the other two parameterizations, butthe errors produced are much larger than the other two parameterizations.This shows that the log-likelihood value might have been smaller, but themodel itself performs poorly and produces estimates that are far from thetrue function values. Overall, we can say that the θ2(unbounded) and θ2β pa-rameterizations have the best performance, error wise, for the problematicsection of high performance values.We are aiming to pick the parameterization that performs the best overall of the problems presented in the full simulation. Currently, we are ableto narrow down the list of contenders to two parameterizations whose per-formances are relatively similar in terms of both the log-likelihood values(performance) as well as the errors (root-mean-squared-error and absolutemaximum error). These parameterizations are the θ2(unbounded) and theθ2β parameterizations. The main differences (and downside) of the θ2β param-eterization is the way in which the starting points of the parameter valuesare generated. The optimization was started at points equivalent to theparameter outputs of the finished optimization of the β run. Time wisethe θ2β run requires that the β parameterization run first and then we canperform the θ2β parameterization, where as, the θ2(unbounded) parameter-ization only requires one run in order to perform the optimization and fitthe full Gaussian process and achieve similar results. Combining what weknow with the performance analysis and error analysis, with the amount oftime required to run the parameterizations, the parameterization of choicefor the best parameterization is the θ2(unbounded) parameterization.405.5. Results and Discussion0.000.250.500.751.00−15 −10 −5 0log10(RMSE)0.000.250.500.751.00−6 −4 −2 0log10(AME)fitBtTh_unThBFigure 5.7: Empirical cumulative distribution functions of RMSE and AMEfor the top thee parameterizations (with respect to performance) where per-formance is greater than -8.41Chapter 6ConclusionsIn this study, we used the numerical procedure of fitting a Gaussian pro-cess for a large set of test problems. In fitting the Gaussian process, weused eleven different parameterization-optimization pairings in order to de-termine which parameterization performs the best so that it may be deemedthe go-to parameterization.Despite the likelihood function having the property of being invariant,numerically we have seen that this is not the case. Optimization producesdifferent results even though the optimizer started with equivalent values.We started by picking an appropriate design using the random latin hy-percube design using n = 7 · d, 10 · d, and 15 · d. We provided six differentparameterizations that have all been discussed in literature. We then usedthe statistical programming language R to optimize the log-likelihood func-tion in order to produce estimates for the parameter values for fitting themodel.The results of each of the eight hundred and forty tests performed pro-vided us with optimal log-likelihood values, and the parameter values usedto achieve the log-likelihood. We looked at three different tools in orderto evaluate which of the parameterizations performed the best: root meansquared error, absolute maximum error, and performance as defined in Sec-tion 5.4.We first visualized the distributions of the root mean squared errors andabsolute maximum errors. Four parameterizations performed better thanthe rest in terms of errors produced in predicting: beta, beta (no gradient),theta2, and theta (beta). We further investigated the performance metricalso using an empirical cumulative distribution function. From the ECDF’swe outlined three parameterizations that outperformed the other seven: beta(no gradient), theta2, and theta (beta). Between the errors and performancewe can see that the top three performing parametrizations are beta (nogradient), theta2 and theta (beta) as they achieve desirable log-likelihoodvalues and produce estimates that have low errors. We analyzed these threefurther in order to determine which is the best.The parameterizations had virtually identical error ECDF’s and there-42Chapter 6. Conclusionsfore the performance of the three was further analyzed. The three parame-terizations have similar performance curves up until -8. The area of differ-ence is in the right hand side of the graph where the solutions are far awayfrom the best solution found. We further analyzed this problematic area bylooking at the error values associated with the observations with larger per-formance values (greater than -8). The beta parameterization may performbetter than the other two performance wise, but the models achieved in theGaussian process produce poor estimates. The theta2 and theta (beta) pa-rameterizations have similar problematic area error curves. The tie-breakerbetween the theta2 and theta (beta) curves is time, specifically the timerequired to fit Gaussian process and obtain a model. The theta2 param-eterization requires only one run through the optimization process to fitand form the model whereas the theta (beta) parameterization requires onerun through the process first for the beta parameterization and then againthrough the theta2 parameterization.This thesis set out to find which parameterization of the correlationfunction in the Gaussian process performs the best. It was our wish todetermine which parameterization should be the recommended parameter-ization for those wanting to fit a Gaussian process. Based on the evidenceprovided in this thesis the parameterization denoted as being the best isthe parameterization in terms of θ2 using the analytical gradient and theunbounded optimizer.The performance of the parameterization rho is surprising. The rho pa-rameterization has had much success in Markov chain Monte Carlo settingsand further investigation into the parameter and why it performs poorlyshould be considered in the future.43Bibliography[1] Arnold, A. and Loeppky, J. L. (2015), “The Problem with AssessingStatistical Methods,” arXiv preprint arXiv:1510.01417, Submitted toThe American Statistician. → pages 32[2] Audet, C. and Dennis Jr, J. E. (2006), “Mesh adaptive direct search al-gorithms for constrained optimization,” SIAM Journal on optimization,17, 188–217. → pages 10[3] Butler, A., Haynes, R., Humphries, T., and Ranjan, P. (2014), “Effi-cient optimization of the likelihood function in Gaussian process mod-elling,” Computational Statistics & Data Analysis, 73, 40–52. → pages1, 16[4] Chen, H., Loeppky, J. L., Sacks, J., and Welch, W. J. (2016), “Anal-ysis Methods for Computer Experiments: How to Assess and WhatCounts?” Statistical Science, 31, 40–60. → pages 5[5] Currin, C., Mitchell, T., Morris, M., and Ylvisaker, D. (1991),“Bayesian Prediction of Deterministic Functions, With Applicationsto the Design and Analysis of Computer Experiments,” Journal of theAmerican Statistical Association, 86, 953–963. → pages 1, 7[6] Dennis, J. and Woods, D. J. (1987), “Optimization on microcomputers:The Nelder-Mead simplex Algorithm,” New computing environments:microcomputers in large-Scale computing, 116–122. → pages 10[7] Handcock, M. S. and Stein, M. L. (1993), “A Bayesian Analysis ofKriging,” Technometrics, 35, 403–410. → pages 6[8] Higdon, D., Gattiker, J., Williams, B., and Rightley, M. (2008), “Com-puter Model Calibration Using High-Dimensional Output,” Journal ofthe American Statistical Association, 103, 570–583. → pages 16, 38[9] Johnson, M. E., Moore, L. M., and Ylvisaker, D. (1990), “Minimaxand Maximin Distance Designs,” Journal of Statistical Planning andInference, 26, 131–148. → pages 1544Bibliography[10] Linkletter, C., Bingham, D., Hengartner, N., Higdon, D., and Ye, K. Q.(2006), “Variable Selection for Gaussian Process Models in ComputerExperiments,” Technometrics, 48, 478–490. → pages 16, 38[11] Loeppky, J. L., Sacks, J., and Welch, W. J. (2009), “Choosing theSample Size of a Computer Experiment: A Practical Guide,” Techno-metrics, 51, 366–376. → pages 29[12] McKay, M. D., Beckman, R. J., and Conover, W. J. (1979), “A Com-parison of Three Methods for Selecting Values of Input Variables inthe Analysis of Output from a Computer Code,” Technometrics, 21,239–245. → pages 13[13] Morris, M. D. and Mitchell, T. J. (1995), “Exploratory Designs forComputational Experiments,” Journal of Statistical Planning and In-ference, 43, 381–402. → pages 13[14] Mu¨hlenbein, H., Schomisch, M., and Born, J. (1991), “The parallelgenetic algorithm as function optimizer,” Parallel computing, 17, 619–632. → pages 10[15] O’Hagan, A. (1978), “Curve Fitting and Optimal Design for Predic-tion,” JRSSB, 40, 1–42. → pages 4[16] Ranjan, P., Haynes, R., and Karsten, R. (2011), “A ComputationallyStable Approach to Gaussian Process Interpolation of DeterministicComputer Simulation Data,” Technometrics, 53, 366–378. → pages 10,16[17] Rasmussen, C. E. and Williams, C. K. (2006), Gaussian processes formachine learning, vol. 1, MIT press Cambridge. → pages 4[18] Rasmussen, C. E. and Williams, C. K. I. (2006), Gaussian Processesfor Machine Learning, Cambridge, MA: The MIT Press. → pages 8[19] Roustant, O., Ginsbourger, D., and Deville, Y. (2012), “DiceKriging,DiceOptim: Two R packages for the analysis of computer experimen-tRouGinDev2012s by kriging-based metamodeling and optimization,”. → pages 8, 16[20] Sacks, J., Schiller, S. B., and Welch, W. J. (1989), “Designs for Com-puter Experiments,” Technometrics, 31, 41–47. → pages 1, 1345Bibliography[21] Sacks, J., Welch, W. J., Mitchell, T. J., and Wynn, H. P. (1989), “De-sign and Analysis of Computer Experiments (with Discussion),” Sta-tistical Science, 4, 409–435. → pages 8[22] Santner, T. J. ., Williams, B. J. ., and Notz, W. I. (2003), The Designand Analysis of Computer Experiments, New York: Springer. → pages1[23] Santner, T. J., Williams, B. J., and Notz, W. I. (2013), The design andanalysis of computer experiments, Springer Science & Business Media.→ pages 13[24] Shanno, D. F. and Kettler, P. C. (1970), “Optimal Conditioning ofQuasi-Newton Methods,” Mathematics of Computation, 24, 657–664.→ pages 10[25] Stein, M. L. (1999), Interpolation of Spatial Data, Some Theory forKriging, New York: Springer. → pages 8[26] Surjanovic, S. and Bingham, D. (2015), “Virtual Library of SimulationExperiments:,” . → pages 47[27] Wang, J. M., Fleet, D. J., and Hertzmann, A. (2008), “Gaussian processdynamical models for human motion,” IEEE transactions on patternanalysis and machine intelligence, 30, 283–298. → pages 8[28] Welch, W. J., Buck, R. J., Sacks, J., Wynn, H. P., Mitchell, T. J., andMorris, M. D. (1992), “Screening, Predicting, and Computer Experi-ments,” Technometrics, 34, 15–25. → pages 7, 13[29] Zehna, P. W. et al. (1966), “Invariance of maximum likelihood estima-tors,” Annals of Mathematical Statistics, 37, 744. → pages 18[30] Zhu, C., Byrd, R. H., Lu, P., and Nocedal, J. (1997), “Algorithm 778:L-BFGS-B: Fortran subroutines for large-scale bound-constrained opti-mization,” ACM Transactions on Mathematical Software (TOMS), 23,550–560. → pages 1046AppendixAppendix ATest FunctionsThe fourteen functions used for the simulation study are from [26]. Thefunctions are as follows:Borehole function (d = 8)f(x) =2pix3(x4 − x6)ln(x2/x1)(1 + 2x7x3ln(x2/x1)x21x8)+ x3x5x1 ∈ [0.05, 0.15]x2 ∈ [100, 50 000]x3 ∈ [63 070, 115 600]x4 ∈ [990, 1 110]x5 ∈ [63.1, 116]x6 ∈ [700, 820]x7 ∈ [1 120, 1 680]x8 ∈ [9 855, 12 045]Branin function (d = 2)f(x) =(x2 − 5.14pi2x21 +5pix1 − 6)2+ 10(1− 18pi)cos(x1) + 10x1 ∈ [−5, 10]x2 ∈ [0, 15]47Currin et al. (1988) Exponential function (d = 2)f(x) =(1− exp{− 12x2})2300x31 + 1900x21 + 2092x1 + 60100x31 + 500x21 + 4x1 + 20xi ∈ [0, 1] ∀i = 1, 2Dette and Pepelyshev (2010) 8-Dimensional function (d = 8)f(x) = 4(x1−2+8x2−8x22)2+(3−4x2)2+16√x3 + 1(2x3−1)2+8∑i=4i ln1 + i∑j=3jxi ∈ [0, 1] ∀i = 1, . . . , 8Dette and Pepelyshev (2010) Curved function (d = 3)f(x) = 4(x1 − 2 + 8x2 − 8x22)2 + (3− 4x2)2 + 16√x3 + 1(2x3 − 1)2xi ∈ [0, 1] ∀i = 1, 2, 3Dette and Pepelyshev (2010) Exponential function (d = 3)f(x) = 100(exp{ −1x1.751}+ exp{ −2x1.52}+ exp{ −2x1.253})xi ∈ [0, 1] ∀i = 1, 2, 3Friedman function (d = 5)f(x) = 10 sin(pix1x2) + 20(x3 − 0.5)2 + 10x4 + 5x5xi ∈ [0, 1] ∀i = 1, . . . , 5Gramacy and Lee (2009) (d = 6)f(x) = exp{sin((0.9(x1 + 0.48))10)}+ x2x3 + x4x5 and x6 are inactive.A small error term ε(0, 0.052) is added.xi ∈ [0, 1] ∀i = 1, . . . , 648Lim et al. (2002) Nonpolynomial function (d = 2)f(x) =16[(30 + 5x1 sin(5x1))(4 + exp{−5x2})− 100]xi ∈ [0, 1] ∀i = 1, 2Linkletter et al. (2006) Decreasing Coefficients function (d = 10)f(x) = 0.2x1 +0.22x2 +0.24x3 +0.28x4 +0.216x5 +0.232x6 +0.264x7 +0.2128x8x9 and x10 are inactive. A small error term ε(0, 0.052) is added.xi ∈ [0, 1] ∀i = 1, . . . , 10Linkletter et al. (2006) Sinusoidal function (d = 10)f(x) = sin(x1) + sin(5x2)x3, . . . , x10 are inactive. A small error term ε(0, 0.052) is added.xi ∈ [0, 1] ∀i = 1, . . . , 10OTL Circuit function (d = 6)f(x) =(12x2x1+x2+ 0.74)x6(x5 + 9)x6(x5 + 9) + x3+11.35x3x6(x5 + 9) + x3+0.74x3x6(x5 + 9)(x6(x5 + 9) + x3)x4x1 ∈ [50, 150]x2 ∈ [25, 70]x3 ∈ [0.5, 3]x4 ∈ [1.2, 2.5]x5 ∈ [0.25, 1.2]x6 ∈ [50, 300]Piston Simulation function (d = 7)f(x) = 2pi√x1x4 + x22x5x3x7x6V 2,49whereV =x22x4(√A2 + 4x4x5x3x7x6 −A)andA = x5x2 + 19.62x1 − x4x3x2x1 ∈ [30, 60]x2 ∈ [0.005, 0.020]x3 ∈ [0.002, 0.010]x4 ∈ [1 000, 5 000]x5 ∈ [90 000, 110 000]x6 ∈ [290, 296]x7 ∈ [340, 360]Wing Weight function (d = 10)f(x) = 0.036x0.7581 x0.00352(x3cos2(x4))0.6x0.0065 x0.046(100x7cos(x4))−0.3(x8x9)0.49+x1x10x1 ∈ [150, 200]x2 ∈ [220, 300]x3 ∈ [6, 10]x4 ∈ [−10, 10]x5 ∈ [16, 45]x6 ∈ [0.5, 1]x7 ∈ [0.08, 0.18]x8 ∈ [2.5, 6]x9 ∈ [1 700, 2 500]x10 ∈ [0.025, 0.08]50
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
Attribution-NoDerivatives 4.0 International
Graduate
A numerical study of the effects of the parameterization of the gaussian process
Text