Gaussian processes are a powerful tool for nonlinear regression models.
Assume that we have predictor variables $\mathbf{X} = \{\mathbf{x_i}\}_{i=1}^N \in \mathbb{R}^d$ and response variables $\mathbf{y}=\{y_i \in \mathbb{R}\}_{i=1}^N$.
The response variables $\mathbf{y}$ are assumed to dependent on the predictors $\mathbf{X}$,
\begin{equation} y_i \sim \mathcal{N}(f(\mathbf{x}_i),\sigma^2), \ i=1,\ldots,n, \end{equation}where $f$ is a mapping function. Treating $f$ as a random function, we assume that the distribution over $f$ is a Gaussian process, $$ f \sim \mathcal{GP}(m(\mathbf{x}),k(\mathbf{x},\mathbf{x}')), $$ where $m(\cdot)$ and $k(\cdot,\cdot)$ are the mean and kernel functions respectively.
We start by simulating some data
using GaussianProcesses
srand(20140430)
# Training data
n=10; #number of training points
x = 2π * rand(n); #predictors
y = sin.(x) + 0.05*randn(n); #regressors
The first step in modelling with Gaussian Processes is to choose mean functions and kernels which describe the process.
Note that all hyperparameters for the mean and kernel functions and $\sigma$ are given on the log scale. This is true for all strictly positive hyperparameters. Gaussian Processes are represented by objects of type 'GP' and constructed from observation data, a mean function and kernel, and optionally the amount of observation noise.
#Select mean and covariance function
mZero = MeanZero() #Zero mean function
kern = SE(0.0,0.0) #Sqaured exponential kernel (note that hyperparameters are on the log scale)
logObsNoise = -1.0 # log standard deviation of observation noise (this is optional)
gp = GP(x,y,mZero,kern,logObsNoise) #Fit the GP
GP Exact object: Dim = 1 Number of observations = 10 Mean function: Type: GaussianProcesses.MeanZero, Params: Float64[] Kernel: Type: GaussianProcesses.SEIso, Params: [0.0, 0.0] Input observations = [4.85461 5.17653 … 1.99412 3.45676] Output observations = [-0.967293, -1.00705, -1.0904, 0.881121, -0.333213, -0.976965, 0.915934, 0.736218, 0.950849, -0.306432] Variance of observation noise = 0.1353352832366127 Marginal Log-Likelihood = -6.335
Once we've fit the GP
function to the data, we can calculate the predicted mean and variance of of the function at unobserved points $\{\mathbf{x}^\ast,y^\ast\}$, conditional on the observed data $\mathcal{D}=\{\mathbf{y},\mathbf{X}\}$. This is done with the predict_y
function.
The predict_y
function returns the mean vector $\mu(\mathbf{x}^\ast)$ and covariance matrix (variance vector if full_cov=false
) $\Sigma(\mathbf{x}^\ast,\mathbf{x}^{\ast^\top})$ of the predictive distribution,
\begin{equation}
y^\ast|\mathbf{x}^\ast,\mathcal{D} \sim \mathcal{N}(\mu(\mathbf{x}^\ast),\Sigma(\mathbf{x}^\ast,\mathbf{x}^{\ast^\top})+\sigma^2\mathbf{I}),
\end{equation}
where
\begin{eqnarray}
\mu(\mathbf{x}^\ast) &=& k(\mathbf{x}^\ast,\mathbf{X})(k(\mathbf{X}, \mathbf{X}) + \sigma_n^2 \mathbf{I})^{-1}\mathbf{y} \\ \ \mbox{and} \
\Sigma(\mathbf{x}^\ast,\mathbf{x}^{\ast}) &=& k(\mathbf{x}^\ast,\mathbf{x}^\ast) -k(\mathbf{x}^\ast,\mathbf{X})(k(\mathbf{X}, \mathbf{X})+ \sigma_n^2 \mathbf{I})^{-1} k(\mathbf{X},\mathbf{x}^\ast).
\end{eqnarray}
Note you can use the predict_f
function to predict the latent function $\mathbf{f}^\ast$.
μ, σ² = predict_y(gp,linspace(0,2π,100));
Plotting GPs is straightforward and utilises the recipes approach to plotting from the Plots.jl package. More information about plotting GPs and the available functionality can be found in this notebook.
The default plot function plot(gp)
outputs the predicted mean and variance of the function (i.e. uses predict_f
in the background), with the uncertainty in the function represented by a confidence ribbon (set to 95% by default). All optional plotting arguments are given after ;
.
using Plots #Load Plots.jl package
plot(gp; xlabel="x", ylabel="y", title="Gaussian process", legend=false, fmt=:png) # Plot the GP
The hyperparameters are optimized using the Optim.jl package. This offers users a range of optimization algorithms which can be applied to estimate the hyperparameters using type II maximum likelihood estimation. Gradients are available for all mean and kernel functions used in the package and therefore it is recommended that the user utilizes gradient based optimization techniques. As a default, the optimize!
function uses the L-BFGS
solver, however, alternative solvers can be applied.
optimize!(gp; method=Optim.BFGS()) # Optimise the hyperparameters
Dict(:mean=>true,:kern=>true,:noise=>true)
Results of Optimization Algorithm * Algorithm: BFGS * Starting Point: [-1.0,0.0,0.0] * Minimizer: [-2.99484893977448,0.4636861184637202, ...] * Minimum: -3.275745e+00 * Iterations: 9 * Convergence: true * |x - x'| < 1.0e-32: false * |f(x) - f(x')| / |f(x)| < 1.0e-32: true * |g(x)| < 1.0e-08: true * f(x) > f(x'): false * Reached Maximum Number of Iterations: false * Objective Function Calls: 42 * Gradient Calls: 42
plot(gp) #Plot the GP after the hyperparameters have been optimised
MCMC can be run on the GPE hyperparameters through the mcmc
function. Priors for hyperparameters of the mean and kernel parameters can be set through the set_priors!
function. The log noise parameter of the GPE is a Uniform(0,1) distribution and currently can't be changed.
using Distributions
set_priors!(kern, [Normal(), Normal()]) # Uniform(0,1) distribution assumed by default if priors not specified
chain = mcmc(gp)
plot(chain', label=["Noise", "SE log length", "SE log scale"])
BasicMCJob: Variable [1]: p (BasicContMuvParameter) GenericModel: 1 variables, 0 dependencies (directed graph) HMC sampler: number of leaps = 10, leap step = 0.1 VanillaMCTuner: period = 100, verbose = false BasicMCRange: number of steps = 1000, burnin = 0, thinning = 1
The regression example above can be easily extended to higher dimensions. For the purpose of visualisation, and without loss of generality, we consider a 2 dimensional regression example.
#Training data
d, n = 2, 50; #Dimension and number of observations
x = 2π * rand(d, n); #Predictors
y = vec(sin.(x[1,:]).*sin.(x[2,:])) + 0.05*rand(n); #Responses
For problems of dimension>1 we can use isotropic (Iso
) kernels or automatic relevance determination (ARD
) kernels. For Iso
kernels, the length scale parameter $\ell$ is the same for all dimensions. For ARD
kernels, each dimension has different length scale parameter.
The Iso
and ARD
kernels are implemented automatically by replacing the single length scale parameter with a vector of parameters. For example, below we use the Matern 5/2 ARD
kernel, if we wanted to use the Iso
alternative then we would set the kernel as kern=Mat(5/2,0.0,0.0)
.
In this example we use a composite kernel represented as the sum of a Matern 5/2 ARD kernel and a Squared Exponential isotropic kernel. This is easily implemented using the +
symbol, or in the case of a product kernel, using the *
symbol (i.e. kern = Mat(5/2,[0.0,0.0],0.0) * SE(0.0,0.0)
).
mZero = MeanZero() # Zero mean function
kern = Matern(5/2,[0.0,0.0],0.0) + SE(0.0,0.0) # Sum kernel with Matern 5/2 ARD kernel
# with parameters [log(ℓ₁), log(ℓ₂)] = [0,0] and log(σ) = 0
# and Squared Exponential Iso kernel with
# parameters log(ℓ) = 0 and log(σ) = 0
Type: GaussianProcesses.SumKernel Type: GaussianProcesses.Mat52Ard, Params: [-0.0, -0.0, 0.0] Type: GaussianProcesses.SEIso, Params: [0.0, 0.0]
Fit the Gaussian process to the data using the prespecfied mean and covariance functions.
gp = GP(x,y,mZero,kern,-2.0) # Fit the GP
GP Exact object: Dim = 2 Number of observations = 50 Mean function: Type: GaussianProcesses.MeanZero, Params: Float64[] Kernel: Type: GaussianProcesses.SumKernel Type: GaussianProcesses.Mat52Ard, Params: [-0.0, -0.0, 0.0] Type: GaussianProcesses.SEIso, Params: [0.0, 0.0] Input observations = [2.4074 4.74658 … 5.90308 3.64732; 1.10184 1.46959 … 1.66664 0.937862] Output observations = [0.634667, -0.965483, 0.823632, 0.0182852, 0.63505, 0.469851, -0.923228, -0.254079, -0.585157, 0.0903145 … -0.250312, 0.278147, 0.396072, -0.26436, 0.943786, 0.84831, -0.691183, 0.51154, -0.323054, -0.343967] Variance of observation noise = 0.01831563888873418 Marginal Log-Likelihood = -24.734
Using the Optim package we have the option to choose from a range of optimize functions including conjugate gradients. It is also possible to fix the hyperparameters in either the mean function, kernel function or observation noise, by settting them to false in optimize!
(e.g. optimize!(...,mean=false)
).
optimize!(gp) # Optimize the hyperparameters
Dict(:mean=>true,:kern=>true,:noise=>true)
Results of Optimization Algorithm * Algorithm: L-BFGS * Starting Point: [-2.0,-0.0,-0.0,0.0,0.0,0.0] * Minimizer: [-4.280711367180729,1.8511934768970428, ...] * Minimum: -5.508004e+01 * Iterations: 39 * Convergence: false * |x - x'| < 1.0e-32: false * |f(x) - f(x')| / |f(x)| < 1.0e-32: false * |g(x)| < 1.0e-08: false * f(x) > f(x'): true * Reached Maximum Number of Iterations: false * Objective Function Calls: 160 * Gradient Calls: 160
A range of plotting options are availbe through the Plots.jl package.
plot(contour(gp) ,heatmap(gp))