{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Simple GP Regression "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Gaussian processes are a powerful tool for nonlinear regression models. \n",
"\n",
"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$.\n",
"\n",
"The response variables $\\mathbf{y}$ are assumed to dependent on the predictors $\\mathbf{X}$,\n",
"\n",
"\\begin{equation}\n",
"y_i \\sim \\mathcal{N}(f(\\mathbf{x}_i),\\sigma^2), \\ i=1,\\ldots,n, \n",
"\\end{equation}\n",
"\n",
"where $f$ is a mapping function. Treating $f$ as a random function, we assume that the distribution over $f$ is a **Gaussian process**,\n",
"$$\n",
"f \\sim \\mathcal{GP}(m(\\mathbf{x}),k(\\mathbf{x},\\mathbf{x}')),\n",
"$$\n",
"where $m(\\cdot)$ and $k(\\cdot,\\cdot)$ are the mean and kernel functions respectively."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1D regression example"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start by simulating some data"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using GaussianProcesses\n",
"\n",
"srand(20140430)\n",
"# Training data\n",
"n=10; #number of training points\n",
"x = 2π * rand(n); #predictors\n",
"y = sin.(x) + 0.05*randn(n); #regressors"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first step in modelling with Gaussian Processes is to choose mean functions and kernels which describe the process. \n",
"\n",
"**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."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"GP Exact object:\n",
" Dim = 1\n",
" Number of observations = 10\n",
" Mean function:\n",
" Type: GaussianProcesses.MeanZero, Params: Float64[]\n",
" Kernel:\n",
" Type: GaussianProcesses.SEIso, Params: [0.0, 0.0]\n",
" Input observations = \n",
"[4.85461 5.17653 … 1.99412 3.45676]\n",
" Output observations = [-0.967293, -1.00705, -1.0904, 0.881121, -0.333213, -0.976965, 0.915934, 0.736218, 0.950849, -0.306432]\n",
" Variance of observation noise = 0.1353352832366127\n",
" Marginal Log-Likelihood = -6.335"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#Select mean and covariance function\n",
"mZero = MeanZero() #Zero mean function\n",
"kern = SE(0.0,0.0) #Sqaured exponential kernel (note that hyperparameters are on the log scale)\n",
"\n",
"logObsNoise = -1.0 # log standard deviation of observation noise (this is optional)\n",
"gp = GP(x,y,mZero,kern,logObsNoise) #Fit the GP"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"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,\n",
"\\begin{equation}\n",
" 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}),\n",
"\\end{equation}\n",
"where \n",
"\\begin{eqnarray}\n",
"\\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} \\\n",
"\\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).\n",
"\\end{eqnarray}\n",
"\n",
"**Note** you can use the `predict_f` function to predict the latent function $\\mathbf{f}^\\ast$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"μ, σ² = predict_y(gp,linspace(0,2π,100));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting GPs is straightforward and utilises the recipes approach to plotting from the [Plots.jl](https://github.com/JuliaPlots/Plots.jl) package. More information about plotting GPs and the available functionality can be found in this [notebook](https://github.com/STOR-i/GaussianProcesses.jl/blob/master/notebooks/Plotting%20GPs.ipynb).\n",
"\n",
"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 `;`."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots #Load Plots.jl package\n",
"\n",
"plot(gp; xlabel=\"x\", ylabel=\"y\", title=\"Gaussian process\", legend=false, fmt=:png) # Plot the GP"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The hyperparameters are optimized using the [Optim.jl](https://github.com/JuliaOpt/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. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Dict(:mean=>true,:kern=>true,:noise=>true)\n"
]
},
{
"data": {
"text/plain": [
"Results of Optimization Algorithm\n",
" * Algorithm: BFGS\n",
" * Starting Point: [-1.0,0.0,0.0]\n",
" * Minimizer: [-2.99484893977448,0.4636861184637202, ...]\n",
" * Minimum: -3.275745e+00\n",
" * Iterations: 9\n",
" * Convergence: true\n",
" * |x - x'| < 1.0e-32: false\n",
" * |f(x) - f(x')| / |f(x)| < 1.0e-32: true\n",
" * |g(x)| < 1.0e-08: true\n",
" * f(x) > f(x'): false\n",
" * Reached Maximum Number of Iterations: false\n",
" * Objective Function Calls: 42\n",
" * Gradient Calls: 42"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"optimize!(gp; method=Optim.BFGS()) # Optimise the hyperparameters"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(gp) #Plot the GP after the hyperparameters have been optimised "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"BasicMCJob:\n",
" Variable [1]: p (BasicContMuvParameter)\n",
" GenericModel: 1 variables, 0 dependencies (directed graph)\n",
" HMC sampler: number of leaps = 10, leap step = 0.1\n",
" VanillaMCTuner: period = 100, verbose = false\n",
" BasicMCRange: number of steps = 1000, burnin = 0, thinning = 1"
]
},
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Distributions\n",
"\n",
"set_priors!(kern, [Normal(), Normal()]) # Uniform(0,1) distribution assumed by default if priors not specified\n",
"chain = mcmc(gp)\n",
"plot(chain', label=[\"Noise\", \"SE log length\", \"SE log scale\"])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Multi-dimensional regression"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"#Training data\n",
"d, n = 2, 50; #Dimension and number of observations\n",
"x = 2π * rand(d, n); #Predictors\n",
"y = vec(sin.(x[1,:]).*sin.(x[2,:])) + 0.05*rand(n); #Responses"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"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)`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)`)."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Type: GaussianProcesses.SumKernel\n",
" Type: GaussianProcesses.Mat52Ard, Params: [-0.0, -0.0, 0.0]\n",
" Type: GaussianProcesses.SEIso, Params: [0.0, 0.0]\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mZero = MeanZero() # Zero mean function\n",
"kern = Matern(5/2,[0.0,0.0],0.0) + SE(0.0,0.0) # Sum kernel with Matern 5/2 ARD kernel \n",
" # with parameters [log(ℓ₁), log(ℓ₂)] = [0,0] and log(σ) = 0\n",
" # and Squared Exponential Iso kernel with\n",
" # parameters log(ℓ) = 0 and log(σ) = 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Fit the Gaussian process to the data using the prespecfied mean and covariance functions."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"GP Exact object:\n",
" Dim = 2\n",
" Number of observations = 50\n",
" Mean function:\n",
" Type: GaussianProcesses.MeanZero, Params: Float64[]\n",
" Kernel:\n",
" Type: GaussianProcesses.SumKernel\n",
" Type: GaussianProcesses.Mat52Ard, Params: [-0.0, -0.0, 0.0]\n",
" Type: GaussianProcesses.SEIso, Params: [0.0, 0.0]\n",
" Input observations = \n",
"[2.4074 4.74658 … 5.90308 3.64732; 1.10184 1.46959 … 1.66664 0.937862]\n",
" 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]\n",
" Variance of observation noise = 0.01831563888873418\n",
" Marginal Log-Likelihood = -24.734"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gp = GP(x,y,mZero,kern,-2.0) # Fit the GP"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the [Optim](https://github.com/JuliaOpt/Optim.jl) 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)`)."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Dict(:mean=>true,:kern=>true,:noise=>true)\n"
]
},
{
"data": {
"text/plain": [
"Results of Optimization Algorithm\n",
" * Algorithm: L-BFGS\n",
" * Starting Point: [-2.0,-0.0,-0.0,0.0,0.0,0.0]\n",
" * Minimizer: [-4.280711367180729,1.8511934768970428, ...]\n",
" * Minimum: -5.508004e+01\n",
" * Iterations: 39\n",
" * Convergence: false\n",
" * |x - x'| < 1.0e-32: false\n",
" * |f(x) - f(x')| / |f(x)| < 1.0e-32: false\n",
" * |g(x)| < 1.0e-08: false\n",
" * f(x) > f(x'): true\n",
" * Reached Maximum Number of Iterations: false\n",
" * Objective Function Calls: 160\n",
" * Gradient Calls: 160"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"optimize!(gp) # Optimize the hyperparameters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A range of plotting options are availbe through the [Plots.jl](https://github.com/JuliaPlots/Plots.jl) package."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(contour(gp) ,heatmap(gp))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.6.0",
"language": "julia",
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.6.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}