# Simple GP Regression 

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.

## 1D regression example

We start by simulating some data


```julia
using GaussianProcesses
using Random

Random.seed!(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.


```julia
#Select mean and covariance function
mZero = MeanZero()                   #Zero mean function
kern = SE(0.0,0.0)                   #Squared 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: MeanZero, Params: Float64[]
      Kernel:
        Type: SEIso{Float64}, Params: [0.0, 0.0]
      Input observations = 
    [4.854610892030431 5.176527683588912 … 1.9941157865846477 3.4567590405289352]
      Output observations = [-0.9672931901680022, -1.0070469610508248, -1.0903962829715461, 0.8811208605233783, -0.3332131916909575, -0.9769651482106089, 0.9159341732111771, 0.7362180710829012, 0.9508490881368652, -0.3064317818816142]
      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 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$.


```julia
μ, σ² = predict_y(gp,range(0,stop=2π,length=100));
```

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).

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 `;`.


```julia
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](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. 


```julia
using Optim
optimize!(gp; method=ConjugateGradient())   # Optimise the hyperparameters
```




     * Status: failure (objective increased between iterations) (line search failed)
    
     * Candidate solution
        Minimizer: [-2.99e+00, 4.64e-01, -2.29e-01]
        Minimum:   -3.275745e+00
    
     * Found with
        Algorithm:     Conjugate Gradient
        Initial Point: [-1.00e+00, 0.00e+00, 0.00e+00]
    
     * Convergence measures
        |x - x'|               = 4.62e-09 ≰ 0.0e+00
        |x - x'|/|x'|          = 1.54e-09 ≰ 0.0e+00
        |f(x) - f(x')|         = 1.35e-13 ≰ 0.0e+00
        |f(x) - f(x')|/|f(x')| = 4.12e-14 ≰ 0.0e+00
        |g(x)|                 = 5.69e-08 ≰ 1.0e-08
    
     * Work counters
        Seconds run:   0  (vs limit Inf)
        Iterations:    27
        f(x) calls:    66
        ∇f(x) calls:   41





```julia
plot(gp; legend=false, fmt=:png)   #Plot the GP after the hyperparameters have been optimised 
```

By default all hyperparameters or optimized. But the function `optimize!` allows also to force hyperparameters to remain constant or optimize them in a box.


```julia
optimize!(gp; kern = false)   # Don't optimize kernel hyperparameters
optimize!(gp; kernbounds = [[-1, -1], [1, 1]]) # Optimize the kernel parameters in a box with lower bounds [-1, -1] and upper bounds [1, 1]
```




     * Status: success
    
     * Candidate solution
        Minimizer: [-2.99e+00, 4.64e-01, -2.29e-01]
        Minimum:   -3.275745e+00
    
     * Found with
        Algorithm:     Fminbox with L-BFGS
        Initial Point: [-2.99e+00, 4.64e-01, -2.29e-01]
    
     * Convergence measures
        |x - x'|               = 0.00e+00 ≤ 0.0e+00
        |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
        |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
        |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
        |g(x)|                 = 6.32e-08 ≰ 1.0e-08
    
     * Work counters
        Seconds run:   1  (vs limit Inf)
        Iterations:    1
        f(x) calls:    3
        ∇f(x) calls:   3




MCMC, specifically a Hamiltonian Monte Carlo in this instance, 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.


```julia
using Distributions

set_priors!(kern, [Normal(), Normal()]) # Uniform(0,1) distribution assumed by default if priors not specified
@time chain = mcmc(gp, nIter=5000);
```

    Number of iterations = 5000, Thinning = 1, Burn-in = 1 
    Step size = 0.100000, Average number of leapfrog steps = 10.131000 
    Number of function calls: 50656
    Acceptance rate: 0.975200 
      2.699465 seconds (3.62 M allocations: 218.559 MiB, 2.83% gc time)


There is additional support for inference to be done in a fully Bayesian fashion through the use of an [elliptical slice sampler](http://proceedings.mlr.press/v9/murray10a/murray10a.pdf). While a HMC sampler can often be shown to be highly efficient, the sampler's efficiency can be highly dependent upon good initial choice of the sampler's hyperparameters. Conversely, ESS has no free parameters and is designed to be highly efficient in tightly correlated Gaussian posteriors: a geometry commonly found in Gaussian process models. Currently, inference via an ESS is only supported when the likelihood is Gaussian. 

For advice on how to fit Gaussian processes with non-Gaussian data, see our notebook on [Poisson regression](https://github.com/STOR-i/GaussianProcesses.jl/blob/master/notebooks/Poisson%20regression.ipynb) or [classification](https://github.com/STOR-i/GaussianProcesses.jl/blob/master/notebooks/Classification.ipynb).


```julia
set_priors!(gp.logNoise, [Normal()])
@time ess_chain = ess(gp, nIter=5000);
```

    Number of iterations = 5000, Thinning = 1, Burn-in = 1 
    Number of function calls: 28804
    Acceptance rate: 0.210049 
      0.353048 seconds (1.59 M allocations: 95.275 MiB, 12.60% gc time)



```julia
p1 = plot(chain', label=["Noise" "SE log length" "SE log scale"]; fmt=:png, size = (1200, 600), title="Hamiltonian Monte-Carlo", ylims=(-4, 2.5))
p2 = plot(ess_chain', label=["Noise" "SE log length" "SE log scale"]; fmt=:png, size = (1200, 600), title="Elliptical Slice Sampler", ylims=(-4, 2.5))
plot(p1, p2, xlabel="Iteration", ylabel="Parameter Value")
```

As can be seen in the timing outputs, Hamiltonian MC is considerably slower than elliptical slice sampling, however, this increase in computational cost is accompanied by chains with more desirable mixing.

## Multi-dimensional regression

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.


```julia
#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)`).


```julia
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: SumKernel{Mat52Ard{Float64},SEIso{Float64}}
      Type: Mat52Ard{Float64}, Params: [-0.0, -0.0, 0.0]  Type: SEIso{Float64}, Params: [0.0, 0.0]



Fit the Gaussian process to the data using the prespecfied mean and covariance functions.


```julia
gp = GP(x,y,mZero,kern,-2.0)          # Fit the GP
```




    GP Exact object:
      Dim = 2
      Number of observations = 50
      Mean function:
        Type: MeanZero, Params: Float64[]
      Kernel:
        Type: SumKernel{Mat52Ard{Float64},SEIso{Float64}}
          Type: Mat52Ard{Float64}, Params: [-0.0, -0.0, 0.0]      Type: SEIso{Float64}, Params: [0.0, 0.0]
      Input observations = 
    [3.05977 4.74752 … 2.82127 5.38224; 2.02102 4.27258 … 6.13114 1.56497]
      Output observations = [0.08509, 0.924505, 0.275745, -0.448035, -0.784758, -0.316803, -0.823483, -0.886726, 0.0059149, 0.414951  …  -0.413905, -0.347505, 0.46108, -0.204102, -0.538689, 0.554203, -0.874479, -0.0506017, -0.0215167, -0.745944]
      Variance of observation noise = 0.01831563888873418
      Marginal Log-Likelihood = -29.547



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!(...,domean=false)`).


```julia
optimize!(gp)                         # Optimize the hyperparameters
```




    Results of Optimization Algorithm
     * Algorithm: L-BFGS
     * Starting Point: [-2.0,-0.0,-0.0,0.0,0.0,0.0]
     * Minimizer: [-4.277211995773057,5.536664077544042, ...]
     * Minimum: -5.023830e+01
     * Iterations: 34
     * Convergence: false
       * |x - x'| ≤ 0.0e+00: false 
         |x - x'| = 4.74e-02 
       * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: false
         |f(x) - f(x')| = -1.36e-14 |f(x)|
       * |g(x)| ≤ 1.0e-08: false 
         |g(x)| = 1.87e-05 
       * Stopped by an increasing objective: true
       * Reached Maximum Number of Iterations: false
     * Objective Calls: 115
     * Gradient Calls: 115



A range of plotting options are availbe through the [Plots.jl](https://github.com/JuliaPlots/Plots.jl) package.


```julia
plot(contour(gp) ,heatmap(gp); fmt=:png)
```

1e15bd364e4ea5fe143c29699d21bef4f7fd5130
ipynb
Jupyter Notebook
notebooks/Regression.ipynb
kykim0/GaussianProcesses.jl
0b790388359a87a91d219b69d53ff475a9eae4f5
MIT
2015-07-31T21:32:34.000Z
2022-03-03T07:08:58.000Z
notebooks/Regression.ipynb
kykim0/GaussianProcesses.jl
0b790388359a87a91d219b69d53ff475a9eae4f5
MIT
2015-05-28T12:04:03.000Z
2022-02-21T06:48:23.000Z
notebooks/Regression.ipynb
kykim0/GaussianProcesses.jl
0b790388359a87a91d219b69d53ff475a9eae4f5
MIT
2015-05-31T12:44:26.000Z
2022-03-29T05:21:38.000Z
Qwen/Qwen-72B
1. YES
2. YES