# gpR: Gaussian processes in R
#
# Copyright (C) 2015 - 2016 Simon Dirmeier
#
# This file is part of gpR.
#
# gpR is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# gpR is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with gpR. If not, see <http://www.gnu.org/licenses/>.
#' @title Predict a new set of data points x.new using Gaussian process regression
#'
#' @noRd
#' @param object an \code{lvgpr.data} object
#' @return an \code{lvgpr.pred} object
#' \item{y.predict }{ the predicted y* values given the \code{lvgpr} object}
#' \item{mean }{ the posterior mean values}
#' \item{cov }{ the posterior covariance/kernel}
#' @references
#' Rasmussen C.E. and Williams C.K.I. (2006), \emph{Gaussian Processes for Machine Learning}, MIT Press \cr
#' \url{http://www.gaussianprocess.org/gpml/} \cr \cr
#' Barber D. (2015), \emph{Bayesian Reasoning and Machine Learning}, Cambridge University Press \cr
#' \url{http://web4.cs.ucl.ac.uk/staff/D.Barber/pmwiki/pmwiki.php?n=Brml.Online}
predict.gaussian <- function(obj, ...)
{
x.train <- obj@x.train
y.train <- obj@y.train
x.new <- obj@x.new
pars <- obj@pars
# small pseudo-count for numerical stability
C <- 1e-5 * diag(length(x.train))
# covariance of training inputs and new inputs
cov.new.train <- covariance.function(x1=x.new, x2=x.train, pars=pars)
# inverse covariance of training points
cov.inv.train.train <- solve(C + covariance.function(x1=x.train,
x2=x.train,
pars=pars))
# posterior mean value (see references: Barber p396, Equation~19.2.5;
# Rasmussen p16, Equation~2.23)
m <- cov.new.train %*% cov.inv.train.train %*% y.train
# posterior covariance of new inputs
# (see references: Barber p396, Equation~19.2.5;
# Rasmussen p16, Equation~2.24)
cov.new.new <- covariance.function(x1=x.new, x2=x.new, pars=pars) -
cov.new.train %*% cov.inv.train.train %*% t(cov.new.train)
# sample from the posterior Gaussian process to get y*
y <- sample.from.gp(x=x.new, m=m, K=cov.new.new)
ret <- list(y.predict=y, cov=cov.new.new, mean=m)
ret$call <- match.call()
class(ret) <- "lvgpr.pred"
ret
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.