R/dynemu_mc.R

In dynemu: Emulation of Dynamic Simulators via One-Step-Ahead Approach

#' Predictive posterior computation via Monte Carlo (MC) approximation for one-steap-ahead Gaussian process (GP) emulators
#'
#' @description The function computes the predictive posterior distribution (mean and variance) for GP emulators
#' using MC method, given uncertain inputs.
#'
#' @details Given a trained set of GP emulators `dynGP` fitted using \code{dynemu_GP}, this function:
#' 1. Computes the MC-approximated predictive posterior mean and variance for each state variable.
#' 2. Incorporates input uncertainty by integrating over the input distribution via MC sampling.
#'
#' Unlike \code{dynemu_exact}, which provides a closed-form solution,
#' this function uses MC sampling to approximate the posterior.
#'
#' @references
#' H. Mohammadi, P. Challenor, and M. Goodfellow (2019). Emulating dynamic non-linear simulators using Gaussian processes.
#' \emph{Computational Statistics & Data Analysis}, 139, 178-196;
#' \doi{doi:10.1016/j.csda.2019.05.006}
#'
#' @param mean numeric vector or matrix specifying the mean of uncertain inputs. The number of columns must match the input dimension of `dynGP`.
#' @param var numeric vector or matrix specifying the variance of uncertain inputs. The number of columns must match the input dimension of `dynGP`.
#' @param dynGP list of trained GP emulators fitted by \code{dynemu_GP}, each corresponding to a state variable.
#' @param mc.sample a number of mc samples generated for MC approximation. Default is \code{1000}.
#'
#' @return A list containing:
#' \itemize{
#'   \item \code{predy}: A single-row matrix of predictive mean values, with each column for corresponding state variable.
#'   \item \code{predsig2}: A single-row matrix of predictive variance values, with each column for corresponding state variable.
#' }
#'
#' @importFrom MASS mvrnorm
#' @usage dynemu_mc(mean, var, dynGP, mc.sample)
#' @export
#'
#' @examples
#' library(lhs)
#' ### Lotka-Volterra equations ###
#' LVmod0D <- function(Time, State, Pars) {
#'   with(as.list(c(State, Pars)), {
#'     IngestC <- rI * P * C
#'     GrowthP <- rG * P * (1 - P/K)
#'     MortC <- rM * C
#'
#'     dP <- GrowthP - IngestC
#'     dC <- IngestC * AE - MortC
#'     return(list(c(dP, dC)))
#'   })
#' }
#'
#' ### Define parameters ###
#' pars <- c(rI = 1.5, rG = 1.5, rM = 2, AE = 1, K = 10)
#'
#' ### Define time sequence ###
#' times <- seq(0, 30, by = 0.01)
#'
#' ### Initial conditions ###
#' set.seed(1)
#' d <- 2
#' n <- 12*d
#' Design <- maximinLHS(n, d) * 5 # Generate LHS samples in range [0,5]
#' colnames(Design) <- c("P", "C")
#'
#' ### Fit GP emulators ###
#' fit.dyn <- dynemu_GP(LVmod0D, times, pars, Design)
#'
#' ### Define uncertain inputs ###
#' xmean <- c(P = 1, C = 2)
#' xvar <- c(P = 1e-5, C = 1e-5)
#'
#' ### Define number of MC samples ###
#' mc.sample <- 1000
#'
#' ### Compute the next point ###
#' dynemu_mc(xmean, xvar, fit.dyn, mc.sample)
#'

dynemu_mc <- function(mean, var, dynGP, mc.sample=1000){
  if (is.null(dim(mean))) mean <- matrix(mean, nrow = 1)
  if (is.null(dim(var))) var <- matrix(var, nrow = 1)
  colnames(mean) <- colnames(var) <- colnames(dynGP[[1]]$X)
  if (ncol(mean) != ncol(var) || ncol(mean) != length(dynGP)){
    stop("The number of columns in mean and var must match the number of GP emulators in dynGP. \n")
  }
  d <- length(dynGP)

  ### fit GPs ###
  pred_mc <- vector("list", 2)
  for (i in 1:2) { # For mu and sig2
    pred_mc[[i]] <- matrix(NA, nrow = 1, ncol = d)
  }

  for (i in 1:d) {
    xt.sample <- mvrnorm(mc.sample, mean, diag(c(var)))
    colnames(xt.sample) <- colnames(dynGP[[1]]$X)

    mc.pred <- pred.GP(dynGP[[i]], xt.sample)

    pred_mc[[1]][1,i] <- mean(mc.pred$mu)
    pred_mc[[2]][1,i] <- (mean(mc.pred$sig2) + var(mc.pred$mu))
  }
  colnames(pred_mc[[1]]) <- colnames(pred_mc[[2]]) <- colnames(dynGP[[1]]$X)

  return(list(predy=pred_mc[[1]], predsig2=pred_mc[[2]]))
}