dynemu_mc {dynemu} | R Documentation |
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.
Usage
dynemu_mc(mean, var, dynGP, mc.sample)
Arguments
mean |
numeric vector or matrix specifying the mean of uncertain inputs. The number of columns must match the input dimension of 'dynGP'. |
var |
numeric vector or matrix specifying the variance of uncertain inputs. The number of columns must match the input dimension of 'dynGP'. |
dynGP |
list of trained GP emulators fitted by |
mc.sample |
a number of mc samples generated for MC approximation. Default is |
Details
Given a trained set of GP emulators 'dynGP' fitted using 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 dynemu_exact
, which provides a closed-form solution,
this function uses MC sampling to approximate the posterior.
Value
A list containing:
-
predy
: A single-row matrix of predictive mean values, with each column for corresponding state variable. -
predsig2
: A single-row matrix of predictive variance values, with each column for corresponding state variable.
References
H. Mohammadi, P. Challenor, and M. Goodfellow (2019). Emulating dynamic non-linear simulators using Gaussian processes. Computational Statistics & Data Analysis, 139, 178-196; doi:10.1016/j.csda.2019.05.006
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)