C_GP_ci {activegp}

CI on Eigenvalues via Monte Carlo/GP

Description

CI on Eigenvalues via Monte Carlo/GP

Usage

C_GP_ci(model, B = 100)

Arguments

model

A homGP model

B

Monte Carlo iterates

Value

A list with elements ci giving 95

Examples

################################################################################
## Example of uncertainty quantification on C estimate
################################################################################
library(hetGP); library(lhs)
set.seed(42)

nvar <- 2
n <- 20
nits <- 20

# theta gives the subspace direction
f <- function(x, theta, nugget = 1e-6){
  if(is.null(dim(x))) x <- matrix(x, 1)
  xact <- cos(theta) * x[,1] - sin(theta) * x[,2]
  return(hetGP::f1d(xact) + rnorm(n = nrow(x), sd = rep(nugget, nrow(x))))
}

theta_dir <- pi/6
act_dir <- c(cos(theta_dir), -sin(theta_dir))

# Create design of experiments and initial GP model
design <- X <- matrix(signif(maximinLHS(n, nvar), 2), ncol = nvar)
response <- Y <- apply(design, 1, f, theta = theta_dir)
model <- mleHomGP(design, response, known = list(beta0 =  0))

res <- C_GP_ci(model)

plot(c(1, 2), log(c(mean(res$eigen_draws[,1]), mean(res$eigen_draws[,2]))),
  ylim = range(log(res$eigen_draws)), ylab = "Eigenvalue", xlab = "Index")
  segments(1, log(res$ci[1,1]), 1, log(res$ci[2,1]))
  segments(2, log(res$ci[1,2]), 2, log(res$ci[2,2]))