predict.ppgasp {RobustGaSP}

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

object

an object of class ppgasp.

testing_input

a matrix containing the inputs where the rgasp is to perform prediction.

testing_trend

a matrix of mean/trend for prediction.

r0

the distance between input and testing input. If the value is NA, it will be calculated later. It can also be specified by the user. If specified by user, it is either a matrix or list. The default value is NA.

interval_data

a boolean value. If T, the interval of the data will be calculated. Otherwise, the interval of the mean of the data will be calculted.

outasS3

a boolean parameter indicating whether the output of the function should be as an S3 object.

loc_index

specified coodinate index of the prediction. The default value is NA and prediction will be computed for all coordinates. If e.g. loc_index=c(3,5), it means the prediction will be computed on only the third and fifth coordinates, corresponding the coordinates of the third and fifth columns of the output matrix.

...

Extra arguments to be passed to the function (not implemented yet).

Value

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

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.

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

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 <[e-mail]>

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