## Function to predict new data points for a GP
#' Get the posterior distribution for f(X)
#'
#' @param gp A gaussianProcess object with data, kernel, etc.
#' @param X An n x d matrix of new observations
#' @return posterior A list with fields mean and covariance
#'
#' @export
predict.gaussianProcess <- function(gp, X) {
if(! is.matrix(X)){
stop("Input is not a matrix")
}
Knew <- gp$kernel(X, gp$data)
# get the posterior mean
means <- gp$meanFunc(X)
postMean <- means + Knew %*% gp$alpha
Kself <- gp$kernel(X, X)
Linv <- forwardsolve(gp$cholesky, diag(dim(gp$cholesky)[1]))
postCov <- Linv %*% t(Knew)
postCov <- (Knew %*% t(Linv)) %*% postCov
postCov <- Kself - postCov
# add prior noise
# postCov <- postCov + gp$noise.var * diag(dim(postCov)[1])
posterior <- list(mean=postMean, covariance=postCov)
return(posterior)
}
## Override the sample method to be able to have an S3 generic (funky stuff)
sample.default = sample
sample <- function(obj, ...) {
UseMethod("sample")
}
#' Sample from the posterior distribution for f(X) (or the posterior predictive)
#'
#' @param gp A gaussianProcess object with data, kernel, etc.
#' @param X An n x d matrix of new observations
#' @param n Number of samples
#' @param pred Boolean, whether to sample from the posterior predictive
#' @return An n dimensional vector of samples from the posterior (predictive)
#'
#' @export
sample.gaussianProcess <- function(gp, X, n, pred=F) {
# get the posterior of f
post <- predict(gp, X)
mu <- post$mean
Sigma <- post$covariance
if(pred) {
Sigma <- Sigma + gp$noise.var * diag(dim(Sigma)[1])
}
# sample
y <- MASS::mvrnorm(n, mu, Sigma)
return(y)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.