# predict.rcalibration: Prediction for the robust calibration model In RobustCalibration: Robust Calibration of Imperfect Mathematical Models

## Description

Function to make prediction on Robust Calibration models after the rcalibration class has been constructed.

## Usage

 1 2 3 4 5 6 ## S4 method for signature 'rcalibration' predict(object, testing_input,X_testing=matrix(0,dim(testing_input)[1],1), n_thinning=10, testing_output_weights=rep(1,dim(testing_input)[1]), interval_est=NULL,interval_data=F, math_model=NULL,...) 

## Arguments

 object an object of class rcalibration. testing_input a matrix containing the inputs where the predict is to perform prediction. X_testing a matrix of mean/trend for prediction. n_thinning number of points thinning the MCMC posterior samples. testing_output_weights the weight of testing outputs. interval_est a vector for the the posterior credible interval. If interval_est is NULL, we do not compute the posterior credible interval. It can be specified as a vector of values ranging from zero to one. E.g. if interval_est=c(0.025, 0.975), the 95 posterior credible interval will be computed. interval_data a bool value to decide whether the experimental noise is included for computing the posterior credible interval. math_model a function for the math model to be calibrated. ... extra arguments to be passed to the function (not implemented yet).

## Value

The returned value is a S4 CLass predictobj.rcalibration.

## Author(s)

Mengyang Gu [aut, cre]

Maintainer: Mengyang Gu <mgu6@jhu.edu>

## References

A. O'Hagan and M. C. Kennedy (2001), Bayesian calibration of computer models, Journal of the Royal Statistical Society: Series B (Statistical Methodology, 63, 425-464.

Bayarri, Maria J and Berger, James O and Paulo, Rui and Sacks, Jerry and Cafeo, John A and Cavendish, James and Lin, Chin-Hsu and Tu, Jian (2007) A framework for validation of computer models. Technometrics. 49, 138–154.

M. Gu (2016), Robust Uncertainty Quantification and Scalable Computation for Computer Models with Massive Output, Ph.D. thesis., Duke University.

M. Gu and L. Wang (2017) Scaled Gaussian Stochastic Process for Computer Model Calibration and Prediction. arXiv preprint arXiv:1707.08215.

M. Gu (2018) Jointly Robust Prior for Gaussian Stochastic Process in Emulation, Calibration and Variable Selection . arXiv preprint 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 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 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 #------------------------------------------------------------------------------ # Example: an example used in Susie Bayarri et. al. 2007 Technometrics paper #------------------------------------------------------------------------------ ##reality test_funct_eg1<-function(x){ 3.5*exp(-1.7*x)+1.5 } ##math model math_model_eg1<-function(x,theta){ 5*exp(-x*theta) } ## noise observations (sampled from reality + independent Gaussian noises) ## each has 3 replicates input=c(rep(.110,3),rep(.432,3),rep(.754,3),rep(1.077,3),rep(1.399,3),rep(1.721,3), rep(2.043,3),rep(2.366,3),rep(2.688,3),rep(3.010,3)) output=c(4.730,4.720,4.234,3.177,2.966,3.653,1.970,2.267,2.084,2.079,2.409,2.371,1.908,1.665,1.685, 1.773,1.603,1.922,1.370,1.661,1.757,1.868,1.505,1.638,1.390,1.275,1.679,1.461,1.157,1.530) ## calculating the average or the stack data n_stack=length(output)/3 output_stack=rep(0,n_stack) input_stack=rep(0,n_stack) for(j in 1:n_stack){ output_stack[j]=mean(output[ ((j-1)*3+1):(3*j)]) input_stack[j]=mean(input[ ((j-1)*3+1):(3*j)]) } output_stack=as.matrix(output_stack) input_stack=as.matrix(input_stack) ## plot the output and stack output #plot(input,output,pch=16,col='red') #lines(input_stack,output_stack,pch=16,col='blue',type='p') ## fit model using S-GaSP for the discrepancy ## one can change S and S_0 for the number of posterior and burn-in samples ## Normallly you may need a larger number of posterior sample ## you can set S=50000 and S_0=5000 ## one may also change the sd of the proposal distribution using sd_proposal model_sgasp=rcalibration(design=input_stack, observations=output_stack, p_theta=1,simul_type=1, math_model=math_model_eg1,theta_range=matrix(c(0,10),1,2), S=10000,S_0=2000,discrepancy_type='S-GaSP') # one can fit the GaSP model for discrepancy function by discrepancy_type='GaSP' # one can fit a model without the discrepancy function by discrepancy_type='no-discrepancy' ## posterior of the calibration parameter #plot(model_sgasp@post_sample[,1],type='l',xlab='num',ylab=expression(theta)) show(model_sgasp) ## ## test data set testing_input=as.matrix(seq(0,6,0.02)) ##perform prediction prediction_sgasp=predict(model_sgasp,testing_input,math_model=math_model_eg1, interval_est=c(0.025,0.975),interval_data=TRUE, n_thinning =20 ) ##real test output testing_output=test_funct_eg1(testing_input) ##the prediction by S-GaSP min_val=min(prediction_sgasp@mean,prediction_sgasp@interval,output,testing_output) max_val=max(prediction_sgasp@mean,prediction_sgasp@interval,output,testing_output) plot(testing_input,prediction_sgasp@mean,type='l',col='blue',xlab='x',ylab='y', ylim=c(min_val,max_val) ) lines(testing_input,prediction_sgasp@interval[,1],col='blue',lty=2) lines(testing_input,prediction_sgasp@interval[,2],col='blue',lty=2) lines(input,output,type='p') lines(testing_input,prediction_sgasp@math_model_mean,col='blue',lty=3) lines(testing_input,testing_output,type='l') legend("topright", legend=c("reality", "predictive mean","95 percent posterior credible interval", "predictive mean of the math model"), col=c("black", "blue","blue","blue"), lty=c(1,1,2,3),cex=.6) ## MSE if the math model and discrepancy are used for prediction mean((testing_output-prediction_sgasp@mean)^2) ## MSE if the math model is used for prediction mean((testing_output-prediction_sgasp@math_model_mean)^2) ################################## #the example with a mean structure ################################## ##now let's fit model with mean model_sgasp_with_mean=rcalibration(design=input_stack, observations=output_stack, p_theta=1,X=matrix(1,dim(input_stack)[1],1), have_trend=TRUE,simul_type=1, math_model=math_model_eg1, theta_range=matrix(c(0,10),1,2), S=10000,S_0=2000, discrepancy_type='S-GaSP') #posterior #plot(model_sgasp_with_mean@post_sample[,1],type='l',xlab='num',ylab=expression(theta)) show(model_sgasp_with_mean) ## test data set testing_input=as.matrix(seq(0,6,0.02)) ##in this way we don't output the interval estimation ## one can have interval using interval_est=c(0.025,0.975),interval_data=TRUE prediction_sgasp_with_mean=predict(model_sgasp_with_mean,testing_input, X_testing=matrix(1,dim(testing_input)[1],1), math_model=math_model_eg1,n_thinning = 20) ##plot for the S-GaSP ##for this example, with a mean structure, it fits much better min_val=min(prediction_sgasp_with_mean@mean,output,testing_output) max_val=max(prediction_sgasp_with_mean@mean,output,testing_output) plot(testing_input,prediction_sgasp_with_mean@mean,type='l',col='blue',xlab='x', ylab='y',ylim=c(min_val,max_val) ) #lines(testing_input,prediction_sgasp_with_mean@interval[,1],col='blue',lty=2) #lines(testing_input,prediction_sgasp_with_mean@interval[,2],col='blue',lty=2) lines(input,output,type='p') lines(testing_input,prediction_sgasp_with_mean@math_model_mean,col='blue',lty=3) lines(testing_input,testing_output,type='l') legend("topright", legend=c("reality", "predictive mean", "predictive mean of the math model"), col=c("black", "blue","blue"), lty=c(1,1,3),cex=.6) ## MSE if the math model and discrepancy are used for prediction mean((testing_output-prediction_sgasp_with_mean@mean)^2) ## MSE if the math model is used for prediction mean((testing_output-prediction_sgasp_with_mean@math_model_mean)^2) ## Not run: #------------------------------------------------------------- #the example with the emulator #------------------------------------------------------------- n_design=80 design_simul=matrix(runif(n_design*2),n_design,2) #library(lhs) #design_simul=maximinLHS(n=n_design,k=2) design_simul[,1]=6*design_simul[,1] ##the first one is the observed input x design_simul[,2]=10*design_simul[,2] ##the second one is the calibration parameter \theta output_simul=math_model_eg1(design_simul[,1],design_simul[,2]) ##this is a little slow compared with the previous model model_sgasp_with_mean_emulator=rcalibration(design=input_stack, observations=output_stack, p_theta=1,simul_type=0, have_trend=T,X=matrix(1,dim(input_stack)[1],1), input_simul=design_simul, output_simul=output_simul, theta_range=matrix(c(0,10),1,2), S=10000,S_0=2000,discrepancy_type='S-GaSP') ##now the output is a list show(model_sgasp_with_mean_emulator) ##here is the plot plot(model_sgasp_with_mean_emulator@post_sample[,4],type='l',xlab='num',ylab=expression(theta)) plot(model_sgasp_with_mean_emulator@post_value,type='l',xlab='num',ylab='posterior value') prediction_sgasp_with_mean_emulator=predict(model_sgasp_with_mean_emulator,testing_input, X_testing=matrix(1,dim(testing_input)[1],1), interval_est=c(0.025,0.975), interval_data=TRUE) ##for this example, with a mean structure, it fits much better min_val=min(prediction_sgasp_with_mean_emulator@mean,output,testing_output, prediction_sgasp_with_mean_emulator@math_model_mean) max_val=max(prediction_sgasp_with_mean_emulator@mean,output,testing_output, prediction_sgasp_with_mean_emulator@math_model_mean) plot(testing_input,prediction_sgasp_with_mean_emulator@mean,type='l',col='blue',xlab='x', ylab='y',ylim=c(min_val,max_val) ) #lines(testing_input,prediction_sgasp_with_mean@interval[,1],col='blue',lty=2) #lines(testing_input,prediction_sgasp_with_mean@interval[,2],col='blue',lty=2) lines(input,output,type='p') lines(testing_input,prediction_sgasp_with_mean_emulator@math_model_mean,col='blue',lty=3) lines(testing_input,testing_output,type='l') legend("topright", legend=c("reality", "predictive mean", "predictive mean of the math model"), col=c("black", "blue","blue"), lty=c(1,1,3),cex=.6) ## MSE if the math model and discrepancy are used for prediction mean((testing_output-prediction_sgasp_with_mean_emulator@mean)^2) ## MSE if the math model is used for prediction mean((testing_output-prediction_sgasp_with_mean_emulator@math_model_mean)^2) ## End(Not run) 

RobustCalibration documentation built on May 2, 2019, 9:36 a.m.