# R/predictGP.R In ebenmichael/gaussianProcess: Fits Gaussian Processes for Regression

#### Documented in predict.gaussianProcesssample.gaussianProcess

```## 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
# 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)
}
```
ebenmichael/gaussianProcess documentation built on May 13, 2017, 10:58 a.m.