Robust Gaussian Stochastic Process Emulation

Description

Robust parameter estimation and prediction of Gaussian stochastic process emulators. It allows for robust parameter estimation and prediction using Gaussian stochastic process emulator. It also implements the parallel partial Gaussian stochastic process emulator for computer model with massive outputs See the reference: Mengyang Gu and Jim Berger, 2016, Annals of Applied Statistics; Mengyang Gu, Xiaojing Wang and Jim Berger, 2018, Annals of Statistics.

Details

The DESCRIPTION file:

Index of help topics:

RobustGaSP-package      Robust Gaussian Stochastic Process Emulation
findInertInputs         find inert inputs with the posterior mode
humanity.X              data from the humanity model
leave_one_out_rgasp     leave-one-out fitted values and standard
                        deviation for robust GaSP model
plot                    Plot for Robust GaSP model
ppgasp                  Setting up the parallel partial GaSP model
ppgasp-class            PP GaSP class
predict.ppgasp          Prediction for PP GaSP model
predict.rgasp           Prediction for Robust GaSP model
predppgasp-class        Predicted PP GaSP class
predrgasp-class         Predictive robust GaSP class
rgasp                   Setting up the robust GaSP model
rgasp-class             Robust GaSP class
show                    Show Robust GaSP object
show.ppgasp             Show parllel partial Gaussian stochastic
                        process (PP GaSP) object
simulate                Sample for Robust GaSP model

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

J.O. Berger, V. De Oliveira and B. Sanso (2001), Objective Bayesian analysis of spatially correlated data, Journal of the American Statistical Association, 96, 1361-1374.

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

M. Gu, X. Wang and J.O. Berger (2018), Robust Gaussian stochastic process emulation, Annals of Statistics, 46(6A), 3038-3066.

M. Gu (2018), Jointly robust prior for Gaussian stochastic process in emulation, calibration and variable selection, arXiv:1804.09329.

R. Paulo (2005), Default priors for Gaussian processes, Annals of statistics, 33(2), 556-582.

J. Sacks, W.J. Welch, T.J. Mitchell, and H.P. Wynn (1989), Design and analysis of computer experiments, Statistical Science, 4, 409-435.

Examples

  #------------------------
  # a 3 dimensional example
  #------------------------
  # dimensional of the inputs
  dim_inputs <- 3    
  # number of the inputs
  num_obs <- 30       
  # uniform samples of design
  input <- matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  ##maximin lhd sample
  
  ####
  # outputs from the 3 dim dettepepel.3.data function
  
  output = matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]<-dettepepel.3.data(input[i,])
  }
  
  # use constant mean basis, with no constraint on optimization
  m1<- rgasp(design = input, response = output, lower_bound=FALSE)
  
  # the following use constraints on optimization
  # m1<- rgasp(design = input, response = output, lower_bound=TRUE)
  
  # the following use a single start on optimization
  # m1<- rgasp(design = input, response = output, lower_bound=FALSE)
  
  # number of points to be predicted 
  num_testing_input <- 5000    
  # generate points to be predicted
  testing_input <- matrix(runif(num_testing_input*dim_inputs),num_testing_input,dim_inputs)
  # Perform prediction
  m1.predict<-predict(m1, testing_input, outasS3 = FALSE)
  # Predictive mean
  #m1.predict@mean  
  
  # The following tests how good the prediction is 
  testing_output <- matrix(0,num_testing_input,1)
  for(i in 1:num_testing_input){
    testing_output[i]<-dettepepel.3.data(testing_input[i,])
  }
  
  # compute the MSE, average coverage and average length
  # out of sample MSE
  MSE_emulator <- sum((m1.predict@mean-testing_output)^2)/(num_testing_input)  
  
  # proportion covered by 95% posterior predictive credible interval
  prop_emulator <- length(which((m1.predict@lower95<=testing_output)
                   &(m1.predict@upper95>=testing_output)))/num_testing_input
  
  # average length of  posterior predictive credible interval
  length_emulator <- sum(m1.predict@upper95-m1.predict@lower95)/num_testing_input
  
  # output of prediction
  MSE_emulator
  prop_emulator
  length_emulator  
  # normalized RMSE
  sqrt(MSE_emulator/mean((testing_output-mean(output))^2 ))

Borehole function

Description

It is an 8-dimensional test function that models water flow through a borehole. Its simplicity and quick evaluation makes it a commonly used function for testing a wide variety of methods in computer experiments.

Usage

borehole(x)

Arguments

Details

For more details, see Worley, B. A. (1987). Deterministic uncertainty analysis (No. CONF-871101-30). Oak Ridge National Lab., TN (USA).

Value

The response is given by

f(x) = \frac{2\pi T_u (H_u - H_l)}{\ln(\frac{r}{r_w})\Big(1+\frac{2LT_u}{\ln(\frac{r}{r_w})r^2_w K_w}+\frac{T_u}{T_l}\Big)}

where x is an 8-dimensional vector with rw <- x[1], r <- x[2], Tu <- x[3], Hu <- x[4], Tl <- x[5], Hl <- x[6], L <- x[7], Kw <- x[8].

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Harper, W. V., & Gupta, S. K. (1983). Sensitivity/uncertainty analysis of a borehole scenario comparing Latin Hypercube Sampling and deterministic sensitivity approaches (No. BMI/ONWI-516). Battelle Memorial Inst., Columbus, OH (USA). Office of Nuclear Waste Isolation.

Worley, B. A. (1987). Deterministic uncertainty analysis (No. CONF-871101-30). Oak Ridge National Lab., TN (USA).

S. Surjanovic, D. Bingham, Virtual Library of Simulation Experiments: Test Functions and Datasets, retrieved March 29, 2016, from http://www.sfu.ca/~ssurjano/borehole.html.

Sample for Robust GaSP model

Description

Function to sample Robust GaSP after the Robust GaSP model has been constructed.

Usage

## S4 method for signature 'rgasp'
Sample(object, testing_input, num_sample=1,
testing_trend= matrix(1,dim(testing_input)[1],1),
r0=NA, rr0=NA, sample_data=T,
...)

Arguments

Value

The returned value is a matrix where each column is a sample on the prespecified inputs.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Mengyang Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

Examples

  #------------------------
  # a 1 dimensional example
  #------------------------
  
###########1dim higdon.1.data 
p1 = 1     ###dimensional of the inputs
dim_inputs1 <- p1
n1 = 15   ###sample size or number of training computer runs you have 
num_obs1 <- n1
input1 = 10*matrix(runif(num_obs1*dim_inputs1), num_obs1,dim_inputs1) ##uniform
#####lhs is better
#library(lhs)
#input1 = 10*maximinLHS(n=num_obs1, k=dim_inputs1)  ##maximin lhd sample
output1 = matrix(0,num_obs1,1)
for(i in 1:num_obs1){
  output1[i]=higdon.1.data (input1[i])
}





m1<- rgasp(design = input1, response = output1, lower_bound=FALSE)

#####locations to samples
testing_input1 = seq(0,10,1/50) 
testing_input1=as.matrix(testing_input1)
#####draw 10 samples
m1_sample=Sample(m1,testing_input1,num_sample=10)

#####plot these samples
matplot(testing_input1,m1_sample, type='l',xlab='input',ylab='output')
lines(input1,output1,type='p')

Convert a rgasp or ppgasp S4 object prediction into a S3 object

Description

This function converts the default S4 object prediction into a S3 object

Usage

as.S3prediction(object, ...)

Arguments

Value

The returned value is a list with

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Examples

  #------------------------
  # a 3 dimensional example
  #------------------------
  # dimensional of the inputs
  dim_inputs <- 3    
  # number of the inputs
  num_obs <- 30       
  # uniform samples of design
  input <- matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  ##maximin lhd sample
  
  # outputs from the 3 dim dettepepel.3.data function
  
  output = matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]<-dettepepel.3.data (input[i,])
  }
  
  # use constant mean basis, with no constraint on optimization
  m1<- rgasp(design = input, response = output, lower_bound=FALSE)
  
  # the following use constraints on optimization
  # m1<- rgasp(design = input, response = output, lower_bound=TRUE)
  
  # the following use a single start on optimization
  # m1<- rgasp(design = input, response = output, lower_bound=FALSE, multiple_starts=FALSE)
  
  # number of points to be predicted 
  num_testing_input <- 5000    
  # generate points to be predicted
  testing_input <- matrix(runif(num_testing_input*dim_inputs),num_testing_input,dim_inputs)
  # Perform prediction
  m1.predict<-predict(m1, testing_input, outasS3 = FALSE)
  
  # The returned object is of predrgasp-class
  str(m1.predict)
  # To have the prediction as a list
  m1.predict.aslist <- as.S3prediction(m1.predict)
  str(m1.predict.aslist)
  

Convert a RobustGaSP S3 object prediction into a S4 object

Description

This function converts a S3 object prediction into the default S4 object prediction

Usage

as.S4prediction(object, ...)

Arguments

Value

The returned value is a S4 object of class predrgasp-class with

call:

call to the function.

mean:

predictive mean for the testing inputs.

lower95:

lower bound of the 95% posterior credible interval.

upper95:

upper bound of the 95% posterior credible interval.

sd:

standard deviation of each testing_input.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Examples

  #------------------------
  # a 3 dimensional example
  #------------------------
  # dimensional of the inputs
  dim_inputs <- 3    
  # number of the inputs
  num_obs <- 30       
  # uniform samples of design
  input <- matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  ##maximin lhd sample
  
  # outputs from the 3 dim dettepepel.3.data function
  
  output = matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]<-dettepepel.3.data (input[i,])
  }
  
  # use constant mean basis, with no constraint on optimization
  m1<- rgasp(design = input, response = output, lower_bound=FALSE)
  
  # the following use constraints on optimization
  # m1<- rgasp(design = input, response = output, lower_bound=TRUE)
  
  # the following use a single start on optimization
  # m1<- rgasp(design = input, response = output, lower_bound=FALSE, multiple_starts=FALSE)
  
  # number of points to be predicted 
  num_testing_input <- 5000    
  # generate points to be predicted
  testing_input <- matrix(runif(num_testing_input*dim_inputs),num_testing_input,dim_inputs)
  # Perform prediction
  m1.predict<-predict(m1, testing_input, outasS3 = FALSE)
  # Notice the call slot of the object
  print(m1.predict@call)
  
  # To convert the prediction to a S3 object 
  m1.predict.aslist <- as.S3prediction(m1.predict)
  # To recover back the prediction as a predrgasp-class object
  m1.predict.aspredgasp <- as.S4prediction.predict(m1.predict.aslist)
  str(m1.predict.aslist)
  # Notice that in this case the @call slot is different than the initial
  print(m1.predict.aspredgasp@call)
  

PP GaSP constructor after estimating the parameters

Description

This function constructs the PP GaSP model if the range and noise-variance ratio parameters are given or have been estimated.

Usage

construct_ppgasp(beta, nu, R0, X, zero_mean, output, kernel_type, alpha)

Arguments

Value

A list. The first element is a lower triangular matrix L, a cholesky decomposition of R, i.e. LL^t=R and R is the correlation matrix. The second element is lower triangular matrix LX a cholesky decomposition of X^tR^{-1}X. The third element is a matrix of theta_hat, the mean (trend) parameters. The last element is a vector of sigma2_hat, which is the estimated variance parameter on each function.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Robust GaSP constructor after estimating the parameters

Description

This function constructs the Robust GaSP model if the range and noise-variance ratio parameters are given or have been estimated.

Usage

construct_rgasp(beta, nu, R0, X, zero_mean, output, kernel_type, alpha)

Arguments

Value

A list. The first element is a lower triangular matrix L, a cholesky decomposition of R, i.e. LL^t=R and R is the correlation matrix. The second element is lower triangular matrix LX a cholesky decomposition of X^tR^{-1}X . The third element is a vector theta_hat, the mean (trend) parameters. The last element is numeric valued sigma2_hat, which is the estimated variance parameter.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Dette & Pepelyshev (2010) Curved Function

Description

A 3-dimensional test function.

Usage

dettepepel.3.data(x)

Arguments

Details

For more details see Dette, H., & Pepelyshev, A. (2010). Generalized Latin hypercube design for computer experiments. Technometrics, 52(4).

Value

A real value equal to the function value evaluated at x.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Dette, H., & Pepelyshev, A. (2010). Generalized Latin hypercube design for computer experiments. Technometrics, 52(4).

S. Surjanovic, D. Bingham, Virtual Library of Simulation Experiments: Test Functions and Datasets, retrieved March 29, 2016, from http://www.sfu.ca/~ssurjano/detpep10curv.html.

Environmental model function

Description

This function is the environmental model where the output is the concentration of the pollutant at the space-time grid used in in Bliznyuk et al. (2008). The physical inputs are 4 dimensions.

Usage

environ.4.data(x, s=c(0.5, 1, 1.5, 2, 2.5), t=seq(from=0.3, to=60, by=0.3))

Arguments

Details

Bliznyuk, N., Ruppert, D., Shoemaker, C., Regis, R., Wild, S., & Mugunthan, P. (2008). Bayesian calibration and uncertainty analysis for computationally expensive models using optimization and radial basis function approximation. Journal of Computational and Graphical Statistics, 17(2), 270-294.

Value

A real value of this function evaluated at x.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

N. Bliznyuk, D. Ruppert, C. Shoemaker, R. Regis, S. Wild, & P. Mugunthan (2008). Bayesian calibration and uncertainty analysis for computationally expensive models using optimization and radial basis function approximation. Journal of Computational and Graphical Statistics, 17(2), 270-294.

S. Surjanovic, D. Bingham, Virtual Library of Simulation Experiments: Test Functions and Datasets, retrieved March 29, 2016, from http://www.sfu.ca/~ssurjano/fried.html.

Euclidean distance matrix between two input matrices

Description

Function to construct the euclidean distance matrix with the two input matrices.

Usage

euclidean_distance(input1, input2)

Arguments

Value

The euclidean distance matrix with the number of rows and the number of columns being the number of rows in the first and second input matrices.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Examples

  # dimensional of the inputs
  dim_inputs <- 8    
  
  # number of the inputs
  num_obs <- 30       
  
  # uniform samples of design
  input <-matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 

  # the Euclidean distance matrix 
  R0=euclidean_distance(input, input)

find inert inputs with the posterior mode

Description

The function tests for inert inputs (inputs that barely affect the outputs) using the posterior mode.

Usage

findInertInputs(object,threshold=0.1)

Arguments

Details

This function utilizes the following quantity

object@p*object@beta_hat*object@CL/sum(object@beta_hat*object@CL)

for each input to identify the inert outputs. The average estimated normalized inverse range parameters will be 1. If the estimated normalized inverse range parameters of an input is close to 0, it means this input might be an inert input.

In this method, a prior that has shrinkage effects is suggested to use, .e.g the jointly robust prior (i.e. one should set prior_choice='ref_approx' in rgasp() to obtain the use rgasp object before using this function). Moreover, one may not add a lower bound of the range parameters to perform this method, i.e. one should set lower_bound=F in rgasp(). For more details see Chapter 4 in the reference below.

Mengyang Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

Value

A vector that has the same dimension of the number of inputs indicating how likely the inputs are inerts. The average value is 1. When a value is very close to zero, it tends to be an inert inputs.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Mengyang Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

Examples

  #-----------------------------------------------
  # test for inert inputs in the Borehole function
  #-----------------------------------------------
# dimensional of the inputs
dim_inputs <- 8    
# number of the inputs
num_obs <- 40       

# uniform samples of design
set.seed(0)
input <-matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
# Following codes use maximin Latin Hypercube Design, which is typically better than uniform
# library(lhs)
# input <- maximinLHS(n=num_obs, k=dim_inputs)  # maximin lhd sample

# rescale the design to the domain
input[,1]<-0.05+(0.15-0.05)*input[,1];
input[,2]<-100+(50000-100)*input[,2];
input[,3]<-63070+(115600-63070)*input[,3];
input[,4]<-990+(1110-990)*input[,4];
input[,5]<-63.1+(116-63.1)*input[,5];
input[,6]<-700+(820-700)*input[,6];
input[,7]<-1120+(1680-1120)*input[,7];
input[,8]<-9855+(12045-9855)*input[,8];

# outputs from the 8 dim Borehole function

output=matrix(0,num_obs,1)
for(i in 1:num_obs){
  output[i]=borehole(input[i,])
}





# use constant mean basis with trend, with no constraint on optimization
m3<- rgasp(design = input, response = output, lower_bound=FALSE)

P=findInertInputs(m3)

Friedman Function

Description

A 5-dimensional test function.

Usage

friedman.5.data(x)

Arguments

Details

For more details see Friedman, J. H. (1991). Multivariate adaptive regression splines. The annals of statistics, 19(1), 1-67.

Value

A real value equal to the function value evaluated at x.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Friedman, J. H. (1991). Multivariate adaptive regression splines. The annals of statistics, 19(1), 1-67.

S. Surjanovic, D. Bingham, Virtual Library of Simulation Experiments: Test Functions and Datasets, retrieved March 29, 2016, from http://www.sfu.ca/~ssurjano/fried.html.

A C++ function to generate predictive mean and cholesky decomposition of the scaled covariance function.

Description

The function computes predictive mean and cholesky decomposition of the scaled covariance function.

Usage

generate_predictive_mean_cov(beta, nu, input, X,zero_mean,output,  
testing_input,X_testing, L, LX, theta_hat, sigma2_hat,rr0,r0, 
kernel_type,alpha,method,sample_data)

Arguments

Value

A list of 2 elements. The first is a vector for predictive mean for testing inputs. The second is a scaled covariance matrix of the predictive distribution.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Mengyang Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

Higdon (2002) Function

Description

Higdon (2002) 1-dimensional test function.

Usage

higdon.1.data(s)

Arguments

Details

For more details, see Higdon D. (2002) Space and Space-Time Modeling using Process Convolutions. In: Anderson C.W., Barnett V., Chatwin P.C., El-Shaarawi A.H. (eds) Quantitative Methods for Current Environmental Issues. Springer, London.

Value

A real number equal to the Higdon (2002) function value at s.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Higdon, D. (2002). Space and space-time modeling using process convolutions. In Quantitative methods for current environmental issues (pp. 37-56). Springer London.

S. Surjanovic, D. Bingham, Virtual Library of Simulation Experiments: Test Functions and Datasets, retrieved March 29, 2016, from https://www.sfu.ca/~ssurjano/hig02.html.

Examples

s <- seq(0,10,0.01)
y <- higdon.1.data(s)
plot(s,y, xlab='s',ylab='y',type='l',main='Higdon (2002) function')

data from the humanity model

Description

This data set provides the training data and testing data from the 'diplomatic and military operations in a non-warfighting domain' (DIAMOND) simulator. It porduces the number of casualties during the second day to sixth day after the earthquake and volcanic eruption in Giarre and Catania. See (Overstall and Woods (2016)) for details.

Usage

	data(humanity_model)

Format

Four data frame with observations on the following variables.

humanity.X

A matrix of the training inputs.

humanity.Y

A matrix of the output of the calsualties from the second to sixth day after the the earthquake and volcanic eruption for each set of training inputs.

humanity.Xt

A matrix of the test inputs.

humanity.Yt

A matrix of the test output of the calsualties.

References

A. M. Overstall and D. C. Woods (2016). Multivariate emulation of computer simulators: model selection and diagnostics with application to a humanitarian relief model. Journal of the Royal Statistical Society: Series C (Applied Statistics), 65(4):483-505.

B. Taylor and A. Lane. Development of a novel family of military campaign simulation models. Journal of the Operational Research Society, 55(4):333-339, 2004.

leave-one-out fitted values and standard deviation for robust GaSP model

Description

A function to calculate leave-one-out fitted values and the standard deviation of the prediction on robust GaSP models after the robust GaSP model has been constructed.

Usage

leave_one_out_rgasp(object)

Arguments

Value

A list of 2 elements with

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Mengyang Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

See Also

rgasp

Examples

library(RobustGaSP)
 #------------------------
  # a 3 dimensional example
  #------------------------
  # dimensional of the inputs
  dim_inputs <- 3    
  # number of the inputs
  num_obs <- 30       
  # uniform samples of design
  input <- matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  ##maximin lhd sample
  
  ####
  # outputs from the 3 dim dettepepel.3.data function
  
  output = matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]<-dettepepel.3.data (input[i,])
  }
  
  # use constant mean basis, with no constraint on optimization
  m1<- rgasp(design = input, response = output, lower_bound=FALSE)
  
  ##leave one out predict
  leave_one_out_m1=leave_one_out_rgasp(m1)
  
  ##predictive mean 
  leave_one_out_m1$mean
  ##standard deviation
  leave_one_out_m1$sd
  ##standardized error
  (leave_one_out_m1$mean-output)/leave_one_out_m1$sd

  

Lim et at. (2002) nonpolynomial function

Description

This function is 2-dimensional test function.

Usage

limetal.2.data(x)

Arguments

Details

For more details, see Lim, Y. B., Sacks, J., Studden, W. J., & Welch, W. J. (2002). Design and analysis of computer experiments when the output is highly correlated over the input space. Canadian Journal of Statistics, 30(1), 109-126.

Welch, W. J., Buck, R. J., Sacks, J., Wynn, H. P., Mitchell, T. J., & Morris, M. D. (1992). Screening, predicting, and computer experiments. Technometrics, 34(1), 15-25.

Value

A real value of this function evaluated at x.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Lim, Y. B., Sacks, J., Studden, W. J., & Welch, W. J. (2002). Design and analysis of computer experiments when the output is highly correlated over the input space. , Canadian Journal of Statistics, 30(1), 109-126.

Welch, W. J., Buck, R. J., Sacks, J., Wynn, H. P., Mitchell, T. J., & Morris, M. D. (1992). Screening, predicting, and computer experiments. Technometrics, 34(1), 15-25.

S. Surjanovic, D. Bingham, Virtual Library of Simulation Experiments: Test Functions and Datasets, retrieved March 29, 2016, from http://www.sfu.ca/~ssurjano/limetal02non.html.

The natural logarithm of the jointly robust prior (up to a normalizing constant)

Description

A function to the natural logarithm of the jointly robust prior (up to a normalizing constant).

Usage

log_approx_ref_prior(param, nugget, nugget_est, CL, a, b)

Arguments

Value

The numerical value of the derivative of the approximate reference prior with regard to inverse-range parameters and the nugget-variance ratio parameter. When the nugget is fixed, the derivative is on inverse-range parameters.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

M. Gu (2018), Jointly robust prior for Gaussian stochastic process in emulation, calibration and variable selection, arXiv:1804.09329.

See Also

rgasp

Examples

# inputs
x<-runif(10);
n<-length(x);

# default prior parameters
a<-0.2
b<-n^{-1}*(a+1)
R0<-as.matrix(abs(outer(x,x, "-")))
CL<- mean(R0[which(R0>0)])

# compute the density of log reference prior up to a normalizing constant
param <- seq(-10,10,0.01)
prior <- rep(0,length(param))
for(i in 1:length(param)){
  prior[i] <- exp(log_approx_ref_prior(param[i],nugget=0,nugget_est=FALSE,CL,a,b) )
}
# plot
plot(param,prior,type='l',
                xlab='Logarithm of inverse range parameters',
                ylab='Prior density up to a normalizing constant')

Derivative of the jointly robust prior

Description

The function computes the derivative of the approximate reference prior with regard to inverse range parameter and the nugget-noise ratio parameter (if not fixed). When the nugget is fixed, it only compute the derivative with regard to the inverse range parameter; otherwise it produces derivative with regard to inverse range parameter and noise ratio parameter.

Usage

log_approx_ref_prior_deriv(param, nugget, nugget_est, CL, a, b)

Arguments

Value

The numerical value of the derivative of the jointly robust prior with regard to beta (the inverse-range parameters) and nugget (the nugget-variance ratio parameter). When the nugget is fixed, the derivative is on inverse-range parameters.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

See Also

log_approx_ref_prior,rgasp

Natural logarithm of marginal likelihood of the robust GaSP model

Description

This function computes the natural logarithm of marginal likelihood after marginalizing out the mean (trend) and variance parameters by the location-scale prior.

Usage

log_marginal_lik(param, nugget, nugget_est, R0, X, zero_mean,output, kernel_type, alpha)

Arguments

Value

The numerical value of natural logarithm of the marginal likelihood.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

Derivative of natural logarithm of the marginal likelihood

Description

The derivative of natural logarithm of marginal likelihood of the Robust GaSP model with regard to inverse range parameters and nugget-variance ratio parameter after marginalizing out the mean (trend) and variance parameters the location-scale prior. When the nugget is fixed, it only computes the derivative with regard to the inverse range parameter; otherwise it produces derivative with regard to inverse range parameter and nugget-variance ratio parameter.

Usage

log_marginal_lik_deriv(param, nugget, nugget_est, R0, X, zero_mean, 
output, kernel_type, alpha)

Arguments

Value

The numerical value of the derivative of natural logarithm of marginal likelihood with regard to range and nugget-variance ratio parameter (if not fixed). When the nugget is fixed, the derivative is on inverse-range parameters.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

See Also

log_marginal_lik.

Derivative of natural logarithm of the marginal likelihood

Description

The derivative of natural logarithm of marginal likelihood of the PP GaSP model with regard to inverse range parameters and nugget-variance ratio parameter after marginalizing out the mean (trend) and variance parameters by the location-scale prior. When the nugget is fixed, it only compute the derivative with regard to the inverse range parameter; otherwise it produces derivative with regard to inverse range parameter and nugget-variance ratio parameter.

Usage

log_marginal_lik_deriv_ppgasp(param, nugget, nugget_est, R0, X, zero_mean, 
output, kernel_type, alpha)

Arguments

Value

The numerical value of the derivative of natural logarithm of marginal likelihood with regard to range and nugget-variance ratio parameter (if not fixed). When the nugget is fixed, the derivative is on inverse-range parameters.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

See Also

log_marginal_lik_ppgasp.

Natural logarithm of marginal likelihood of the PP GaSP model

Description

This function computes the natural logarithm of marginal likelihood of the PP GaSP model after marginalizing out the mean (trend) and variance parameters by the location-scale prior.

Usage

log_marginal_lik_ppgasp(param, nugget, nugget_est, R0, X, zero_mean,output,
kernel_type, alpha)

Arguments

Value

The numerical value of natural logarithm of the marginal likelihood.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

Natural logarithm of profile likelihood of the robust GaSP model

Description

This function computes the natural logarithm of profile likelihood after plugging in the maximum likelihood estimator of the mean (trend) and variance parameters.

Usage

log_profile_lik(param, nugget, nugget_est, R0, X, zero_mean,output, kernel_type, alpha)

Arguments

Value

The numerical value of natural logarithm of the profile likelihood.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

Derivative of natural logarithm of the profile likelihood

Description

The derivative of natural logarithm of profile likelihood of the Robust GaSP model with regard to inverse range parameters and nugget-variance ratio parameter after plugging in the maximum likelihood estimator of the mean (trend) and variance parameters. When the nugget parameter is fixed, it only computes the derivative with regard to the inverse range parameter; otherwise it produces derivative with regard to inverse range parameter and nugget-variance ratio parameter.

Usage

log_profile_lik_deriv(param, nugget, nugget_est, R0, X, zero_mean, 
output, kernel_type, alpha)

Arguments

Value

The numerical value of the derivative of natural logarithm of marginal likelihood with regard to range and nugget-variance ratio parameter (if not fixed). When the nugget is fixed, the derivative is on inverse-range parameters.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

See Also

log_profile_lik.

Derivative of natural logarithm of the profile likelihood

Description

The derivative of natural logarithm of profile likelihood of the PP GaSP model with regard to inverse range parameters and nugget-variance ratio parameter after plugging in the maximum likelihood estimator of the mean (trend) and variance parameters. When the nugget is fixed, it only compute the derivative with regard to the inverse range parameter; otherwise it produces derivative with regard to inverse range parameter and nugget-variance ratio parameter.

Usage

log_profile_lik_deriv_ppgasp(param, nugget, nugget_est, R0, X, zero_mean, 
output, kernel_type, alpha)

Arguments

Value

The numerical value of the derivative of natural logarithm of profile likelihood with regard to range and nugget-variance ratio parameter (if not fixed). When the nugget is fixed, the derivative is on inverse-range parameters.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

See Also

log_profile_lik_ppgasp.

Natural logarithm of profile likelihood of the PP GaSP model

Description

This function computes the natural logarithm of profile likelihood of the PP GaSP model after plugging in the maximum likelihood estimator of the mean (trend) and variance parameters.

Usage

log_profile_lik_ppgasp(param, nugget, nugget_est, R0, X, zero_mean
,output,kernel_type, alpha)

Arguments

Value

The numerical value of natural logarithm of the profile likelihood.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

Natural logarithm of reference marginal posterior density of the robust GaSP model

Description

This function computes the natural logarithm of marginal posterior density with reference prior of inverse range parameter (beta parameterization) after marginalizing out the mean (trend) and variance parameters by the location-scale prior.

Usage

log_ref_marginal_post(param, nugget, nugget_est, R0, X, zero_mean,
output, kernel_type, alpha)

Arguments

Value

The natural logarithm of marginal posterior density with reference prior of inverse range parameter (beta parameterization).

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Mengyang Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

Natural logarithm of reference marginal posterior density of the PP GaSP model

Description

This function computes the natural logarithm of marginal posterior density of the PP GaSP model with reference prior of inverse range parameter (beta parameterization) after marginalizing out the mean (trend) and variance parameters by the location-scale prior.

Usage

log_ref_marginal_post_ppgasp(param, nugget, nugget_est, R0, X, zero_mean,
output, kernel_type, alpha)

Arguments

Value

The natural logarithm of marginal posterior density with reference prior of inverse range parameter (beta parameterization).

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

The derivative of matern correlation function (roughness parameter equal to 3/2) with regard to the inverse range parameter

Description

This function computes the derivative of a correlation matrix using 1-dimensional matern correlation function (roughness parameter equal to 3/2), with regard to the inverse range parameter.

Usage

matern_3_2_deriv(R0_i, R, beta_i)

Arguments

Value

A matrix in which each element is the derivative of 1-dimensional matern correlation function with regard to the inverse range parameter.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

See Also

matern_3_2_funct.

Matern correlation function with roughness parameter equal to 3/2.

Description

This function computes the values of Matern correlation function with roughness parameter equal to 3/2.

Usage

matern_3_2_funct(d, beta_i)

Arguments

Value

A matrix in which each element is the value of the Matern correlation function with roughness parameter equal to 1.5 evaluated at that location.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

The derivative of matern correlation function with regard to the inverse range parameter

Description

This function computes the derivative of a correlation matrix using 1-dimensional matern correlation function, with regard to the inverse range parameter.

Usage

matern_5_2_deriv(R0_i, R, beta_i)

Arguments

Value

A matrix in which each element is the derivative of 1-dimensional matern correlation function with regard to the inverse range parameter.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

See Also

matern_5_2_funct.

Matern correlation function with roughness parameter equal to 5/2.

Description

This function computes the values of Matern correlation function with roughness parameter equal to 5/2.

Usage

matern_5_2_funct(d, beta_i)

Arguments

Value

A matrix in which each element is the value of the Matern correlation function with roughness parameter equal to 2.5 evaluated at that location.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Natural logarithm of approximate reference marginal posterior density of the robust GaSP model

Description

Natural logarithm of marginal posterior density with the approximate reference prior of inverse range parameter after marginalizing out the mean (trend) and variance parameters by the location-scale prior.

Usage

neg_log_marginal_post_approx_ref(param, nugget, nugget.est
,R0, X, zero_mean,output, CL, a, b,kernel_type, alpha)

Arguments

Value

The natural logarithm of the marginal posterior density with approximate reference prior of inverse range parameter (beta parameterization).

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Mengyang Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

Derivative of negative natural logarithm of approximate reference marginal posterior density

Description

The function computes the derivative (with regard to log of inverse range parameter) of natural logarithm of marginal posterior density with the jointly robust prior prior after marginalizing out the mean (trend) and variance parameters by the location-scale prior.

Usage

neg_log_marginal_post_approx_ref_deriv(param, nugget, nugget.est, 
                                       R0, X, zero_mean,output, CL, a, b, 
                                      kernel_type, alpha)

Arguments

Value

The derivative of natural logarithm of marginal posterior density with jointly robust prior prior.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Mengyang Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

M. Gu (2018), Jointly Robust Prior for Gaussian Stochastic Process in Emulation, Calibration and Variable Selection, arXiv:1804.09329.

See Also

neg_log_marginal_post_approx_ref.

Derivative of negative natural logarithm of approximate reference marginal posterior density of the PP GaSP model

Description

The function computes the derivative (with regard to log of inverse range parameter) of natural logarithm of marginal posterior density of the PP GaSP model with the jointly robust prior after marginalizing out the mean (trend) and variance parameters by the location-scale prior.

Usage

neg_log_marginal_post_approx_ref_deriv_ppgasp(param, nugget, nugget.est, 
                                       R0, X, zero_mean,output, CL, a, b, 
                                      kernel_type, alpha)

Arguments

Value

The derivative of natural logarithm of marginal posterior density with the jointly robust prior.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

M. Gu (2018), Jointly robust prior for Gaussian stochastic process in emulation, calibration and variable selection, arXiv:1804.09329.

See Also

neg_log_marginal_post_approx_ref_ppgasp.

Natural logarithm of approximate reference marginal posterior density of the PP GaSP model

Description

Natural logarithm of marginal posterior density with the jointly robust prior of inverse range parameter of the PP GaSP model after marginalizing out the mean (trend) and variance parameters by the location-scale prior.

Usage

neg_log_marginal_post_approx_ref_ppgasp(param, nugget, 
nugget.est, R0, X, zero_mean,output, CL, a, b,kernel_type, alpha)

Arguments

Value

The natural logarithm of the marginal posterior density with the jointly robust prior prior of inverse range parameter (beta parameterization).

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

M. Gu (2018), Jointly Robust Prior for Gaussian Stochastic Process in Emulation, Calibration and Variable Selection, arXiv:1804.09329.

Negative natural logarithm of reference marginal posterior density of the robust GaSP model with regard to a specific parameterization.

Description

Negative natural logarithm of marginal posterior density (with regard to a specific parameterization) with reference prior of inverse range parameter (beta parameterization) after marginalizing out the mean (trend) and variance parameters by the location-scale prior.

Usage

neg_log_marginal_post_ref(param, nugget, nugget.est, 
                          R0, X, zero_mean,output, prior_choice, 
                          kernel_type, alpha)

Arguments

Value

The negative natural logarithm of marginal posterior density with reference prior with regard to a specific parameterization.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Mengyang Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

Negative natural logarithm of reference marginal posterior density of the PP GaSP model with regard to a specific parameterization.

Description

Negative natural logarithm of marginal posterior density (with regard to a specific parameterization) of the PP GaSP model with reference prior of inverse range parameter (beta parameterization) after marginalizing out the mean (trend) and variance parameters by the location-scale prior.

Usage

neg_log_marginal_post_ref_ppgasp(param, nugget, nugget.est
,R0, X, zero_mean,output, prior_choice,kernel_type, alpha)

Arguments

Value

The negative natural logarithm of marginal posterior density with reference prior with regard to a specific parameterization.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

M. Gu (2018), Jointly Robust Prior for Gaussian Stochastic Process in Emulation, Calibration and Variable Selection, arXiv:1804.09329.

The derivative of periodic exponential correlation function with regard to the inverse range parameter

Description

The function computes the derivative of a correlation matrix parameterized by the periodic exponential correlation function, with regard to the inverse range parameter.

Usage

periodic_exp_deriv(R0_i, R, beta_i)

Arguments

the inverse range parameter.

Value

A matrix in which each element is the derivative of periodic exponential correlation function with regard to the inverse range parameter.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

See Also

periodic_exp_funct.

Periodic folding of the exponential correlation function .

Description

This function computes the values of periodic folding of the exponential correlation function (i.e. the power expoential correlation with roughness parameter being 2).

Usage

periodic_exp_funct(d, beta_i)

Arguments

Value

A matrix in which each element is the value of the periodic folding of the exponential correlation function evaluated at that location.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

The derivative of periodic Gaussian correlation function with regard to the inverse range parameter

Description

The function computes the derivative of a correlation matrix parameterized by the periodic Gaussian correlation function, with regard to the inverse range parameter.

Usage

periodic_gauss_deriv(R0_i, R, beta_i)

Arguments

the inverse range parameter.

Value

A matrix in which each element is the derivative of periodic Gaussian correlation function with regard to the inverse range parameter.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

See Also

periodic_gauss_funct.

Periodic folding of the Gaussian correlation function .

Description

This function computes the values of periodic folding of the Gaussian correlation function (i.e. the power expoential correlation with roughness parameter being 2).

Usage

periodic_gauss_funct(d, beta_i)

Arguments

Value

A matrix in which each element is the value of the periodic folding of the Gaussian correlation function evaluated at that location.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Plot for Robust GaSP model

Description

Function to make plots on Robust GaSP models after the Robust GaSP model has been constructed.

Usage

## S4 method for signature 'rgasp'
plot(x,y, ...)

Arguments

Value

Three plots: the leave-one-out fitted values versus exact values, standardized residuals and QQ plot.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

Examples

 library(RobustGaSP)
  #------------------------
  # a 3 dimensional example
  #------------------------
  # dimensional of the inputs
  dim_inputs <- 3    
  # number of the inputs
  num_obs <- 30       
  # uniform samples of design
  input <- matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  ##maximin lhd sample
  
  # outputs from the 3 dim dettepepel.3.data function
  
  output = matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]<-dettepepel.3.data (input[i,])
  }
  
  # use constant mean basis, with no constraint on optimization
  m1<- rgasp(design = input, response = output, lower_bound=FALSE)
 
  # plot
  plot(m1)

The derivative of power exponential correlation function with regard to the inverse range parameter

Description

The function computes the derivative of a correlation matrix using 1-dimensional power exponential correlation function, with regard to the inverse range parameter.

Usage

pow_exp_deriv(R0_i, R, beta_i, alpha_i)

Arguments

the inverse range parameter.

Value

A matrix in which each element is the derivative of 1-dimensional power exponential correlation function with regard to the inverse range parameter.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

See Also

pow_exp_funct.

Pow exponential correlation function.

Description

This function computes the values of power exponential correlation function.

Usage

pow_exp_funct(d, beta_i, alpha_i)

Arguments

Value

A matrix in which each element is the value of the power exponential function evaluated at that location.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Setting up the parallel partial GaSP model

Description

Setting up the parallel partial GaSP model for estimating the parameters (if the parameters are not given).

Usage

   ppgasp(design,response,trend=matrix(1,dim(response)[1],1),zero.mean="No",nugget=0,
    nugget.est=F,range.par=NA,method='post_mode',prior_choice='ref_approx',a=0.2,
    b=1/(length(response))^{1/dim(as.matrix(design))[2]}*(a+dim(as.matrix(design))[2]),
    kernel_type='matern_5_2',isotropic=F,R0=NA,
    optimization='lbfgs',
    alpha=rep(1.9,dim(as.matrix(design))[2]),lower_bound=T,
    max_eval=max(30,20+5*dim(design)[2]),initial_values=NA,num_initial_values=2)
 

Arguments

Value

ppgasp returns a S4 object of class ppgasp (see ppgasp-class).

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu, X. Wang and J.O. Berger (2018), Robust Gaussian stochastic process emulation, Annals of Statistics, 46(6A), 3038-3066.

M. Gu (2018), Jointly robust prior for Gaussian stochastic process in emulation, calibration and variable selection, arXiv:1804.09329.

J. Nocedal (1980), Updating quasi-Newton matrices with limited storage, Math. Comput., 35, 773-782.

D. C. Liu and J. Nocedal (1989), On the limited memory BFGS method for large scale optimization, Math. Programming, 45, p. 503-528.

Brent, R. (1973), Algorithms for Minimization without Derivatives. Englewood Cliffs N.J.: Prentice-Hall.

Examples

  library(RobustGaSP)
  
  ###parallel partial Gaussian stochastic process (PP GaSP) model 
  ##for the humanity model
  data(humanity_model)
  ##120 runs. The input has 13 variables and output is 5 dimensional.
  ##PP GaSP Emulator
  m.ppgasp=ppgasp(design=humanity.X,response=humanity.Y,nugget.est= TRUE)
  show(m.ppgasp)

  ##make predictions
  m_pred=predict(m.ppgasp,humanity.Xt)
  sqrt(mean((m_pred$mean-humanity.Yt)^2))
  mean(m_pred$upper95>humanity.Yt & humanity.Yt>m_pred$lower95)
  mean(m_pred$upper95-m_pred$lower95)
  sqrt( mean( (mean(humanity.Y)-humanity.Yt)^2 ))

  ##with a linear trend on the selected input performs better
  ## Not run: 
    ###PP GaSP Emulation with a linear trend for the humanity model
    data(humanity_model)
    ##pp gasp with trend
    n<-dim(humanity.Y)[1]
    n_testing=dim(humanity.Yt)[1]
    H=cbind(matrix(1,n,1),humanity.X$foodC)
    H_testing=cbind(matrix(1,n_testing,1),humanity.Xt$foodC)
    m.ppgasp_trend=ppgasp(design=humanity.X,response=humanity.Y,trend=H, 
    nugget.est= TRUE)
    
    show(m.ppgasp_trend)
    
    ##make predictions
    m_pred_trend=predict(m.ppgasp_trend,humanity.Xt,testing_trend=H_testing)
    sqrt(mean((m_pred_trend$mean-humanity.Yt)^2))
    mean(m_pred_trend$upper95>humanity.Yt & humanity.Yt>m_pred_trend$lower95)
    mean(m_pred_trend$upper95-m_pred_trend$lower95)
  
## End(Not run)

PP GaSP class

Description

S4 class for PP GaSP model if the range and noise-variance ratio parameters are given and/or have been estimated.

Objects from the Class

Objects of this class are created and initialized with the function ppgasp that computes the calculations needed for setting up the analysis.

Slots

p:

Object of class integer. The dimensions of the inputs.

num_obs:

Object of class integer. The number of observations.

k:

Object of class integer. The number of outputs in each computer model run.

input:

Object of class matrix with dimension n x p. The design of experiments.

output:

Object of class matrix with dimension n x k. Each row denotes a output vector in each run of the computer model.

X:

Object of class matrix of with dimension n x q. The mean basis function, i.e. the trend function.

zero_mean:

A character to specify whether the mean is zero or not. "Yes" means it has zero mean and "No"" means the mean is not zero.

q:

Object of class integer. The number of mean basis.

LB:

Object of class vector with dimension p x 1. The lower bound for inverse range parameters beta.

beta_initial:

Object of class vector with the initial values of inverse range parameters p x 1.

beta_hat:

Object of class vector with dimension p x 1. The inverse-range parameters.

log_post:

Object of class numeric with the logarithm of marginal posterior.

R0:

Object of class list of matrices where the j-th matrix is an absolute difference matrix of the j-th input vector.

theta_hat:

Object of class vector with dimension q x 1. The the mean (trend) parameter.

L:

Object of class matrix with dimension n x n. The Cholesky decomposition of the correlation matrix R, i.e.

L\%*\%t(L)=R

sigma2_hat:

Object of the class matrix. The estimated variance parameter of each output.

LX:

Object of the class matrix with dimension q x q. The Cholesky decomposition of the correlation matrix

t(X)\%*\%R^{-1}\%*\%X

CL:

Object of the class vector used for the lower bound and the prior.

nugget:

A numeric object used for the noise-variance ratio parameter.

nugget.est:

A logical object of whether the nugget is estimated (T) or fixed (F).

kernel_type:

A vector of character to specify the type of kernel to use.

alpha:

Object of class vector with dimension p x 1 for the roughness parameters in the kernel.

method:

Object of class character to specify the method of parameter estimation. There are three values: post_mode, mle and mmle.

isotropic:

Object of class logical to specify whether the kernel is isotropic.

call:

The call to ppgasp function to create the object.

Methods

show

Prints the main slots of the object.

predict

See predict.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

See Also

RobustGaSP for more details about how to create a RobustGaSP object.

Prediction for robust GaSP model

Description

A function to make prediction on robust GaSP models after the robust GaSP model has been constructed.

Usage

pred_rgasp(beta, nu, input, X, zero_mean,output, testing_input,
           X_testing, L, LX, theta_hat, sigma2_hat, 
           q_025, q_975, r0, kernel_type, alpha,method,interval_data)

Arguments

Value

A list of 4 elements. The first is a vector for predictive mean for testing inputs. The second is a vector for lower quantile for 95% posterior credible interval and the third is the upper quantile for 95% posterior credible interval for these testing inputs. The last is a vector of standard deviation of each testing inputs.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Mengyang Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

See Also

predict.rgasp

Prediction for PP GaSP model

Description

Function to make prediction on the PP GaSP model after the PP GaSP model has been constructed.

Usage

## S4 method for signature 'ppgasp'
predict(object, testing_input, 
testing_trend= matrix(1,dim(testing_input)[1],1),r0=NA, 
interval_data=T,
outasS3 = T,loc_index=NA, ...)

Arguments

Value

If the parameter outasS3=F, then the returned value is a S4 object of class predppgasp-class with

If the parameter outasS3=T, then the returned value is a list with

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

Examples

  library(RobustGaSP)
  #----------------------------------
  # an example of environmental model
  #----------------------------------
  
  set.seed(1)
  #Here the sample size is very small. Consider to use more observations 
  n=80
  p=4
  ##using the latin hypercube will be better
  #library(lhs)
  #input_samples=maximinLHS(n,p)
  input_samples=matrix(runif(n*p),n,p)
  input=matrix(0,n,p)
  input[,1]=7+input_samples[,1]*6
  input[,2]=0.02+input_samples[,2]*1
  input[,3]=0.01+input_samples[,3]*2.99
  input[,4]=30.01+input_samples[,4]*0.285
  
  k=300
  output=matrix(0,n,k)
  ##environ.4.data is an environmental model on a spatial-time vector
  ##? environ.4.data
  for(i in 1:n){
    output[i,]=environ.4.data(input[i,],s=seq(0.15,3,0.15),t=seq(4,60,4)  )
  }
  
  ##samples some test inputs
  n_star=1000
  sample_unif=matrix(runif(n_star*p),n_star,p)
  
  testing_input=matrix(0,n_star,p)
  testing_input[,1]=7+sample_unif[,1]*6
  testing_input[,2]=0.02+sample_unif[,2]*1
  testing_input[,3]=0.01+sample_unif[,3]*2.99
  testing_input[,4]=30.01+sample_unif[,4]*0.285
  
  
  testing_output=matrix(0,n_star,k)
  
  s=seq(0.15,3,0.15)
  t=seq(4,60,4) 
  
  for(i in 1:n_star){
    testing_output[i,]=environ.4.data(testing_input[i,],s=s,t=t )
  }
  
  ##we do a transformation of the output 
  ##one can change the number of initial values to test
  log_output_1=log(output+1)
  #since we have lots of output, we use 'nelder-mead' for optimization
  m.ppgasp=ppgasp(design=input,response=log_output_1,kernel_type
                  ='pow_exp',num_initial_values=2,optimization='nelder-mead')
  
  m_pred.ppgasp=predict(m.ppgasp,testing_input)
  ##we transform back for the prediction
  m_pred_ppgasp_median=exp(m_pred.ppgasp$mean)-1
  ##mean squared error
  mean( (m_pred_ppgasp_median-testing_output)^2)
  ##variance of the testing outputs
  var(as.numeric(testing_output))
  
  ##makes plots for the testing 
  par(mfrow=c(1,2))
  testing_plot_1=matrix(testing_output[1,],  length(t), length(s) )
  
  max_testing_plot_1=max(testing_plot_1)
  min_testing_plot_1=min(testing_plot_1)
  
  image(x=t,y=s,testing_plot_1,  col = hcl.colors(100, "terrain"),main='test outputs')
  contour(x=t,y=s,testing_plot_1, levels = seq(min_testing_plot_1, max_testing_plot_1,
                                               by = (max_testing_plot_1-min_testing_plot_1)/5),
          add = TRUE, col = "brown")
  
  ppgasp_plot_1=matrix(m_pred_ppgasp_median[1,],  length(t), length(s) )
  max_ppgasp_plot_1=max(ppgasp_plot_1)
  min_ppgasp_plot_1=min(ppgasp_plot_1)
  
  image(x=t,y=s,ppgasp_plot_1,  col = hcl.colors(100, "terrain"),main='prediction')
  contour(x=t,y=s,ppgasp_plot_1, levels = seq(min_testing_plot_1, max_ppgasp_plot_1,
                                              by = (max_ppgasp_plot_1-min_ppgasp_plot_1)/5),
          add = TRUE, col = "brown")
  dev.off()
  

Prediction for Robust GaSP model

Description

Function to make prediction on the robust GaSP model after the robust GaSP model has been constructed.

Usage

## S4 method for signature 'rgasp'
predict(object,testing_input,testing_trend= matrix(1,dim(testing_input)[1],1),
r0=NA,interval_data=T,
outasS3 = T,...)

Arguments

Value

If the parameter outasS3=F, then the returned value is a S4 object of class predrgasp-class with

call:

call to predict.rgasp function where the returned object has been created.

mean:

predictive mean for the testing inputs.

lower95:

lower bound of the 95% posterior credible interval.

upper95:

upper bound of the 95% posterior credible interval.

sd:

standard deviation of each testing_input.

If the parameter outasS3=T, then the returned value is a list with

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

M. Gu, X. Wang and J.O. Berger (2018), Robust Gaussian Stochastic Process Emulation, Annals of Statistics, 46(6A), 3038-3066.

M. Gu (2018), Jointly Robust Prior for Gaussian Stochastic Process in Emulation, Calibration and Variable Selection, arXiv:1804.09329.

Examples

  #------------------------
  # a 3 dimensional example
  #------------------------
  # dimensional of the inputs
  dim_inputs <- 3    
  # number of the inputs
  num_obs <- 30       
  # uniform samples of design
  input <- matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  ##maximin lhd sample
  
  # outputs from the 3 dim dettepepel.3.data function
  
  output = matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]<-dettepepel.3.data (input[i,])
  }
  
  # use constant mean basis, with no constraint on optimization
  m1<- rgasp(design = input, response = output, lower_bound=FALSE)
  
  # the following use constraints on optimization
  # m1<- rgasp(design = input, response = output, lower_bound=TRUE)
  
  # the following use a single start on optimization
  # m1<- rgasp(design = input, response = output, lower_bound=FALS)
  
  # number of points to be predicted 
  num_testing_input <- 5000    
  # generate points to be predicted
  testing_input <- matrix(runif(num_testing_input*dim_inputs),num_testing_input,dim_inputs)
  # Perform prediction
  m1.predict<-predict(m1, testing_input)
  # Predictive mean
  # m1.predict$mean  
  
  # The following tests how good the prediction is 
  testing_output <- matrix(0,num_testing_input,1)
  for(i in 1:num_testing_input){
    testing_output[i]<-dettepepel.3.data(testing_input[i,])
  }
  
  # compute the MSE, average coverage and average length
  # out of sample MSE
  MSE_emulator <- sum((m1.predict$mean-testing_output)^2)/(num_testing_input)  
  
  # proportion covered by 95% posterior predictive credible interval
  prop_emulator <- length(which((m1.predict$lower95<=testing_output)
                   &(m1.predict$upper95>=testing_output)))/num_testing_input
  
  # average length of  posterior predictive credible interval
  length_emulator <- sum(m1.predict$upper95-m1.predict$lower95)/num_testing_input
  
  # output of prediction
  MSE_emulator
  prop_emulator
  length_emulator  
  # normalized RMSE
  sqrt(MSE_emulator/mean((testing_output-mean(output))^2 ))


  #-----------------------------------
  # a 2 dimensional example with trend
  #-----------------------------------
  # dimensional of the inputs
  dim_inputs <- 2    
  # number of the inputs
  num_obs <- 20       
  
  # uniform samples of design
  input <-matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  ##maximin lhd sample
  
  # outputs from the 2 dim Brainin function
  
  output <- matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]<-limetal.2.data (input[i,])
  }
  
  #mean basis (trend)
  X<-cbind(rep(1,num_obs), input )
  
  
  # use constant mean basis with trend, with no constraint on optimization
  m2<- rgasp(design = input, response = output,trend =X,  lower_bound=FALSE)
  
  
  # number of points to be predicted 
  num_testing_input <- 5000    
  # generate points to be predicted
  testing_input <- matrix(runif(num_testing_input*dim_inputs),num_testing_input,dim_inputs)
  
  # trend of testing
  testing_X<-cbind(rep(1,num_testing_input), testing_input )
  
  
  # Perform prediction
  m2.predict<-predict(m2, testing_input,testing_trend=testing_X)
  # Predictive mean
  #m2.predict$mean  
  
  # The following tests how good the prediction is 
  testing_output <- matrix(0,num_testing_input,1)
  for(i in 1:num_testing_input){
    testing_output[i]<-limetal.2.data(testing_input[i,])
  }
  
  # compute the MSE, average coverage and average length
  # out of sample MSE
  MSE_emulator <- sum((m2.predict$mean-testing_output)^2)/(num_testing_input)  
  
  # proportion covered by 95% posterior predictive credible interval
  prop_emulator <- length(which((m2.predict$lower95<=testing_output)
                   &(m2.predict$upper95>=testing_output)))/num_testing_input
  
  # average length of  posterior predictive credible interval
  length_emulator <- sum(m2.predict$upper95-m2.predict$lower95)/num_testing_input
  
  # output of prediction
  MSE_emulator
  prop_emulator
  length_emulator  
  # normalized RMSE
  sqrt(MSE_emulator/mean((testing_output-mean(output))^2 ))


    ###here try the isotropic kernel (a function of Euclidean distance)
  m2_isotropic<- rgasp(design = input, response = output,trend =X,  
             lower_bound=FALSE,isotropic=TRUE)
  
  m2_isotropic.predict<-predict(m2_isotropic, testing_input,testing_trend=testing_X)
  
  # compute the MSE, average coverage and average length
  # out of sample MSE
  MSE_emulator_isotropic <- sum((m2_isotropic.predict$mean-testing_output)^2)/(num_testing_input)
  
  # proportion covered by 95% posterior predictive credible interval
  prop_emulator_isotropic <- length(which((m2_isotropic.predict$lower95<=testing_output)
                                &(m2_isotropic.predict$upper95>=testing_output)))/num_testing_input
  
  # average length of  posterior predictive credible interval
  length_emulator_isotropic <- sum(m2_isotropic.predict$upper95-
  m2_isotropic.predict$lower95)/num_testing_input
  
  MSE_emulator_isotropic
  prop_emulator_isotropic
  length_emulator_isotropic
  ##the result of isotropic kernel is not as good as the product kernel for this example


  #--------------------------------------------------------------------------------------
  # an 8 dimensional example using only a subset inputs and a noise with unknown variance
  #--------------------------------------------------------------------------------------
  set.seed(1)
  # dimensional of the inputs
  dim_inputs <- 8    
  # number of the inputs
  num_obs <- 50       
  
  # uniform samples of design
  input <-matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  # maximin lhd sample
  
  # rescale the design to the domain
  input[,1]<-0.05+(0.15-0.05)*input[,1];
  input[,2]<-100+(50000-100)*input[,2];
  input[,3]<-63070+(115600-63070)*input[,3];
  input[,4]<-990+(1110-990)*input[,4];
  input[,5]<-63.1+(116-63.1)*input[,5];
  input[,6]<-700+(820-700)*input[,6];
  input[,7]<-1120+(1680-1120)*input[,7];
  input[,8]<-9855+(12045-9855)*input[,8];
  
  # outputs from the 8 dim Borehole function
  
  output=matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]=borehole(input[i,])
  }
  
  
    
    
  
  # use constant mean basis with trend, with no constraint on optimization
  m3<- rgasp(design = input[,c(1,4,6,7,8)], response = output,
             nugget.est=TRUE, lower_bound=FALSE)
  
  
  # number of points to be predicted 
  num_testing_input <- 5000    
  # generate points to be predicted
  testing_input <- matrix(runif(num_testing_input*dim_inputs),num_testing_input,dim_inputs)
  
  # resale the points to the region to be predict
  testing_input[,1]<-0.05+(0.15-0.05)*testing_input[,1];
  testing_input[,2]<-100+(50000-100)*testing_input[,2];
  testing_input[,3]<-63070+(115600-63070)*testing_input[,3];
  testing_input[,4]<-990+(1110-990)*testing_input[,4];
  testing_input[,5]<-63.1+(116-63.1)*testing_input[,5];
  testing_input[,6]<-700+(820-700)*testing_input[,6];
  testing_input[,7]<-1120+(1680-1120)*testing_input[,7];
  testing_input[,8]<-9855+(12045-9855)*testing_input[,8];
  
  
  # Perform prediction
  m3.predict<-predict(m3, testing_input[,c(1,4,6,7,8)])
  # Predictive mean
  #m3.predict$mean  
  
  # The following tests how good the prediction is 
  testing_output <- matrix(0,num_testing_input,1)
  for(i in 1:num_testing_input){
    testing_output[i]<-borehole(testing_input[i,])
  }
  
  # compute the MSE, average coverage and average length
  # out of sample MSE
  MSE_emulator <- sum((m3.predict$mean-testing_output)^2)/(num_testing_input)  
  
  # proportion covered by 95% posterior predictive credible interval
  prop_emulator <- length(which((m3.predict$lower95<=testing_output)
                   &(m3.predict$upper95>=testing_output)))/num_testing_input
  
  # average length of  posterior predictive credible interval
  length_emulator <- sum(m3.predict$upper95-m3.predict$lower95)/num_testing_input
  
  # output of sample prediction
  MSE_emulator
  prop_emulator
  length_emulator  
  # normalized RMSE
  sqrt(MSE_emulator/mean((testing_output-mean(output))^2 ))

Predicted PP GaSP class

Description

S4 class for the prediction of a PP GaSP model

Objects from the Class

Objects of this class are created and initialized with the function predict.ppgasp that computes the prediction on the PP GaSP model after the PP GaSP model has been constructed.

Slots

call:

call to predict.ppgasp function where the returned object has been created.

mean:

predictive mean for the testing inputs.

lower95:

lower bound of the 95% posterior credible interval.

upper95:

upper bound of the 95% posterior credible interval.

sd:

standard deviation of each testing_input.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

See Also

predict.ppgasp for more details about how to make predictions based on a ppgasp object.

Predictive robust GaSP class

Description

S4 class for the prediction of a Robust GaSP

Objects from the Class

Objects of this class are created and initialized with the function predict.rgasp that computes the prediction on Robust GaSP models after the Robust GaSP model has been constructed.

Slots

call:

call to predict.rgasp function where the returned object has been created.

mean:

predictive mean for the testing inputs.

lower95:

lower bound of the 95% posterior credible interval.

upper95:

upper bound of the 95% posterior credible interval.

sd:

standard deviation of each testing_input.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

See Also

predict.rgasp for more details about how to make predictions based on a rgasp object.

Setting up the robust GaSP model

Description

Setting up the robust GaSP model for estimating the parameters (if the parameters are not given).

Usage

  rgasp(design, response,trend=matrix(1,length(response),1),zero.mean="No",nugget=0,
    nugget.est=F,range.par=NA,method='post_mode',prior_choice='ref_approx',a=0.2,
    b=1/(length(response))^{1/dim(as.matrix(design))[2]}*(a+dim(as.matrix(design))[2]),
    kernel_type='matern_5_2',isotropic=F,R0=NA, 
    optimization='lbfgs', alpha=rep(1.9,dim(as.matrix(design))[2]),
    lower_bound=T,max_eval=max(30,20+5*dim(design)[2]),
    initial_values=NA,num_initial_values=2)
 

Arguments

Value

rgasp returns a S4 object of class rgasp (see rgasp-class).

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu, X. Wang and J.O. Berger (2018), Robust Gaussian stochastic process emulation, Annals of Statistics, 46(6A), 3038-3066.

M. Gu (2018), Jointly robust prior for Gaussian stochastic process in emulation, calibration and variable selection, arXiv:1804.09329.

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

M. Gu. and J.O. Berger (2016). Parallel partial Gaussian process emulation for computer models with massive output. Annals of Applied Statistics, 10(3), 1317-1347.

E.T. Spiller, M.J. Bayarri, J.O. Berger and E.S. Calder and A.K. Patra and E.B. Pitman, and R.L. Wolpert (2014), Automating emulator construction for geophysical hazard maps. SIAM/ASA Journal on Uncertainty Quantification, 2(1), 126-152.

J. Nocedal (1980), Updating quasi-Newton matrices with limited storage, Math. Comput., 35, 773-782.

D. C. Liu and J. Nocedal (1989), On the limited memory BFGS method for large scale optimization, Math. Programming, 45, p. 503-528.

Brent, R. (1973), Algorithms for Minimization without Derivatives. Englewood Cliffs N.J.: Prentice-Hall.

Examples

  library(RobustGaSP)
  #------------------------
  # a 3 dimensional example
  #------------------------
  # dimensional of the inputs
  dim_inputs <- 3    
  # number of the inputs
  num_obs <- 50       
  # uniform samples of design
  input <- matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  ##maximin lhd sample
  
  ####
  # outputs from the 3 dim dettepepel.3.data function
  
  output = matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]<-dettepepel.3.data (input[i,])
  }
  
  # use constant mean basis, with no constraint on optimization
  # and marginal posterior mode estimation
  m1<- rgasp(design = input, response = output, lower_bound=FALSE)
  
  # you can use specify the estimation as maximum likelihood estimation (MLE)
  m2<- rgasp(design = input, response = output, method='mle',lower_bound=FALSE)
  
  ##let's do some comparison on prediction
  n_testing=1000
  testing_input=matrix(runif(n_testing*dim_inputs),n_testing,dim_inputs)
  
  m1_pred=predict(m1,testing_input=testing_input)
  m2_pred=predict(m2,testing_input=testing_input)
  
  
  ##root of mean square error and interval
  test_output = matrix(0,n_testing,1)
  for(i in 1:n_testing){
    test_output[i]<-dettepepel.3.data (testing_input[i,])
  }
  
  ##root of mean square error
  sqrt(mean( (m1_pred$mean-test_output)^2))
  sqrt(mean( (m2_pred$mean-test_output)^2))
  sd(test_output)
  #---------------------------------------
  # a 1 dimensional example with zero mean
  #---------------------------------------


  input=10*seq(0,1,1/14)
  output<-higdon.1.data(input)
  #the following code fit a GaSP with zero mean by setting zero.mean="Yes"
  model<- rgasp(design = input, response = output, zero.mean="Yes")
  model
  
  testing_input = as.matrix(seq(0,10,1/100))
  model.predict<-predict(model,testing_input)
  names(model.predict)
  
  #########plot predictive distribution
  testing_output=higdon.1.data(testing_input)
  plot(testing_input,model.predict$mean,type='l',col='blue',
       xlab='input',ylab='output')
  polygon( c(testing_input,rev(testing_input)),c(model.predict$lower95,
        rev(model.predict$upper95)),col =  "grey80", border = FALSE)
  lines(testing_input, testing_output)
  lines(testing_input,model.predict$mean,type='l',col='blue')
  lines(input, output,type='p')
  
  ## mean square erros
  mean((model.predict$mean-testing_output)^2)


  #-----------------------------------
  # a 2 dimensional example with trend
  #-----------------------------------
  # dimensional of the inputs
  dim_inputs <- 2    
  # number of the inputs
  num_obs <- 20       
  
  # uniform samples of design
  input <-matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  # maximin lhd sample
  
  # outputs from a 2 dim function
  
  output <- matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]<-limetal.2.data (input[i,])
  }
  
  ####trend or mean basis
  X<-cbind(rep(1,num_obs), input )
  
  
  # use constant mean basis with trend, with no constraint on optimization
  m2<- rgasp(design = input, response = output,trend =X,  lower_bound=FALSE)

  show(m2)      # show this rgasp object 
  
  m2@beta_hat       # estimated inverse range parameters
  m2@theta_hat      # estimated trend parameters

  #--------------------------------------------------------------------------------------
  # an 8 dimensional example using only a subset inputs and a noise with unknown variance
  #--------------------------------------------------------------------------------------
  set.seed(1)
  # dimensional of the inputs
  dim_inputs <- 8    
  # number of the inputs
  num_obs <- 50       
  
  # uniform samples of design
  input <-matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  # maximin lhd sample
  
  # rescale the design to the domain
  input[,1]<-0.05+(0.15-0.05)*input[,1];
  input[,2]<-100+(50000-100)*input[,2];
  input[,3]<-63070+(115600-63070)*input[,3];
  input[,4]<-990+(1110-990)*input[,4];
  input[,5]<-63.1+(116-63.1)*input[,5];
  input[,6]<-700+(820-700)*input[,6];
  input[,7]<-1120+(1680-1120)*input[,7];
  input[,8]<-9855+(12045-9855)*input[,8];
  
  # outputs from the 8 dim Borehole function
  
  output=matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]=borehole(input[i,])
  }
  
  
    
    
  
  # use constant mean basis with trend, with no constraint on optimization
  m3<- rgasp(design = input[,c(1,4,6,7,8)], response = output, 
            nugget.est=TRUE, lower_bound=FALSE)

  m3@beta_hat       # estimated inverse range parameters
  m3@nugget     

Robust GaSP class

Description

S4 class for Robust GaSP if the range and noise-variance ratio parameters are given and/or have been estimated.

Objects from the Class

Objects of this class are created and initialized with the function rgasp that computes the calculations needed for setting up the analysis.

Slots

p:

Object of class integer. The dimensions of the inputs.

num_obs:

Object of class integer. The number of observations.

input:

Object of class matrix with dimension n x p. The design of experiments.

output:

Object of class matrix with dimension n x 1. The Observations or output vector.

X:

Object of class matrix of with dimension n x q. The mean basis function, i.e. the trend function.

zero_mean:

A character to specify whether the mean is zero or not. "Yes" means it has zero mean and "No"" means the mean is not zero.

q:

Object of class integer. The number of mean basis.

LB:

Object of class vector with dimension p x 1. The lower bound for inverse range parameters beta.

beta_initial:

Object of class vector with the initial values of inverse range parameters p x 1.

beta_hat:

Object of class vector with dimension p x 1. The inverse-range parameters.

log_post:

Object of class numeric with the logarithm of marginal posterior.

R0:

Object of class list of matrices where the j-th matrix is an absolute difference matrix of the j-th input vector.

theta_hat:

Object of class vector with dimension q x 1. The the mean (trend) parameter.

L:

Object of class matrix with dimension n x n. The Cholesky decomposition of the correlation matrix R, i.e.

L\%*\%t(L)=R

sigma2_hat:

Object of the class numeric. The estimated variance parameter.

LX:

Object of the class matrix with dimension q x q. The Cholesky decomposition of the correlation matrix

t(X)\%*\%R^{-1}\%*\%X

CL:

Object of the class vector used for the lower bound and the prior.

nugget:

A numeric object used for the noise-variance ratio parameter.

nugget.est:

A logical object of whether the nugget is estimated (T) or fixed (F).

kernel_type:

A vector of character to specify the type of kernel to use.

alpha:

Object of class vector with dimension p x 1 for the roughness parameters in the kernel.

method:

Object of class character to specify the method of parameter estimation. There are three values: post_mode, mle and mmle.

isotropic:

Object of class logical to specify whether the kernel is isotropic.

call:

The call to rgasp function to create the object.

Methods

show

Prints the main slots of the object.

predict

See predict.

Note

The response output must have one dimension. The number of observations in input must be equal to the number of experiments output.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

See Also

RobustGaSP for more details about how to create a RobustGaSP object.

Search for the default lower bound of range parameters.

Description

Function to find the values to construct the default lower bound of range parameters.

Usage

search_LB_prob(param, R0, COND_NUM_UB, p, kernel_type, alpha, nugget)

Arguments

Value

A vector of values used in constructing the default lower bound of range parameters.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

Mengyang Gu. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

Product correlation matrix with the product form

Description

Function to construct the product correlation matrix with the product form.

Usage

separable_kernel(R0, beta, kernel_type, alpha)

Arguments

Value

The product correlation matrix with the product form.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Product correlation matrix with the product form

Description

Function to construct the product correlation matrix with the product form. The kernel can be different for each coordinate of the input.

Usage

separable_multi_kernel(R0, beta, kernel_type, alpha)

Arguments

Value

The product correlation matrix with the product form.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Show Robust GaSP object

Description

Function to print Robust GaSP models after the Robust GaSP model has been constructed.

Usage

## S4 method for signature 'rgasp'
show(object)

Arguments

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Examples

  #------------------------
  # a 3 dimensional example
  #------------------------
  # dimensional of the inputs
  dim_inputs <- 3    
  # number of the inputs
  num_obs <- 30       
  # uniform samples of design
  input <- matrix(runif(num_obs*dim_inputs), num_obs,dim_inputs) 
  
  # Following codes use maximin Latin Hypercube Design, which is typically better than uniform
  # library(lhs)
  # input <- maximinLHS(n=num_obs, k=dim_inputs)  ##maximin lhd sample
  
  ####
  # outputs from the 3 dim dettepepel.3.data function
  
  output = matrix(0,num_obs,1)
  for(i in 1:num_obs){
    output[i]<-dettepepel.3.data (input[i,])
  }
  
  # use constant mean basis, with no constraint on optimization
  m1<- rgasp(design = input, response = output, lower_bound=FALSE)
  
  # the following use constraints on optimization
  # m1<- rgasp(design = input, response = output, lower_bound=TRUE)
  
  # the following use a single start on optimization
  # m1<- rgasp(design = input, response = output, lower_bound=FALSE)
  
  
  show(m1)

Show parllel partial Gaussian stochastic process (PP GaSP) object

Description

Function to print the PP GaSP model after the PP GaSP model has been constructed.

Usage

## S4 method for signature 'ppgasp'
show(object)

Arguments

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

Examples

  library(RobustGaSP)
  
  ###PP GaSP model for the humanity model
  data(humanity_model)
  ##pp gasp
  m.ppgasp=ppgasp(design=humanity.X,response=humanity.Y,nugget.est= TRUE)
  show(m.ppgasp)

Sample for Robust GaSP model

Description

Function to sample Robust GaSP after the Robust GaSP model has been constructed.

Usage

## S4 method for signature 'rgasp'
simulate(object, testing_input, num_sample=1,
testing_trend= matrix(1,dim(testing_input)[1],1),
r0=NA,rr0=NA,sample_data=T,...)

Arguments

Value

The returned value is a matrix where each column is a sample on the prespecified inputs.

Author(s)

Mengyang Gu [aut, cre], Jesus Palomo [aut], James Berger [aut]

Maintainer: Mengyang Gu <mengyang@pstat.ucsb.edu>

References

M. Gu. (2016). Robust uncertainty quantification and scalable computation for computer models with massive output. Ph.D. thesis. Duke University.

Examples

  #------------------------
  # a 1 dimensional example
  #------------------------
  
###########1dim higdon.1.data 
p1 = 1     ###dimensional of the inputs
dim_inputs1 <- p1
n1 = 15   ###sample size or number of training computer runs you have 
num_obs1 <- n1
input1 = 10*matrix(runif(num_obs1*dim_inputs1), num_obs1,dim_inputs1) ##uniform
#####lhs is better
#library(lhs)
#input1 = 10*maximinLHS(n=num_obs1, k=dim_inputs1)  ##maximin lhd sample
output1 = matrix(0,num_obs1,1)
for(i in 1:num_obs1){
  output1[i]=higdon.1.data (input1[i])
}




m1<- rgasp(design = input1, response = output1, lower_bound=FALSE)

#####locations to samples
testing_input1 = seq(0,10,1/50) 
testing_input1=as.matrix(testing_input1)
#####draw 10 samples
m1_sample=simulate(m1,testing_input1,num_sample=10)

#####plot these samples
matplot(testing_input1,m1_sample, type='l',xlab='input',ylab='output')
lines(input1,output1,type='p')