Bayesian treed Gaussian process models with an application to computer modeling

Much of modern engineering design is now done through computer modeling, which is both faster and more cost-effective than building small-scale models, particularly in the earlier stages of design when more scope for changes is desired. As design proceeds, increasingly sophisticated simulators may be created. Our work here was motivated by a simulator of a proposed rocket booster. NASA relies heavily on simulators for design, as wind tunnel experiments are quite expensive and still not fully realistic of the range of flight experiences. In particular, one of the highly critical points in time for a rocket booster is the moment that it re-enters the atmosphere. Such conditions are difficult to recreate in a wind tunnel, and it is obviously impossible to run a standard physical experiment. Thus to learn about the behavior of the proposed rocket booster, NASA uses computer simulation.

Simulators can involve large amounts of physical modeling, requiring the numerical solution of complex systems of differential equations. The NASA simulator was no exception, typically requiring between five and twenty hours for a single run. Thus, NASA was interested in creating a statistical model of the simulator itself, an emulator or surrogate model, in the terminology of computer modeling. The standard approach in the literature for emulation is to model the simulator output with a stationary smooth Gaussian process (GP) (Sacks et al., 1989;Kennedy and O’Hagan, 2001;Santner et al., 2003). However, this approach proved to be inadequate for the NASA data. In particular, we were faced with problems of nonstationarity, heteroscedasticity, and the size of the dataset. Thus we introduce here an expansion of GPs, based on the idea of Bayesian partition models (Chipman et al., 2002;Denison et al., 2002), which is able to address these issues.

GPs are conceptually straightforward, can easily accommodate prior knowledge in the form of covariance functions, and can return estimates of predictive confidence, which were desired by NASA. However, we highlight three disadvantages of the standard form of a GP which we had to confront on this dataset, and expect to encounter on a wide range of other applications. First, inference on the GP scales poorly with the number of data points, N, typically requiring computing time in O(N 3 ) for calculating inverses of N × N covariance matrices. Second, GP models are usually stationary in that the same covariance structure is used throughout the entire input space, which may be too strong of a modeling assumption. Third, the estimated predictive error (as opposed to the predictive mean value) of a stationary model does not directly depend on the locally observed response values. Rather, the predictive error at a point depends only on the locations of the nearby observations and on a global measure of error that uses all of the discrepancies between observations and predictions without regard for their distance from the point of interest (because of the stationarity assumption). (Section 4.3 provides more details, in particular note that Eq. ( 12) does not depend on z.) In many real-world spatial and stochastic problems, such a uniform modeling of uncertainty will not be desirable. Instead, some regions of the space will tend to exhibit larger variability than others. On the other hand, fully nonstationary Bayesian GP models (e.g., Higdon et al., 1999;Schmidt and O’Hagan, 2003) can be difficult to fit, and are not computationally tractable for more than a relatively small number of datapoints. Further discussion of nonstationary models is deferred until the end of Section 3.2.

All of these shortcomings can be addressed by partitioning the input space into regions, and fitting separate stationary GP models within each region (e.g., Kim et al., 2005). Partitioning provides a relatively straightforward mechanism for creating a nonstationary model, and can ameliorate some of the computational demands by fitting models to less data. A Bayesian model averaging approach allows for the explicit estimation of predictive uncertainty, which can now vary beyond the constraints of a stationary model. Finally, an R package with implementations of all of the models dis-cussed in this paper is available at http://www.cran.r-project.org/web/packages/tgp/index.html . We note that by partitioning, we do not have any theoretical guarantee of continuity in the fitted function. However, as we demonstrate in several examples, Bayesian model averaging yields mean fitted functions that are typically quite smooth in practice, giving fits that are indistinguishable from continuous functions except when the data call for the contrary.

Indeed the ability to accurately model possible discontinuities is a side benefit of this approach.

The rest of the paper is organized as follows. Section 2 describes the motivating data in further detail. Section 3 provides some background material. Section 4 combines stationary GPs and treed

…(Full text truncated)…

This content is AI-processed based on ArXiv data.