ppgasp: Setting up the parallel partial GaSP model

Description Usage Arguments Value Author(s) References Examples

Description

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

Usage

1
2
3
4
5
6
7
   ppgasp(design,response,trend=matrix(1,dim(response)[1],1),zero.mean="No",nugget=0,
    nugget.est=F,range.par=NA,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',
    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

design

a matrix of inputs.

response

a matrix of outputs where each row is one runs of the computer model output.

trend

the mean/trend matrix of inputs. The default value is a vector of ones.

zero.mean

it has zero mean or not. The default value is FALSE meaning the mean is not zero. TRUE means the mean is zero.

nugget

numerical value of the nugget variance ratio. If nugget is equal to 0, it means there is either no nugget or the nugget is estimated. If the nugget is not equal to 0, it means a fixed nugget. The default value is 0.

nugget.est

boolean value. T means nugget should be estimated and F means nugget is fixed or not estimated. The default value is F F.

range.par

either NA or a vector. If it is NA, it means range parameters are estimated; otherwise range parameters are given. The default value is NA.

prior_choice

the choice of prior for range parameters and noise-variance parameters. ref_xi and ref_gamma means the reference prior with reference prior with the log of inverse range parameterization ξ or range parameterization γ. ref_approx uses the jointly robust prior to approximate the reference prior. The default choice is ref_approx.

a

prior parameters in the jointly robust prior. The default value is 0.2.

b

prior parameters in the jointly robust prior. The default value is n^{-1/p}(a+p) where n is the number of runs and p is the dimension of the input vector.

kernel_type

A vector specifying the type of kernels of each coordinate of the input. matern_3_2 and matern_5_2 are Matern correlation with roughness parameter 3/2 and 5/2 respectively. pow_exp is power exponential correlation with roughness parameter alpha. If pow_exp is to be used, one needs to specify its roughness parameter alpha. The default choice is matern_5_2.

alpha

roughness parameters in the pow_exp correlation functions. The default choice is a vector with each entry being 1.9.

lower_bound

boolean value. T means the default lower bounds of the inverse range parameters are used to constrained the optimization and F means the optimization is unconstrained. The default value is T and we also suggest to use F in various scenarios.

max_eval

the maximum number of steps to estimate the range and nugget parameters.

initial_values

a matrix of initial values of the kernel parameters to be optimized numerically, where each row of the matrix contains a set of the log inverse range parameters and the log nugget parameter.

num_initial_values

the number of initial values of the kernel parameters in optimization.

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 <[email protected]>

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.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
  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)

  ##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)

RobustGaSP documentation built on June 6, 2019, 1:02 a.m.