Gaussian Process Regression

A simulation is seen as a function \(f(x)+\epsilon\) (\(x \in \mathbb{R}^D\)) with additional random error \(\epsilon \sim \mathcal{N}(0,v_t)\). This optional error is due to the stochastic nature of most simulations.

Note

This error does not depend on \(x\). An extension to model an error depending on \(x\) will be added in the future.

The GaussianProcess module uses regression to model the simulation as a Gaussian process. We refer to Girard’s thesis [1] for a really good explanation of GP regression.

“A Gaussian process is a collection of random variables, any finite number of which have (consistent) joint Gaussian distributions.” – Rasmussen, C. E. Gaussian Processes in Machine Learning, Advanced Lectures on Machine Learning, Springer Berlin Heidelberg, 2004, 3176, 63-71

“Once the mean and covariance functions are defined, everything else about GPs follows from the basic rules of probability applied to multivariate Gaussians.” – Zoubin Ghahramani

Note

The following explanation is an excerpt from my paper [2].

\(m(x) = E[G(x)]\) (\(x\in\mathbb{R}^d\)) is the mean function of the Gaussian process and \(C(x_i,x_j) = \mathrm{Cov}[G(x_i),G(x_j)]\) is the covariance function. For any given set of inputs \(\{x_1,\dots,x_n\}\), \((G(x_1),\dots,G(x_n))\) is a random vector with an \(n\)-dimensional normal distribution \(\mathcal{N}({\mu},\Sigma)\). \({\mu}\) is the vector of mean values \((m(x_1),\dots,m(x_n))\) and \(\Sigma\) is the covariance matrix with \(\Sigma_{ij} = C(x_i,x_j)\)

For modeling \(f(x)\), we use the Gaussian squared exponential covariance function \(C(x_i,x_j) = v \exp \left[-\tfrac{1}{2}(x_i-x_j)^T {W}^{-1}(x_i-x_j)\right]\), as it performs well regarding accuracy. Without loss of generality, we use \(m(x)=m\).

\({W}^{-1} = \mathrm{diag}(w_1,\dots,w_d)\) and \(v\) are hyperparameters of the Gaussian process. We have a set of observed simulation results \(\mathcal{B} = \{(x_i,t_i)\}_{i=1}^N\) as supporting points with \(x_i \in \mathbb{R}^d\) and \(t_i = f(x_i)+\epsilon_i\) with white noise \(\epsilon_i \sim \mathcal{N}(0,v_t)\). This noise represents the aleatory uncertainty of the simulation that calculates \(f(x)\).

To model the noisy observations \(t_i\) we use \(K=\Sigma + v_t I\) as covariance matrix for regression. Therefore, we have an additional hyperparameter \(v_t\). The hyperparameters \(W^{-1}\), \(v_t\), and \(v\) can be estimated from a given dataset with a maximum likelihood approach. In the following, we use \({\beta}={K}^{-1}{t}\) with \({t} = (t_1,\dots,t_N)^T\) to simplify notation.

Now, the Gaussian process can be used as a surrogate for the simulation model. The prediction \(\mu(x)\) and code plus aleatory uncertainty \(\sigma^2 (x)\) at a new input \(x\) is:

\[\begin{split}\mu (x) &= \sum_{i=1}^N \beta_i C(x,x_i) \\ \sigma^2 (x) &= C(x,x) - \sum_{i,j=1}^N K_{ij}^{-1}C(x,x_i)C(x,x_j) + v_t\end{split}\]

1-d Example:

Usage: See Getting Started

References

[1]	Girard, A. Approximate Methods for Propagation of Uncertainty with Gaussian Process Models, University of Glasgow, 2004

[2]	Baumgaertel, P.; Endler, G.; Wahl, A. M. & Lenz, R. Inverse Uncertainty Propagation for Demand Driven Data Acquisition, Proceedings of the 2014 Winter Simulation Conference, IEEE Press, 2014, 710-721 (https://www6.cs.fau.de/publications/public/2014/WinterSim2014_baumgaertel.pdf)