R/GPSample.R In OakleyJ/MUCM: Gaussian process emulator methods based on the MUCM toolkit

#' @title GP Sample

#' @description Function to sample from posterior Gaussian process

#' @param n the number of samples required
#' @param mean a vector giving the means of the variables
#' @param var a positive-definite symmetric matrix specifying the covariance matrix of the variables

#' @author Jeremy Oakley.
#'
#' @importFrom mvtnorm rmvnorm
#' @export
GPSample <- function(n, mean, var) {

# Use pivoted Cholesky to decide which outputs should be sampled
# Use the mean of the remaining outputs, conditional on the sampled values
U <- chol(var, pivot = T)
index.in <- attr(U, "pivot")[1:attr(U, "rank")]  # index of which points will be sampled
index.out <- setdiff(1:length(mean), index.in)  # index of which points will not be sampled

V1 <- var[index.in, index.in]
m1 <- mean[index.in, ]
m2 <- mean[index.out, ]

new.y1 <- t(rmvnorm(n, m1, V1))

# Calculate conditional mean of unsampled points, given the sampled points
new.y2 <- matrix(t(m2), length(m2), n) + var[index.out, index.in] %*% solve(V1, new.y1 - matrix(t(m1), length(m1), n))

f.sample <- matrix(0, length(mean), n)
f.sample[index.in, ] <- new.y1
f.sample[index.out, ] <- new.y2
f.sample
}
OakleyJ/MUCM documentation built on May 7, 2019, 9:01 p.m.