GPdoemd is an open-source python package for design of experiments for model discrimination that uses Gaussian process surrogate models to approximate and maximise the divergence between marginal predictive distributions of rival mechanistic models. GPdoemd uses the divergence prediction to suggest a maximally informative next experiment

Misener, R

Olofsson, S

English

Spiral - Imperial College Digital Repository

arXiv:1810.02561v1  [cs.MS]  5 Oct 2018GPdoemd: a python package for design of experiments formodel discrimination using Gaussian process surrogatesSimon Olofsson simon.olofsson15@imperial.ac.ukRuth Misener r.misener@imperial.ac.ukDepartment of Computing, Imperial College, London, SW7 2AZ, United Kingdom–AbstractGPdoemd is an open-source python package for design of experiments for model dis-crimination that uses Gaussian process surrogate models to approximate and maximisethe divergence between marginal predictive distributions of rival mechanistic models. GP-doemd uses the divergence prediction to suggest a maximally informative next experiment.Keywords: Software, Gaussian processes, design of experiments, model discriminationExplicit parametric mechanistic models are common in science and engineering, e.g. eco-nomics (Black and Scholes, 1973), biology (Mehrian et al., 2018), and control theory (Bem-porad et al., 2002). Researchers develop multiple rival mechanistic models, correspondingto different hypotheses about underlying system mechanisms, but typically lack sufficientdata to discriminate between the models (Box and Hill, 1967; Buzzi-Ferraris and Forzatti,1983). To discard inaccurate models, GPdoemd suggests gathering more data through ad-ditional experiments. To minimise cost, e.g. money and time, GPdoemd designs optimalnew experiments, i.e. experiments yielding maximally informative results.1. BackgroundGPdoemd assumes multiple rival mechanistic models fi, e.g. systems of ODEs, which take asinputs some design variables x, e.g. temperature, pressure, flow rate, and model parametersθi, e.g. chemical reaction rate. Typically, the model parameters θi are tuned to make themodel predictions fit the experimental data D. However, there will be model parameteruncertainty, i.e. θi ∼ p(θi | D), which propagates through to uncertainty in the modelpredictions. To account for the model parameter uncertainty, we approximate the models’marginal predictive distributionsp(fi(x) | D) =∫p(fi(x) |θi)p(θi | D)dθi . (1)The optimal next experiment x∗ is found by maximising a design utility functionx∗ = arg maxxU {fi(x) | D ; i = 1, . . . ,M} , (2)where U{·} is a divergence measure between the rival models’ marginal predictive distribu-tions. The idea is to find the experiment x∗ for which the model predictions differ the most,within the error margin imposed by the model parameter uncertainty and measurementnoise variance. Equation (2) is typically intractable, so the solution must be approximated.c©2018 Simon Olofsson and Ruth Misener.Olofsson and MisenerTypical methods for design of experiments for model discrimination can be divided intotwo different approaches (Olofsson et al., 2018). First, the analytical approach, e.g. Box andHill (1967), Buzzi-Ferraris et al. (1990) and Michalik et al. (2010), assumes that the mecha-nistic models are (approximately) linear in the model parameters, and the model parameterposterior p(θi | D) is a multivariate Gaussian distribution. The analytical approach is com-putationally efficient but requires gradient information that may not be readily available incomplex models or legacy code. The second approach is the data-driven, i.e. Monte Carlo-based, approach, e.g. Vanlier et al. (2014) and Ryan et al. (2015), where samples are drawnfrom the model parameter posterior in order to compute the optimal next experiment. Nogradient information is required, but the data-driven approach can be computationally veryexpensive due to the number of samples required for large x and θi search spaces.Olofsson et al. (2018) hybridise the analytical and data-driven approaches by trainingGaussian process (GP) surrogate models on sampled data and using the surrogate models inan analytical fashion. The GPdoemd open-source software package implements this hybridapproach. The hybrid approach computes p(fi(x) | D) = N (µ̆, Σ̆) by placing a GP priorfi ∼ GP(0, kxkθ) on the original mechanistic models fi, assuming a Gaussian model param-eter distribution θ | D ∼ N (θ∗,Σθ), and then using closed-form approximations to computethe first two moments µ̆ and Σ̆ of (1). Design utility functions U{·} from classical literature,which are closed-form for Gaussian distributions, find the optimal next experiment.2. ImplementationGPdoemd, available from the GitHub repository https://github.com/cog-imperial/GPdoemd,is a python package implementing the Olofsson et al. (2018) hybrid approach. It uses func-tionality from the GPy (since 2012) python package for GP inference. The only otherdependencies are on the standard numpy (v1.7+) and scipy (v0.17+) packages. GPdoemdis regularly tested for python v3.4 and up. On GitHub, we provide documentation forinstalling and using GPdoemd via a PDF and Jupyter notebook demonstrations.The modular toolbox, illustrated in Figure 1, offers a choice between different GP kernelfunctions, inference methods, methods to approximate the marginal predictive distributions,design utility functions and model discrimination criteria. New modules can be easilyimplemented and added to the GPdoemd toolbox. The toolbox currently comes with theTable 1 case studies which are mostly from literature, although collaborators developedmixing. Researchers and engineers may try the Olofsson et al. (2018) hybrid approach andcompare the performance to competing methods.3. Syntax and Supported FeaturesAssuming the rival models have been proposed, GPdoemd assists in model discrimination.Model type A model object is initialised using a python dictionary containing the modelname (name), the model function f(x,θi) handle (call), the design variable and model pa-rameter dimensions (dim x and dim p), the number of target dimensions (num outputs), modelparameter bounds (p bounds), experimental measurement noise (co)variance (meas noise var),and a list of discrete design variable dimensions (binary variables). This dictionary is passedto one of the implemented model types (Box 1 in Figure 1).2GPdoemd: a python package for design of experiments for model discriminationModel dictionary 1. ModelAnalyticNumericalGPModelSparseGPModelGPGriefModel3. GP kernelsRBFExponentialMatern32Matern52CosineRatQuad2. Param. estimationdiff evolleast squares4. Param. covariancelaplace approximation5. Approximativemarginalisationfirst order taylorsecond order taylor6. Design criterionHRBHBFAW7. Discriminationcriteriongaussian posteriorchi2akaikeFigure 1: The modular structure of the GPdoemd open-source software package.Table 1: GPdoemd case studies, with the number of design variables |x| (continuous ordiscrete), model parameters |θi|, target dimensions |y| and rival models M .Name Reference |x|, (∈ Z) |θi| |y| M Typebff1983 Buzzi-Ferraris and Forzatti (1983) 3, (0) 5 1 5 Analyticbffeh1984 Buzzi-Ferraris et al. (1984) 2, (0) 4 2 4 Analyticbffc1990a Buzzi-Ferraris et al. (1990) 3, (0) 2–6 1 4 Analyticmixing - 3, (1) 1 1 5 Analyticmsm2010 Michalik et al. (2010) 3, (0) 1 1 10 Analyticvthr2014linear Vanlier et al. (2014) 1, (0) 2–4 1 4 Analyticvthr2014ode Vanlier et al. (2014) 3, (2) 14 1 4 Black-boxtandogan2017 Tandogan et al. (2017) 4, (0) 8–14 2 3 Black-boxParameter estimation Given experimental data Ydata for designs Xdata, GPdoemd opti-mises the model parameter values θ∗ using a prediction error minimisation methods (Box 2in Figure 1): differential evolution (diff evol) or least squares with finite difference gradientapproximation (least squares). Both optimisation methods call functions in scipy.GP kernels The GP surrogate models require a choice of GP kernel functions kx and kθfor the GP prior GP(0, kxkθ). GPdoemd currently supports six different kernel functions(Box 3 in Figure 1) taken from the GPy package and extended with function calls for thesecond derivatives ∂2k(r)/∂r2 with respect to the distance measure r.Model parameter covariance (Σθ) GPdoemd uses a Gaussian approximation N (θ∗,Σθ)of the model parameter distribution and implements a Laplacian approximation of Σθ.Approximating marginal predictive distributions The Olofsson et al. (2018) hybridapproach approximates the marginal predictive distribution in (1) with a Gaussian distribu-tion. The GPdoemd implements two different methods of computing the first two momentsof (1): first- and second-order Taylor approximations (Box 5 of Figure 1), using the firstand second derivatives of the GP surrogate’s predictive mean and variance.3Olofsson and Misenerxp(fi(x)|D)f1 f2 f3xU{fi(x)|D}BH BF AWFigure 2: (Left) Model predictive distributions over designs x for rival models f1, f2, f3.(Right) Three design utility functions for the model predictive distributions on theleft. Note that the different utilities, i.e. Eq. (2), prefer different next experiments.Design criterion GPdoemd provides five different utility functions (Box 6 in Figure 1)for designing the next experiment: HR (Hunter and Reiner, 1965), BH (Box and Hill, 1967),BF (Buzzi-Ferraris et al., 1990) and AW (Michalik et al., 2010). Figure 2 illustrates thedifference between a few of these design criteria.Discrimination criterion GPdoemd provides three different criteria (Box 7 in Figure 1)for model discrimination: normalised Gaussian posteriors (Box and Hill, 1967), χ2 test(Buzzi-Ferraris and Forzatti, 1983), and the Akaike information criterion weights (Michaliket al., 2010). Olofsson et al. (2018) discuss trade-offs between the discrimination criteria.3.1 ExampleAssume a list dlist of model dictionaries, experimental data Xdata, Ydata with measurementnoise variance measvar, and lists X, P and Y of surrogate model training data (design, modelparameters and predictions, respectively) are given. We wish to select the optimal nextexperiment from candidates Xnew.N = Xnew.shape[0] # Number of test pointsM = len( dlist ) # Number of rival modelsE = Ydata.shape[1] # Number of target dimensionsmu, s2 = np.zeros(( N, M, E )), np.zeros(( N, M, E, E ))for i,d in enumerate( dlist ):m = GPdoemd.models.GPModel(d) # Initialise surrogate model# Estimate model parameter valuesm.param estim(Xdata , Ydata , GPdoemd.param estim.least squares , m.p bounds )# Set−up surrogate modelRBF = GPdoemd.kernels.RBFm.gp surrogate (Z=np.c [X[i], P[i]], Y=Y[i], kern x=RBF, kern p=RBF)m.gp optimise ()# Approximate model parameter covariancem.Sigma = GPdoemd.param covar.laplace approximation( m, Xdata )# Approximate marginal predictive distribution at test pointsmu[:,i], s2[:,i] = GPdoemd.marginal.taylor first order( m, Xnew )dc = GPdoemd.design criteria.AW(mu, s2, measvar) # Design criterion at test pointsxnext = Xnew[ np.argmax(dc) ] # Optimal next experiment4GPdoemd: a python package for design of experiments for model discriminationAcknowledgmentsThis work has received funding from the European Union’s Horizon 2020 research andinnovation programme under the Marie Sk lodowska-Curie grant agreement no.675251, andan EPSRC Research Fellowship to R.M. (EP/P016871/1).ReferencesA. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos. The explicit linear quadraticregulator for constrained systems. Automatica, 38:3–20, 2002.F. Black and M. Scholes. The pricing of options and corporate liabilities. J Political Econ,81(3):637–354, 1973.G. E. P. Box and W. J. Hill. Discrimination among mechanistic models. Technometrics, 9(1):57–71, 1967.G. Buzzi-Ferraris and P. Forzatti. A new sequential experimental design procedure fordiscriminating among rival models. Chem Eng Sci, 38(2):225–232, 1983.G. Buzzi-Ferraris, P. Forzatti, G. Emig, and H. Hofmann. Sequential experimental designfor model discrimination in the case of multiple responses. Chem Eng Sci, 39(1):81–85,1984.G. Buzzi-Ferraris, P. Forzatti, and P. Canu. An improved version of a sequential designcriterion for discriminating among rival multiresponse models. Chem Eng Sci, 45(2):477–481, 1990.GPy. GPy: A Gaussian process framework in python.http://github.com/SheffieldML/GPy, since 2012.W. G. Hunter and A. M. Reiner. Designs for discriminating between two rival models.Technometrics, 7(3):307–323, 1965.M. Mehrian, Y. Guyot, I. Papantoniou, S. Olofsson, M. Sonnaert, R. Misener, and L. Geris.Maximizing neotissue growth kinetics in a perfusion bioreactor: An in silico strategyusing model reduction and Bayesian optimization. Biotechnol Bioeng, 115(3):617–629,2018.C. Michalik, M. Stuckert, and W. Marquardt. Optimal experimental design for discriminat-ing numerous model candidates: The AWDC criterion. Ind Eng Chem Res, 49:913–919,2010.S. Olofsson, M. P. Deisenroth, and R. Misener. Design of experiments for model discrimi-nation hybridising analytical and data-driven approaches. In ICML ’18: Proceedings ofthe International Conference on Machine Learning, Stockholm, Sweden, PMLR 80, 2018.E. G. Ryan, C. C. Drovandi, and A. N. Pettitt. Fully Bayesian experimental design forpharmacokinetic studies. Entropy, 17:1063–1089, 2015.N. Tandogan, S. Garcia-Muñoz, M. Sen, T. M. Wilson, J. Y. Buser, S. P. Kolis, I. V.Borkar, and C. A. Alt. Use of model discrimination method in drug substance processdevelopments. In Proceedings of the AIChE Annual Meeting, pages X–Y, 2017.J. Vanlier, C. A. Tiemann, P. A. J. Hilbers, and N. A. W. van Riel. Optimal experimentdesign for model selection in biochemical networks. BMC Syst Biol, 8(20), 2014.5

GPdoemd: a python package for design of experiments for model discrimination

http://spiral.imperial.ac.uk/bitstream/10044/1/63850/2/1810.02561v1.pdf

GPdoemd: a python package for design of experiments for model discrimination

Abstract

Similar works

Full text

Available Versions

Spiral - Imperial College Digital Repository