dynemu_pred {dynemu} | R Documentation |
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 |
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 |
mc.sample |
a number of mc samples generated for MC approximation. Required only if |
trace |
logical indicating whether to print the progress bar. If |
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", 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