R/dynemu_exact.R

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

#' Predictive posterior computation via exact closed-form expression for one-steap-ahead Gaussian process (GP) emulators
#'
#' @description The function computes the predictive posterior distribution (mean and variance) for GP emulators
#' using closed-form expression, given uncertain inputs.
#'
#' @details Given a trained set of GP emulators `dynGP` fitted using \code{dynemu_GP}, this function:
#' 1. Computes the closed-form predictive posterior mean and variance for each state variable.
#' 2. Incorporates input uncertainty by integrating over the input distribution via exact computation.
#'
#' The computation follows a closed-form approach, leveraging the posterior distributions of Linked GP.
#'
#' @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.
#'
#' @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.
#' }
#'
#' @usage dynemu_exact(mean, var, dynGP)
#' @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)
#'
#' ### Compute the next point ###
#' dynemu_exact(xmean, xvar, fit.dyn)
#'

dynemu_exact <- function(mean, var, dynGP){
  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)
  n <- length(dynGP[[1]]$y)
  X <- dynGP[[1]]$X

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

  for (i in 1:d) {
    ### calculate the closed form ###
    const.mean <- cbind(1, X) %*% dynGP[[i]]$beta
    theta <- dynGP[[i]]$theta
    tau2hat <- dynGP[[i]]$tau2hat
    test.mean <- cbind(1, mean) %*% dynGP[[i]]$beta

    Ci <- dynGP[[i]]$Ki
    a <- Ci %*% (dynGP[[i]]$y - const.mean)

    # component calculation
    predy_component <- rep(1,n)
    for (j in 1:d) {
      diff_j <- mean[j] - X[,j]  # Vector of differences
      predy_component <- predy_component * (
        1 / sqrt(1 + 2 * var[j] / theta[j]) *
          exp(-(diff_j^2) / (theta[j] + 2 * var[j]))
      )
    }
    predsig2_component <- matrix(1,n,n)
    for (j in 1:d) {
      Xj <- X[, j]
      mean_j <- mean[j]
      var_j <- var[j]
      theta_j <- theta[j]

      outer_sum <- outer(Xj, Xj, "+") / 2 - mean_j
      outer_diff <- outer(Xj, Xj, "-")

      predsig2_component <- predsig2_component * (
        1 / sqrt(1 + 4 * var_j / theta_j) *
          exp(-outer_sum^2 / (theta_j / 2 + 2 * var_j)) * exp(-outer_diff^2 / (2 * theta_j))
      )
    }

    # mean
    predy <- test.mean + sum(predy_component * a)
    # var
    predsig2 <- pmax(0, tau2hat - (predy - test.mean)^2 + sum((drop(a %o% a) - (Ci * tau2hat)) * predsig2_component))

    pred_exact[[1]][1,i] <- predy
    pred_exact[[2]][1,i] <- predsig2
  }
  colnames(pred_exact[[1]]) <- colnames(pred_exact[[2]]) <- colnames(X)

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