The 'gplite' package.

Description

gplite implements some of the most common Gaussian process (GP) models. The package offers tools for integrating out the latent values analytically using Laplace or expectation propagation (EP) approximation and for estimating the hyperparameters based on maximizing the (approximate) marginal likelihood or posterior. The package also implements some common sparse approximations for larger datasets.

Functions

Here's a list of the most important functions:

gp_init

Set up the GP model.

cf, lik, method, approx

Choose the covariance functions, likelihood (observation model), type of the GP (full or some sparse approximation) and the latent function approximation method (Laplace, EP).

gp_optim, gp_fit

Optimize the model hyperparameters, or just fit the model with the current hyperparameter values.

gp_pred, gp_draw

Make predictions with the fitted model. Can also be used before fitting to obtain prior predictive distribution or draws.

gp_loo, gp_compare

Model assessment and comparison using leave-one-out (LOO) cross-validation.

Approximations to the posterior of the latent values

Description

Functions for initializing the approximation for the latent values, which can then be passed to gp_init. The supported methods are:

approx_laplace

Laplace's method, that is, based on local second order approximation to the log likelihood. For Gaussian likelihood, this means exact inference (no approximation).

approx_ep

Expectation propagation, EP. Approximates the likelihood by introducing Gaussian pseudo-data so that the posterior marginals match to the so called tilted distributions (leave-one-out posterior times the true likelihood factor) as closely as possible. Typically more accurate than Laplace, but slower.

Usage

approx_laplace(maxiter = 30, tol = 1e-04)

approx_ep(damping = 0.9, quad_order = 11, maxiter = 100)

Arguments

Value

The approximation object.

References

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian processes for machine learning. MIT Press.

Examples

# Basic usage
gp <- gp_init(
  cfs = cf_sexp(),
  lik = lik_bernoulli(),
  method = method_fitc(num_inducing = 100),
  approx = approx_ep()
)

Initialize covariance function

Description

Functions for initializing the covariance functions which can then be passed to gp_init. See section Details for explanation of what covariance function is what.

Usage

cf_const(magn = 1, prior_magn = prior_logunif())

cf_lin(vars = NULL, magn = 1, prior_magn = prior_logunif(), normalize = FALSE)

cf_sexp(
  vars = NULL,
  lscale = 0.3,
  magn = 1,
  prior_lscale = prior_logunif(),
  prior_magn = prior_logunif(),
  normalize = FALSE
)

cf_matern32(
  vars = NULL,
  lscale = 0.3,
  magn = 1,
  prior_lscale = prior_logunif(),
  prior_magn = prior_logunif(),
  normalize = FALSE
)

cf_matern52(
  vars = NULL,
  lscale = 0.3,
  magn = 1,
  prior_lscale = prior_logunif(),
  prior_magn = prior_logunif(),
  normalize = FALSE
)

cf_nn(
  vars = NULL,
  sigma0 = 1,
  sigma = 3,
  magn = 1,
  prior_sigma0 = prior_half_t(),
  prior_sigma = prior_half_t(),
  prior_magn = prior_logunif(),
  normalize = TRUE
)

cf_periodic(
  vars = NULL,
  period = 1,
  cf_base = cf_sexp(),
  prior_period = prior_logunif()
)

cf_prod(...)

## S3 method for class 'cf'
cf1 * cf2

Arguments

Details

The supported covariance functions are (see Rasmussen and Williams, 2006):

cf_const

Constant covariance function. Can be used to model the intercept.

cf_lin

Linear covariance function. Produces linear functions.

cf_sexp

Squared exponential (or exponentiated quadratic, or Gaussian) covariance function.

cf_matern32

Matern nu=3/2 covariance function.

cf_matern52

Matern nu=5/2 covariance function.

cf_nn

Neural network covariance function.

cf_periodic

Periodic covariance function. The periodicity is achieved by mapping the original inputs through sine and cosine functions, and then applying the base kernel in this new space.

cf_prod

Product of two or more covariance functions.

Value

The covariance function object.

References

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian processes for machine learning. MIT Press.

Examples

# Generate some toy data
set.seed(1242)
n <- 50
x <- matrix(rnorm(n * 3), nrow = n)
f <- sin(x[, 1]) + 0.5 * x[, 2]^2 + x[, 3]
y <- f + 0.5 * rnorm(n)
x <- data.frame(x1 = x[, 1], x2 = x[, 2], x3 = x[, 3])

# Basic usage (single covariance function)
cf <- cf_sexp()
lik <- lik_gaussian()
gp <- gp_init(cf, lik)
gp <- gp_optim(gp, x, y)
plot(gp_pred(gp, x)$mean, y)

# More than one covariance function; one for x1 and x2, and another one for x3
cf1 <- cf_sexp(c("x1", "x2"))
cf2 <- cf_lin("x3")
cfs <- list(cf1, cf2)
lik <- lik_gaussian()
gp <- gp_init(cfs, lik)
gp <- gp_optim(gp, x, y, maxiter = 500)
plot(gp_pred(gp, x)$mean, y)
plot(x[, 3], gp_pred(gp, x, cfind = 2)$mean) # plot effect w.r.t x3 only

Make predictions with a GP model

Description

Function gp_pred can be used to make analytic predictions for the latent function values at test points, whereas gp_draw can be used to draw from the predictive distribution (or from the prior if the GP has not been fitted yet.)

Usage

gp_draw(
  gp,
  xnew,
  draws = NULL,
  transform = TRUE,
  target = FALSE,
  marginal = FALSE,
  cfind = NULL,
  jitter = NULL,
  seed = NULL,
  ...
)

gp_pred(
  gp,
  xnew,
  var = FALSE,
  quantiles = NULL,
  transform = FALSE,
  cfind = NULL,
  jitter = NULL,
  quad_order = 15,
  ...
)

Arguments

Value

gp_pred returns a list with fields giving the predictive mean, variance and quantiles (the last two are computed only if requested). gp_draw returns an N-by-draws matrix of random draws from the predictive distribution, where N is the number of test points.

References

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian processes for machine learning. MIT Press.

Examples

# Generate some toy data
set.seed(1242)
n <- 50
x <- matrix(rnorm(n * 3), nrow = n)
f <- sin(x[, 1]) + 0.5 * x[, 2]^2 + x[, 3]
y <- f + 0.5 * rnorm(n)
x <- data.frame(x1 = x[, 1], x2 = x[, 2], x3 = x[, 3])

# More than one covariance function; one for x1 and x2, and another one for x3
cf1 <- cf_nn(c("x1", "x2"), prior_sigma0 = prior_half_t(df = 4, scale = 2))
cf2 <- cf_sexp("x3")
cfs <- list(cf1, cf2)
lik <- lik_gaussian()
gp <- gp_init(cfs, lik)
gp <- gp_optim(gp, x, y, maxiter = 500)

# plot the predictions with respect to x1, when x2 = x3 = 0
xt <- cbind(x1 = seq(-3, 3, len = 50), x2 = 0, x3 = 0)
pred <- gp_pred(gp, xt)
plot(xt[, "x1"], pred$mean, type = "l")

# draw from the predictive distribution
xt <- cbind(x1 = seq(-3, 3, len = 50), x2 = 0, x3 = 0)
draws <- gp_draw(gp, xt, draws = 100)
plot(xt[, "x1"], draws[, 1], type = "l")
for (i in 2:50) {
  lines(xt[, "x1"], draws[, i])
}

# plot effect with respect to x3 only
xt <- cbind("x3" = seq(-3, 3, len = 50))
pred <- gp_pred(gp, xt, cfind = 2)
plot(xt, pred$mean, type = "l")

Energy of a GP model

Description

Returns the energy (negative log marginal likelihood) of a fitted GP model with the current hyperparameters. The result is exact for the Gaussian likelihood and dependent on the approx for other cases.

Usage

gp_energy(gp, include_prior = TRUE)

Arguments

Value

The energy value (negative log marginal likelihood).

References

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian processes for machine learning. MIT Press.

Examples

# Generate some toy data
set.seed(1242)
n <- 500
x <- matrix(rnorm(n * 3), nrow = n)
f <- sin(x[, 1]) + 0.5 * x[, 2]^2 + x[, 3]
y <- f + 0.5 * rnorm(n)
x <- data.frame(x1 = x[, 1], x2 = x[, 2], x3 = x[, 3])

# Basic usage
gp <- gp_init(cf_sexp(), lik_gaussian())
gp <- gp_fit(gp, x, y)
e <- gp_energy(gp)

Fit a GP model

Description

Function gp_fit fits a GP model with the current hyperparameters. Notice that this function does not optimize the hyperparameters in any way, but only finds the analytical posterior approximation (depending on chosen approx) for the latent values with the current hyperparameters. For optimizing the hyperparameter values, see gp_optim.

Usage

gp_fit(gp, x, y, trials = NULL, offset = NULL, jitter = NULL, ...)

Arguments

Value

An updated GP model object.

References

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian processes for machine learning. MIT Press.

Examples

# Generate some toy data
set.seed(32004)
n <- 150
sigma <- 0.1
x <- rnorm(n)
ycont <- sin(3 * x) * exp(-abs(x)) + rnorm(n) * sigma
y <- rep(0, n)
y[ycont > 0] <- 1
trials <- rep(1, n)

# Fit the model using Laplace approximation (with the specified hyperparameters)
cf <- cf_sexp(lscale = 0.3, magn = 3)
gp <- gp_init(cf, lik_binomial())
gp <- gp_fit(gp, x, y, trials = trials)

Initialize a GP model

Description

Initializes a GP model with given covariance function(s) and likelihood. The model can then be fitted using gp_fit. For hyperparameter optimization, see gp_optim

Usage

gp_init(
  cfs = cf_sexp(),
  lik = lik_gaussian(),
  method = method_full(),
  approx = approx_laplace()
)

Arguments

Value

A GP model object that can be passed to other functions, for example when optimizing the hyperparameters or making predictions.

References

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian processes for machine learning. MIT Press.

Examples

# Full exact GP with Gaussian likelihood
gp <- gp_init(
  cfs = cf_sexp(),
  lik = lik_gaussian(),
  method = method_full()
)

# Binary classification model with EP approximation for the latent values
# and FITC sparse approximation to facilitate large datasets
gp <- gp_init(
  cfs = cf_sexp(),
  lik = lik_bernoulli(),
  approx = approx_ep(),
  method = method_fitc(num_inducing = 100)
)

Model assessment and comparison

Description

Function gp_loo computes the approximate leave-one-out (LOO) cross-validation statistics for the given GP model with the current hyperparameters. Function gp_compare estimates the difference in the expected predictive accuracy of two or more GP models given their LOO statistics.

Usage

gp_loo(
  gp,
  x,
  y,
  quadrature = TRUE,
  quad_order = 11,
  draws = 4000,
  jitter = NULL,
  seed = NULL,
  ...
)

gp_compare(..., ref = NULL)

Arguments

Value

gp_loo returns a list with LOO statistics. gp_compare returns a matrix with comparison statistics (LOO differences and stardard errors in the estimates).

References

Vehtari A., Mononen T., Tolvanen V., Sivula T. and Winther O. (2016). Bayesian Leave-One-Out Cross-Validation Approximations for Gaussian Latent Variable Models. Journal of Machine Learning Research 17(103):1-38.

Examples

# Generate some toy data
set.seed(32004)
n <- 50
sigma <- 0.1
x <- rnorm(n)
ycont <- sin(3 * x) * exp(-abs(x)) + rnorm(n) * sigma
y <- rep(0, n)
y[ycont > 0] <- 1
trials <- rep(1, n)

# Set up two models
gp1 <- gp_init(cf_sexp(), lik_binomial())
gp2 <- gp_init(cf_matern32(), lik_binomial())

# Optimize
gp1 <- gp_optim(gp1, x, y, trials = trials)
gp2 <- gp_optim(gp2, x, y, trials = trials)

# Compare
loo1 <- gp_loo(gp1, x, y, trials = trials)
loo2 <- gp_loo(gp2, x, y, trials = trials)
gp_compare(loo1, loo2)

Optimize hyperparameters of a GP model

Description

This function can be used to optimize the hyperparameters of the model to the maximum marginal likelihood (or maximum marginal posterior if priors are used), using Nelder-Mead algorithm.

Usage

gp_optim(
  gp,
  x,
  y,
  tol = 1e-04,
  tol_param = 0.1,
  maxiter = 500,
  restarts = 1,
  verbose = TRUE,
  warnings = TRUE,
  ...
)

Arguments

Value

An updated GP model object.

References

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian processes for machine learning. MIT Press.

Examples

# Generate some toy data
set.seed(1242)
n <- 50
x <- matrix(rnorm(n * 3), nrow = n)
f <- sin(x[, 1]) + 0.5 * x[, 2]^2 + x[, 3]
y <- f + 0.5 * rnorm(n)
x <- data.frame(x1 = x[, 1], x2 = x[, 2], x3 = x[, 3])

# Basic usage
cf <- cf_sexp()
lik <- lik_gaussian()
gp <- gp_init(cf, lik)
gp <- gp_optim(gp, x, y)

Save and load a GP model

Description

Convenience functions for saving and loading GP models.

Usage

gp_save(gp, filename)

gp_load(filename)

Arguments

Value

gp_load returns the loaded GP model object.

Examples

gp <- gp_init()

# fit the model (skipped here)

# save the model
filename <- file.path(tempdir(), 'gp.rda')
gp_save(gp, filename)

# load the model and remove the file
gp <- gp_load(filename)
file.remove(filename)

Initialize likelihood

Description

Functions for initializing the likelihood (observation model) which can then be passed to gp_init.

Usage

lik_gaussian(sigma = 0.5, prior_sigma = prior_logunif())

lik_bernoulli(link = "logit")

lik_binomial(link = "logit")

lik_betabinom(link = "logit", phi = 1, prior_phi = prior_logunif())

lik_poisson(link = "log")

Arguments

Details

The supported likelihoods are:

lik_gaussian

Gaussian likelihood. Has no links (uses identity link).

lik_bernoulli

Bernoulli likelihood. Possible links: 'logit' or 'probit'.

lik_binomial

Binomial likelihood. Possible links: 'logit' or 'probit'.

lik_betabinom

Beta binomial likelihood. Possible links: 'logit' or 'probit'.

lik_poisson

Poisson likelihood. Possible links: 'log'.

Value

The likelihood object.

Examples

# Basic usage
cf <- cf_sexp()
lik <- lik_binomial()
gp <- gp_init(cf, lik)

Initialize method or type of the model

Description

Functions for initializing the method or type of the model, which can then be passed to gp_init. The supported methods are:

method_full

Full GP, so full exact covariance function is used, meaning that the inference will be for the n latent function values (fitting time scales cubicly in n).

method_fitc

Fully independent training (and test) conditional, or FITC, approximation (see Quiñonero-Candela and Rasmussen, 2005; Snelson and Ghahramani, 2006). The fitting time scales O(n*m^2), where n is the number of data points and m the number of inducing points num_inducing. The inducing point locations are chosen using the k-means algorithm.

method_rf

Random features, that is, linearized GP. Uses random features (or basis functions) for approximating the covariance function, which means the inference time scales cubicly in the number of approximating basis functions num_basis. For stationary covariance functions random Fourier features (Rahimi and Recht, 2007) is used, and for non-stationary kernels using case specific method when possible (for example, drawing the hidden layer parameters randomly for cf_nn). For cf_const and cf_lin this means using standard linear model, and the inference is performed on the weight space (not in the function space). Thus if the model is linear (only cf_const and cf_lin are used), this will give a potentially huge speed-up if the number of features is considerably smaller than the number of data points.

Usage

method_full()

method_fitc(
  inducing = NULL,
  num_inducing = 100,
  bin_along = NULL,
  bin_count = 10,
  seed = 12345
)

method_rf(num_basis = 400, seed = 12345)

Arguments

Value

The method object.

References

Rahimi, A. and Recht, B. (2008). Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems 20.

Quiñonero-Candela, J. and Rasmussen, C. E (2005). A unifying view of sparse approximate Gaussian process regression. Journal of Machine Learning Research 6:1939-1959.

Snelson, E. and Ghahramani, Z. (2006). Sparse Gaussian processes using pseudo-inputs. In Advances in Neural Information Processing Systems 18.

Examples

#' # Generate some toy data
# NOTE: this is so small dataset that in reality there would be no point
# use sparse approximation here; we use this small dataset only to make this
# example run fast
set.seed(1242)
n <- 50
x <- matrix(rnorm(n * 3), nrow = n)
f <- sin(x[, 1]) + 0.5 * x[, 2]^2 + x[, 3]
y <- f + 0.5 * rnorm(n)
x <- data.frame(x1 = x[, 1], x2 = x[, 2], x3 = x[, 3])

# Full exact GP with Gaussian likelihood
gp <- gp_init(cf_sexp())
gp <- gp_optim(gp, x, y)

# Approximate solution using random features (here we use a very small 
# number of random features only to make this example run fast)
gp <- gp_init(cf_sexp(), method = method_rf(num_basis = 30))
gp <- gp_optim(gp, x, y)

# Approximate solution using FITC (here we use a very small 
# number of incuding points only to make this example run fast)
gp <- gp_init(cf_sexp(), method = method_fitc(num_inducing = 10))
gp <- gp_optim(gp, x, y)

Get or set GP model parameters

Description

get_param returns the current hyperparameters of the GP model in a vector. set_param can be used to set the parameters. Note that these functions are intended mainly for internal usage, and there is typically no need to use these functions directly but instead create a new GP model using gp_init.

Usage

get_param(object, ...)

set_param(object, param, ...)

Arguments

Value

get_param returns the current hyperparameters and set_param the GP model structure with the new parameter values.

Examples

# Set up some model
gp <- gp_init(cf = cf_sexp(), lik = lik_gaussian())

# print out to see the parameter ordering
param <- get_param(gp)
print(param)

# set some new values
param_new <- log(c(0.1, 0.8, 0.3))
names(param_new) <- names(param)
gp <- set_param(gp, param_new)

# check the result
print(get_param(gp))

Initialize prior for hyperparameter

Description

Functions for initializing hyperparameter priors which can then be passed to gp_init. See section Details for the prior explanations.

Usage

prior_fixed()

prior_logunif()

prior_lognormal(loc = 0, scale = 1)

prior_half_t(df = 1, scale = 1)

Arguments

Details

The supported priors are:

prior_fixed

The hyperparameter is fixed to its initial value, and is not optimized by gp_optim.

prior_logunif

Improper uniform prior on the log of the parameter.

prior_lognormal

Log-normal prior (Gaussian prior on the logarithm of the parameter).

prior_half_t

Half Student-t prior for a positive parameter.

Value

The hyperprior object.

References

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian processes for machine learning. MIT Press.

Examples

# Quasi-periodic covariance function, with fixed period
cf1 <- cf_periodic(
  period = 5,
  prior_period = prior_fixed(),
  cf_base = cf_sexp(lscale = 2)
)
cf2 <- cf_sexp(lscale = 40)
cf <- cf1 * cf2
gp <- gp_init(cf)

# draw from the prior
set.seed(104930)
xt <- seq(-10, 10, len = 500)
plot(xt, gp_draw(gp, xt), type = "l")