dynemu_pred: Predictive posterior computation for dynamic systems using...

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

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)

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".

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.

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 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", mc.sample=mc.sample, 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.0007686341
sqrt(mean((pred.exact.mu[,2]-test$C)^2)) # 0.0005255164

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

dynemu documentation built on Aug. 19, 2025, 1:11 a.m.