simul.CRNGP: Fast conditional simulation for a CRNGP model

In hetGP: Heteroskedastic Gaussian Process Modeling and Design under Replication

Fast conditional simulation for a CRNGP model

Description

Fast conditional simulation for a CRNGP model

Usage

## S3 method for class 'CRNGP'
simul(
  object,
  Xgrid,
  ids = NULL,
  nsim = 1,
  eps = sqrt(.Machine$double.eps),
  seqseeds = TRUE,
  check = TRUE
)

Arguments

object	a CRNGP model obtained with mleCRNGP
Xgrid	matrix of (x, seed) locations where the simulation is performed. The last column MUST correspond to seeds values. Xgrid must also contain the evaluated designs (e.g., in object$X0). All design locations are matched with all seed values, either by increasing seed values or repeating the seed sequence.
ids	vector of indices corresponding to observed values in Xgrid
nsim	number of simulations to return
eps	jitter used in the Cholesky decomposition of the covariance matrix for numerical stability
seqseeds	is the seed sequence repeated (e.g., 1 2 3 1 2 3), else it is assumed to be ordered (e.g., 1 1 2 2 3 3)
check	if TRUE, check that Xgrid has the proper structure (slower)

Value

A matrix of size nrow(Xgrid) x nsim.

References

Chiles, J. P., & Delfiner, P. (2012). Geostatistics: modeling spatial uncertainty (Vol. 713). John Wiley & Sons.

Chevalier, C.; Emery, X.; Ginsbourger, D. Fast Update of Conditional Simulation Ensembles Mathematical Geosciences, 2014

Examples

## Not run: 
##------------------------------------------------------------
## Example: Homoskedastic GP modeling on 2d sims
##------------------------------------------------------------
set.seed(2)
nx <- 31
ns <- 5
d <- 2
x <- as.matrix(expand.grid(seq(0,1, length.out = nx), seq(0,1, length.out = nx)))
s <- matrix(seq(1, ns, length.out = ns))
Xgrid <- as.matrix(expand.grid(seq(1, ns, length.out = ns), seq(0,1, length.out = nx), 
                               seq(0,1, length.out = nx)))
Xgrid <- Xgrid[,c(2, 3, 1)]
g <- 1e-6
theta <- c(0.2, 0.5)
KX <- cov_gen(x, theta = theta)
rho <- 0.33
KS <- matrix(rho, ns, ns)
diag(KS) <- 1

YY <- MASS::mvrnorm(n = 1, mu = rep(0, nx*nx*ns), Sigma = kronecker(KX, KS) + g * diag(nx*nx*ns))
YYmat <- matrix(YY, ns, nx*nx)
filled.contour(matrix(YYmat[1,], nx))
filled.contour(matrix(YYmat[2,], nx))

ids <- sample(1:nrow(Xgrid), 80)
X0 <- Xgrid[ids,]
Y0 <-  YY[ids]

# For 3d visualization
# library(rgl)
# plot3d(Xgrid[,1], Xgrid[,2], YY, col = 1 + (Xgrid[,3] - 1) %% 6)
# points3d(X0[,1], X0[,2], Y0, size = 10, col = 1 + ((X0[,3] - 1) %% 6))

model <- mleCRNGP(X0, Y0, known = list(g = 1e-6))

preds <- predict(model, x = Xgrid, xprime = Xgrid)
# surface3d(unique(Xgrid[1:nx^2,1]),unique(Xgrid[,2]), matrix(YY[Xgrid[,3]==1], nx), 
#   front = "lines", back = "lines")
# aspect3d(1, 1, 1)
# surface3d(unique(Xgrid[1:nx^2,1]),unique(Xgrid[,2]), matrix(preds$mean[Xgrid[,3]==1], nx), 
#   front = "lines", back = "lines", col = "red")

# Conditional realizations (classical way)
set.seed(2)
t0 <- Sys.time()
SigmaCond <- 1/2 * (preds$cov + t(preds$cov))
sims <- t(chol(SigmaCond + diag(sqrt(.Machine$double.eps), nrow(Xgrid)))) %*% rnorm(nrow(Xgrid))
sims <- sims + preds$mean
print(difftime(Sys.time(), t0))
# sims <- MASS::mvrnorm(n = 1, mu = preds$mean, Sigma = 1/2 * (preds$cov + t(preds$cov)))
# plot3d(X0[,1], X0[,2], Y0, size = 10, col = 1 + ((X0[,3] - 1) %% 6))
# surface3d(unique(x[,1]), unique(x[,2]), matrix(sims[Xgrid[,3] == 1], nx), col = 1, 
#   front = "lines", back = "lines")
# surface3d(unique(x[,1]), unique(x[,2]), matrix(sims[Xgrid[,3] == 2], nx), col = 2, 
#   front = "lines", back = "lines")

# Alternative for conditional realizations 
# (note: here the design points are part of the simulation points)
set.seed(2)
t0 <- Sys.time()
condreas <- simul(model, Xgrid, ids = ids)
print(difftime(Sys.time(), t0))
# plot3d(X0[,1], X0[,2], Y0, size = 10, col = 1 + ((X0[,3] - 1) %% 6))
# surface3d(unique(x[,1]), unique(x[,2]), matrix(condreas[Xgrid[,3] == 1], nx), col = 1, 
#   front = "lines", back = "lines")
# surface3d(unique(x[,1]), unique(x[,2]), matrix(condreas[Xgrid[,3] == 2], nx), col = 2, 
#   front = "lines", back = "lines")

# Alternative using ordered seeds:
Xgrid2 <- as.matrix(expand.grid(seq(0,1, length.out = nx), 
  seq(0,1, length.out = nx), seq(1, ns, length.out = ns)))
condreas2 <- simul(model, Xgrid2, ids = ids, seqseeds = FALSE)

## Check that values at X0 are coherent:
# condreas[ids,1] - Y0
# sims[ids,1] - Y0

## Check that the empirical mean/covariance is correct
condreas2 <- simul(model, Xgrid, ids = ids, nsim = 1000)
print(range(rowMeans(condreas2) - preds$mean))
print(range(cov(t(condreas2)) - preds$cov))

## End(Not run)

hetGP documentation built on Oct. 3, 2023, 1:07 a.m.