dynemu (version 1.0.0)

dynemu_pred: Predictive posterior computation for dynamic systems using one-steap-ahead approach by Gaussian process (GP) emulators

Description

The function computes the predictive posterior distribution (mean and variance) at all time steps for dynamic systems modeled by GP emulators. It can use either the exact computation or Monte Carlo (MC) approximation to compute the predictive posterior.

Usage

dynemu_pred(dynGP, times, xnew, method="Exact", mc.sample=NULL, trace=TRUE)

Value

A list containing:

• times: The time sequence vector.

• mu: A matrix of predictive mean values for each state variable. Each row corresponds to a time step, and each column corresponds to a state variable.

• sig2: A matrix of predictive variance values for each state variable. Each row corresponds to a time step, and each column corresponds to a state variable.

Arguments

dynGP

list of trained GP emulators fitted by dynemu_GP, each corresponding to a state variable.

times

numeric vector specifying the time sequence at which predictions are to be made.

xnew

numeric vector or single-row matrix specifying the initial state values at time 0. The number of columns must match the input dimension of the GP emulators.

method

character specifying the method for prediction, to be chosen between "Exact"(exact closed-form solution), or "MC"(MC approximation). Default is "Exact".

mc.sample

a number of mc samples generated for MC approximation. Required only if method="MC".

trace

logical indicating whether to print the progress bar. If trace=TRUE, the progress bar is printed. Default is TRUE.

Details

Given a trained set of GP emulators `dynGP` fitted using dynemu_GP, this function: 1. Computes the predictive posterior mean and variance for each state variable at all time steps. 2. The function can either: - Use the exact computation for a closed-form solution by method="Exact". - Use the MC approximation method to generate samples and estimate the posterior by method="MC".

Examples

# \donttest{
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 initial test values ###
state <- c(P = 1, C = 2)

### Generate test set ###
test <- dyn_solve(LVmod0D, times, pars, state)

### Define number of MC samples ###
mc.sample <- 1000

### Prediction ###
pred.exact.dyn <- dynemu_pred(fit.dyn, times, state, method="Exact")
pred.exact.mu <- pred.exact.dyn$mu
pred.exact.sig2 <- pred.exact.dyn$sig2

pred.mc.dyn <- dynemu_pred(fit.dyn, times, state, method="MC", trace=TRUE)
pred.mc.mu <- pred.mc.dyn$mu
pred.mc.sig2 <- pred.mc.dyn$sig2

### Compute RMSE ###
sqrt(mean((pred.exact.mu[,1]-test$P)^2)) # 0.03002421
sqrt(mean((pred.exact.mu[,2]-test$C)^2)) # 0.02291742

sqrt(mean((pred.mc.mu[,1]-test$P)^2)) # 0.03050849
sqrt(mean((pred.mc.mu[,2]-test$C)^2)) # 0.02330542
# }