model: Generates 'model.class' objects.

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/functions.R

Description

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

Usage

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

Arguments

code

the computational code (function of X and theta)

X

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)

Details

The different statistical models are:

where for i in [1,…,n] ε(x_i)~N(0,σ^2), F(.,.)~PG(m_1(.,.),c_1{(.,.),(.,.)}) and δ(.)~PG(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:

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

Value

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

Author(s)

M. Carmassi

References

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

See Also

prior, calibrate, forecast, sequentialDesign

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


## End(Not run)

mathieucarmassi/calibrationcode documentation built on Aug. 14, 2019, 12:35 a.m.