# rcalibration: Setting up the robust Calibration model In RobustCalibration: Robust Calibration of Imperfect Mathematical Models

## Description

Setting up the Calibration model for estimating the parameters via MCMC with or without a discrepancy function.

## Usage

 ```1 2 3 4 5 6 7 8``` ``` rcalibration(design, observations, p_theta=NULL,X=matrix(0,dim(design),1), have_trend=FALSE,simul_type=0, input_simul=NULL, output_simul=NULL,simul_nug=FALSE, math_model=NULL,theta_range=NULL, sd_proposal=c(rep(0.05,p_theta), rep(0.25,dim(design)),0.25),S=10000,S_0=1000, discrepancy_type='S-GaSP', kernel_type='matern_5_2', tilde_lambda=1/2, a=1/2-dim(design),b=1, alpha=rep(1.9,dim(design)), output_weights=rep(1,dim(design))) ```

## Arguments

 `design` a matrix of observed inputs. `observations` a vector of experimental data. `p_theta` an integer about the number of parameters, which should be specified by the user. `X` a matrix of the mean/trend discrepancy between the reality and math model. `have_trend` a bool value meaning whether we assume a mean/trend discrepancy function. `simul_type` an integer about the math model/simulator. If the simul_type is 0, it means we use RobustGaSP R package to build an emulator for emulation. If the simul_type is 1, it means the function of the math model is given by the user. When simul_type is 2 or 3, the mathematical model is the geophyiscal model for Kilauea Volcano. If the simul_type is 2, it means it is for the ascending mode InSAR data; if the simul_type is 3, it means it is for the descending mode InSAR data. `input_simul` an D x (p_x+p_theta) matrix of design for emulating the math model. It is only useful if simul_type is 0, meaning that we emulate the output of the math model. `output_simul` a D dimensional vector of the math model runs on the design (input_simul). It is only useful if simul_type is 0, meaning that we emulate the output of the math model. `simul_nug` a bool value meaning whether we have a nugget for emulating the math model/simulator. If the math model is stochastic, we often need a nugget. If simul_Nug is TRUE, it means we have a nugget for the emulator. If simul_Nug is FALSE, it means we do not have a nugget for the emulator. `math_model` a function of the math model provided by the user. It is only useful if simul_type is 1, meaning that we know the math model and it can be computed fast. If the evaluation the math model is computationally slow, one should set simul_type to be 0 to emulate the math model. `theta_range` a p_theta x 2 matrix of the range of the calibration parameters. The first column is the lower bound and the second column is the upper bound. It should be specified by the user if the simul_type is 0. `sd_proposal` a vector of the standard deviation of the proposal distribution in MCMC. `S` an integer about about how many posterior samples to run. `S_0` an integer about about the number of burn-in samples. `discrepancy_type` characters about the type of the discrepancy. If it is 'no-discrepancy', it means no discrepancy function. If it is 'GaSP', it means the GaSP model for the discrepancy function. If it is 'S-GaSP', it means the S-GaSP model for the discrepancy function. `kernel_type` characters about the type of the discrepancy.type of kernel. `matern_3_2` and `matern_5_2` are `Matern kernel` with roughness parameter 3/2 and 5/2 respectively. `pow_exp` is power exponential kernel with roughness parameter alpha. If `pow_exp` is to be used, one needs to specify its roughness parameter alpha. `tilde_lambda` a numeric value about how close the math model to the reality in squared distance when the S-GaSP model is used for modeling the discrepancy. `a` a scalar of the prior parameter. `b` a scalar of the prior parameter. `alpha` a numeric parameter for the roughness in the kernel. `output_weights` a vector of the weights of the outputs.

## Value

`rcalibration` returns an S4 object of class `rcalibration` (see `rcalibration-class`).

## 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.

K. R. Anderson and M. P. Poland (2016), Bayesian estimation of magma supply, storage, and eroption rates using a multiphysical volcano model: Kilauea volcano, 2000-2012.. Eath and Planetary Science Letters, 447, 161-171.

K. R. Anderson and M. P. Poland (2017), Abundant carbon in the mantle beneath Hawaii. Nature Geoscience, 10, 704-708.

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``` ```library(RobustCalibration) #------------------------------------------------------------- # an example with multiple local maximum of minimum in L2 loss #------------------------------------------------------------- ## the reality test_funct_eg1<-function(x){ x*cos(3/2*x)+x } ## obtain 25 data from the reality plus a noise set.seed(1) ## 10 data points are very small, one may want to add more data n=15 input=seq(0,5,5/(n-1)) input=as.matrix(input) output=test_funct_eg1(input)+rnorm(length(input),mean=0,sd=0.1) num_obs=n=length(output) ## the math model math_model_eg1<-function(x,theta){ sin(theta*x)+x } ##fit the S-GaSP model for the discrepancy ##one can choose the discrepancy_type to GaSP, S-GaSP or no discrepancy ##p_theta is the number of parameters to calibrate and user needs to specifiy ##one may also want to change the number of posterior samples by change S and S_0 p_theta=1 model_sgasp=rcalibration(design=input, observations=output, p_theta=p_theta,simul_type=1, math_model=math_model_eg1,theta_range=matrix(c(0,3),1,2) ,S=10000,S_0=2000,discrepancy_type='S-GaSP') ##if the acceptance rate is too low or two high, one can adjust sd_proposal, e.g. #model_sgasp=rcalibration(design=input, observations=output, p_theta=1,simul_type=1, # sd_proposal=c(rep(0.02,p_theta),rep(0.2,dim(input)),0.2) # math_model=math_model_eg1,theta_range=matrix(c(0,3),1,2) # ,S=10000,S_0=2000,discrepancy_type='S-GaSP') ##posterior samples of calibration parameter and value plot(model_sgasp@post_sample[,1],type='l',xlab='num',ylab=expression(theta)) plot(model_sgasp@post_value,type='l',xlab='num',ylab='posterior value') show(model_sgasp) ##one may want to fit a a model with an estimated baseline mean discrepancy by setting ##X=matrix(1,dim(input_stack),1),have_trend=TRUE model_sgasp_with_mean=rcalibration(design=input, observations=output, p_theta=1,simul_type=1, X=matrix(1,dim(input),1),have_trend=TRUE, math_model=math_model_eg1,theta_range=matrix(c(0,3),1,2), S=10000,S_0=2000,discrepancy_type='S-GaSP') show(model_sgasp_with_mean) ##posterior samples of calibration parameter and value plot(model_sgasp_with_mean@post_sample[,1],type='l',xlab='num',ylab=expression(theta)) plot(model_sgasp_with_mean@post_value,type='l',xlab='num',ylab='posterior value') ## Not run: #------------------------------------------------------------- # an example with multiple local maximum of minimum in L2 loss # for combing the emulator #------------------------------------------------------------- ## the reality test_funct_eg1<-function(x){ x*cos(3/2*x)+x } ## obtain 20 data from the reality plus a noise set.seed(1) n=20 input=seq(0,5,5/(n-1)) input=as.matrix(input) output=test_funct_eg1(input)+rnorm(length(input),mean=0,sd=0.05) num_obs=n=length(output) ## the math model math_model_eg1<-function(x,theta){ sin(theta*x)+x } ##let's build an emulator for the case if the math model is too slow # let's say we can only run the math model n_design times n_design=80 design_simul=matrix(runif(n_design*2),n_design,2) design_simul[,1]=5*design_simul[,1] ##the first one is the observed input x design_simul[,2]=3*design_simul[,2] ##the second one is the calibration parameter output_simul=math_model_eg1(design_simul[,1],design_simul[,2]) ##this is a little slow compared with the previous model model_sgasp_emulator=rcalibration(design=input, observations=output, p_theta=1,simul_type=0, input_simul=design_simul, output_simul=output_simul, theta_range=matrix(c(0,3),1,2), S=10000,S_0=2000,discrepancy_type='S-GaSP') ##now the output is a list show(model_sgasp_emulator) ##here is the plot plot(model_sgasp_emulator@post_sample[,1],type='l',xlab='num',ylab=expression(theta)) plot(model_sgasp_emulator@post_value,type='l',xlab='num',ylab='posterior value') ## End(Not run) ```

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