predict.ppgasp {RobustGaSP} | R Documentation |
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 |
testing_input |
a matrix containing the inputs where the |
testing_trend |
a matrix of mean/trend for prediction. |
r0 |
the distance between input and testing input. If the value
is |
interval_data |
a boolean value. If |
outasS3 |
a boolean parameter indicating whether the output of the function should be as an |
loc_index |
specified coodinate index of the prediction. The default value is |
... |
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: |
|
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 |
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 |
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()