ppgasp {RobustGaSP} | R Documentation |
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
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. |
range.par |
either |
method |
method of parameter estimation. |
prior_choice |
the choice of prior for range parameters and noise-variance parameters. |
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 |
kernel_type |
A vector specifying the type of kernels of each coordinate of the input. |
isotropic |
a boolean value. |
R0 |
the distance between inputs. If the value is |
optimization |
the method for numerically optimization of the kernel parameters. Currently three methods are implemented. |
alpha |
roughness parameters in the |
lower_bound |
boolean value. |
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 <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)