MLoptimConstrained: Constrained MLE Optimization

Description

Constrained MLE optimization for kernels defined by cliques using constrOptim

Usage

MLoptimConstrained(x, y, n.initialtries = 50, limits = NULL, eps.R = 1e-08, Cl, covtype = "gauss", eps.Var = 1e-06, MAXIT = 1000, iso = FALSE)

Arguments

a design matrix of input variables, number of columns should be number of variables

a vector of output variables of the same length as the columns of L

n.initialtries

number of random initial parameters for optimization, defaults to 50

limits

a list with items lower, upper containing boundaries for the covariance parameter vector theta, if NULL suitable bounds are computed from the range of x

eps.R

small positive number indicating the nugget effect added to the covariance matrix diagonalk, defaults to eps.R = 1e-08

list of cliques, can be obtained by function threshold

covtype

an optional character string specifying the covariance structure to be used, to be chosen between "gauss", "matern5_2", "matern3_2", "exp" or "powexp" (see DiceKriging

eps.Var

small positive number providing the limits for the alpha parameters in order to guarantee strict inequalities (0+eps.Var <= alpha="" <="1-esp.Var)," defaults="" to="" eps.Var = 1e-06

MAXIT

maximum number of iterations for optimization, defaults to MAXIT=1000

iso

boolean vector indicating for each clique if it is isotropic (TRUE) or anisotropic (FALSE), defaults to iso = FALSE (all cliques anisotropic)

Value

list of estimated parameter 'alpha' and 'theta' corresponding to the clique structure in 'Cl'

References

Muehlenstaedt, T.; Roustant, O.; Carraro, L.; Kuhnt, S. (2011) Data-driven Kriging models based on FANOVA-decomposition, Statistics and Computing.

Examples

Run this code

### example for ishigami function with cliques {1,3} and {2}
d <- 3
x <- matrix(runif(100*d,-pi,pi),nc=d)
y <- ishigami.fun(x)

Cl <- list(c(2), c(1,3))

# constrained ML optimation with kernel defined by the cliques
parameter <- MLoptimConstrained(x, y, Cl = Cl)

# prediction with the new model
xpred <- matrix(runif(500 * d,-pi,pi), ncol = d)
ypred <- yhat(xpred, x, y, parameter, Cl=Cl)
yexact <- ishigami.fun(xpred)

# rmse
sqrt(mean((ypred[,1]- yexact)^2))

# scatterplot
par(mfrow=c(1,1))
plot(yexact, ypred[,1], asp = 1)
abline(0, 1)

### compare to one single clique {1,2,3}
Cl <- list(c(1,2,3))

# constrained ML optimation with kernel defined by the cliques
parameter <- MLoptimConstrained(x, y, Cl = Cl)

# prediction with the new model
ypred <- yhat(xpred, x, y, parameter, Cl=Cl)

# rmse
sqrt(mean((ypred$mean- yexact)^2))

# scatterplot
par(mfrow=c(1,1))
plot(yexact, ypred$mean, asp = 1)
abline(0, 1)

### isotropic cliques

Cl <- list(c(2),c(1,3))
parameter <- MLoptimConstrained(x, y, Cl = Cl, iso=c(FALSE,TRUE))
ypred <- yhat(xpred, x, y, parameter, Cl=Cl, iso=c(FALSE,TRUE))
sqrt(mean((ypred$mean- yexact)^2))

# the same since first clique has length 1
parameter <- MLoptimConstrained(x, y, Cl = Cl, iso=c(TRUE,TRUE))
ypred <- yhat(xpred, x, y, parameter, Cl=Cl, iso=c(TRUE,TRUE))
sqrt(mean((ypred$mean- yexact)^2))