model: Generates `model.class` objects.

Description

model is a function that generates a calibration model and the associated likelihood.

Usage

model(code, X, Yexp, model = "model1", ...)

Arguments

code

the computational code (function of X and theta)

the matrix of the forced variables

Yexp

the vector of the experiments

model

string of the model chosen ("model1","model2","model3","model4") by default "model1" is chosen. See details for further clarifications

...

additional options (see details section)

Value

model returns a model.class object. This class contains two main methods:

plot(model,x) this method generates the plot for a new $Θ$ , $Θ_{D}$ (for model3 and model4), $σ^{2}$ . The parameter values need to be added to the model with the pipe %<%. The argument x represents the x-axis and have to be specified to produce a plot.
print(model) this method presents several pieces of information about the model.

Details

The different statistical models are:

Model2: $f o r i i n [1, . . ., n] Y e x p_{i} = F (x_{i}, Θ) + ϵ (x_{i})$
Model3: $f o r i i n [1, . . ., n] Y e x p_{i} = f (x_{i}, Θ) + δ (x_{i}) + ϵ (x_{i})$
Model4: $f o r i i n [1, . . ., n] Y e x p_{i} = F (x_{i}, Θ) + δ (x_{i}) + ϵ (x_{i})$

where $f o r i i n [1, \dots, n] ϵ (x_{i}) N (0, σ^{2})$ , $F (., .) P G (m_{1} (., .), c_{1} (., .), (., .))$ and $δ (.) P G (m_{2} (.), c_{2} (., .))$ . There is four kind of models in calibration. They are properly defined in [1].

To establish a Gaussian process three options are available:

opt.gp is an option list containing the parameters to establish the surrogate (only for model2 and model4).
- type type of the chosen kernel (value by default "matern5_2") from km function
- DOE design of experiments for the surrogate (default value NULL). If NULL the DOE is automatically generated with the opt.emul option.
opt.emul is an option list containing characteristics about emulation option (only for model2 and model4).
- p the number of parameter in the model (default value 1)
- n.emul the number of points for contituing the Design Of Experiments (DOE) (default value 100)
- binf the lower bound of the parameter vector (default value 0)
- bsup the upper bound of the parameter vector (default value 1)
opt.sim is an option list containing the design and corresponding outputs of the code, in the case where no numerical code is available (only for model2 and model4).
- Ysim Output of the code
- DOEsim DOE corresponding to the output of the code

To add a discrepancy in the model, the option opt.disc must be added:

opt.disc is an option list containing characteristics on the discrepancy (only for model3 and model4)
- kernel.type see kernel.fun for further details

References

[1] CARMASSI, Mathieu, BARBILLON, Pierre, KELLER, Merlin, et al. Bayesian calibration of a numerical code for prediction. arXiv preprint arXiv:1801.01810, 2018.

Examples

Run this code

# NOT RUN {
###### The code to calibrate
X <- cbind(seq(0,1,length.out=5),seq(0,1,length.out=5))
code <- function(X,theta)
{
  return((6*X[,1]*theta[2]-2)^2*theta[1]*sin(theta[3]*X[,2]-4))
}
Yexp <- code(X,c(1,1,11))+rnorm(5,0,0.1)

###### For the first model
### Generate the model
model1 <- model(code,X,Yexp,"model1")
### Plot the results with the first column of X
model1 %<% list(theta=c(1,1,11),var=0.01)
plot(model1,X[,1],CI="err")

### Summury of the foo class generated
print(model1)

###### For the second model
### code function is available, no DOE generated upstream
binf <- c(0.9,0.9,10.5)
bsup <- c(1.1,1.1,11.5)
opt.gp <- list(type="matern5_2", DOE=NULL)
opt.emul <- list(p=3,n.emul=150,binf=binf,bsup=bsup,type="maximinLHS")
model2 <- model(code,X,Yexp,"model2",
                opt.gp=opt.gp,
                opt.emul=opt.emul)
model2 %<% list(theta=c(1,1,11),var=0.1)
### Plot the model
plot(model2,X[,1])

### code function is available and use a specific DOE
DOE <- DiceDesign::lhsDesign(20,5)$design
DOE[,3:5] <- unscale(DOE[,3:5],binf,bsup)

opt.gp <- list(type="matern5_2", DOE=DOE)
model2 <- model(code,X,Yexp,"model2",
                opt.gp=opt.gp)
model2 %<% list(theta=c(1,1,11),var=0.1)
plot(model2,X[,1])

### no code function but DOE and code output available
Ysim <- matrix(nr=20,nc=1)
for (i in 1:20)
{
  covariates <- as.matrix(DOE[i,1:2])
  dim(covariates) <- c(1,2)
  Ysim[i] <- code(covariates,DOE[i,3:5])
}

opt.sim <- list(Ysim=Ysim,DOEsim=DOE)
opt.gp <- list(type="matern5_2", DOE=NULL)
model2 <- model(code=NULL,X,Yexp,"model2",
                opt.gp=opt.gp,
                opt.sim=opt.sim)
model2 %<% list(theta=c(1,1,11),var=0.1)
plot(model2,X[,1])

###### For the third model
model3 <- model(code,X,Yexp,"model3",opt.disc=list(kernel.type="gauss"))
model3 %<% list(theta=c(1,1,11),thetaD=c(20,0.5),var=0.1)
plot(model3,X[,1],CI="err")
print(model3)

# }
# NOT RUN {
# }

Run the code above in your browser using DataLab

Get 50% off unlimited learning