predict.ppgasp: Prediction for PP GaSP model

Description Usage Arguments Value Author(s) References Examples

Description

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

Usage

1
2
3
## S4 method for signature 'ppgasp'
predict(object, testing_input,  
                   testing_trend= matrix(1,dim(testing_input)[1],1),outasS3 = T, ...)

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.

outasS3

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

...

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 <[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. (2016). Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output. Ph.D. thesis. Duke University.

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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
  library(RobustGaSP)
  #----------------------------------
  # an example of environmental model
  #----------------------------------
  
  set.seed(1)
  n=100
  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=400
  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(3,60,3)  )
  }
  
  ##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)
  
  for(i in 1:n_star){
    testing_output[i,]=environ.4.data(testing_input[i,],s=seq(0.15,3,0.15
    ),t=seq(3,60,3) )
  }
  
  ##we do a transformation of the output 
  ##one can change the number of initial values to test
  log_output_1=log(output+1)
  m.ppgasp=ppgasp(design=input,response=log_output_1,kernel_type
                  ='pow_exp',num_initial_values=2)
  m_pred.ppgasp=predict(m.ppgasp,testing_input)
  ##we transform back for the prediction
  m_pred_ppgasp_mean=exp(m_pred.ppgasp$mean)-1
  ##mean squared error
  mean( (m_pred_ppgasp_mean-testing_output)^2)
  ##variance of the testing outputs
  var(as.numeric(testing_output))
  
  ##makes plots for the testing 
  par(mfrow=c(1,2))
  t=seq(3,60,3)
  s=seq(0.15,3,0.15)
  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_mean[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()

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