diff --git a/R/SCAD.R b/R/SCAD.R old mode 100755 new mode 100644 diff --git a/R/SCAD.derivative.R b/R/SCAD.derivative.R old mode 100755 new mode 100644 diff --git a/R/branin.R b/R/branin.R old mode 100755 new mode 100644 diff --git a/R/camelback.R b/R/camelback.R old mode 100755 new mode 100644 diff --git a/R/computeAuxVariables.R b/R/computeAuxVariables.R old mode 100755 new mode 100644 diff --git a/R/covStruct.create.R b/R/covStruct.create.R old mode 100755 new mode 100644 diff --git a/R/covStruct_XtraClass_Union.R b/R/covStruct_XtraClass_Union.R old mode 100755 new mode 100644 diff --git a/R/covStruct_XtraMethods.R b/R/covStruct_XtraMethods.R old mode 100755 new mode 100644 diff --git a/R/drop.reponse.R b/R/drop.reponse.R old mode 100755 new mode 100644 diff --git a/R/goldsteinPrice.R b/R/goldsteinPrice.R old mode 100755 new mode 100644 diff --git a/R/hartman3.R b/R/hartman3.R old mode 100755 new mode 100644 diff --git a/R/hartman6.R b/R/hartman6.R old mode 100755 new mode 100644 diff --git a/R/km.R b/R/km.R old mode 100755 new mode 100644 index 2d1268f..b1e373a --- a/R/km.R +++ b/R/km.R @@ -112,8 +112,14 @@ } model@param.estim <- TRUE + ## + optim.method <- tolower(optim.method) + if(!(optim.method %in% c("bfgs", "gen"))) { + stop("'optim.method' must be either 'BFGS' or 'gen'. Currently it is '", + optim.method, "'") + } model@optim.method <- optim.method - + ## if ((length(lower) == 0) || (length(upper) == 0)) { bounds <- covParametersBounds(model@covariance, design) if (length(lower) == 0) lower <- bounds$lower @@ -129,7 +135,7 @@ model@upper <- as.numeric(upper) model@parinit <- as.numeric(parinit) - if (optim.method == "BFGS") { + if (optim.method == "bfgs") { if (length(control$pop.size) == 0) control$pop.size <- 20 control$pop.size <- max(control$pop.size, multistart) if (identical(control$trace, FALSE)) control$trace <- 0 @@ -220,7 +226,7 @@ model@penalty$value <- lambda model <- f(model, envir = envir.logLik) } else { - model <- f(model, envir = envir.logLik) + model <- f(model, envir = envir.logLik) } return(model) diff --git a/R/km1Nugget.init.R b/R/km1Nugget.init.R old mode 100755 new mode 100644 diff --git a/R/kmData.R b/R/kmData.R old mode 100755 new mode 100644 diff --git a/R/kmEstimate.R b/R/kmEstimate.R old mode 100755 new mode 100644 index a270b41..b98d504 --- a/R/kmEstimate.R +++ b/R/kmEstimate.R @@ -92,13 +92,13 @@ cat(" - variance bounds : ", c(lower[lp], upper[lp]), "\n") } cat(" - best initial criterion value(s) : ", initList$value, "\n") - if (model@optim.method=="BFGS") cat("\n") + if (model@optim.method=="bfgs") cat("\n") } # end printing # optimization - if (model@optim.method=="BFGS") { - BFGSargs <- c("trace", "parscale", "ndeps", "maxit", "abstol", "reltol", "REPORT", "lnm", "factr", "pgtol") + if (model@optim.method=="bfgs") { + BFGSargs <- c("trace", "parscale", "ndeps", "maxit", "abstol", "reltol", "REPORT", "lnm", "factr", "pgtol") commonNames <- intersect(BFGSargs, names(control)) controlChecked <- control[commonNames] if (length(control$REPORT)==0) { @@ -110,8 +110,7 @@ # multistart in parallel with foreach multistart <- control$multistart - - if (multistart==1){ + if (multistart==1) { model@parinit <- parinit <- as.numeric(parinit[, 1]) model@covariance <- initList$cov[[1]] o <- optim(par = parinit, fn = fn, gr = gr, @@ -120,16 +119,17 @@ model@control$convergence <- o$convergence } else { # multistart with foreach - if (requireNamespace("foreach", quietly = TRUE)){ - olist <- foreach::"%dopar%"(foreach::foreach(i=1:multistart, - .errorhandling='remove'), { - model@covariance <- initList$cov[[i]] - optim(par = parinit[, i], fn = fn, gr = gr, - method = "L-BFGS-B", lower = lower, upper = upper, - control = controlChecked, hessian = FALSE, model, envir=envir) + if (requireNamespace("foreach", quietly = TRUE)) { + olist <- foreach::"%dopar%"(foreach::foreach(i=seq(length=multistart), + .errorhandling='stop'), { + model@covariance <- initList$cov[[i]] + optim(par = parinit[, i], fn = fn, gr = gr, + method = "L-BFGS-B", lower = lower, upper = upper, + control = controlChecked, hessian = FALSE, model, envir=envir) }) + ## Note: if errorhandling drops the errors, the olist will be + ## shorter. I find it useful to have errors displayed. } - # get the best result bestValue <- Inf bestIndex <- NA @@ -142,8 +142,7 @@ bestValue <- currentValue bestIndex <- i } - } # end multistart - + } # end multistart model@covariance <- initList$cov[[bestIndex]] parinit <- parinit[, bestIndex] model@parinit <- as.numeric(parinit) @@ -195,11 +194,10 @@ lower = lower, upper = upper, control = control, model=model, envir=envir) } - model@logLik <- as.numeric(o$value) - if (model@method=="LOO"){ + if (model@method=="LOO") { model@covariance <- vect2covparam(model@covariance, o$par) diff --git a/R/kmNoNugget.init.R b/R/kmNoNugget.init.R old mode 100755 new mode 100644 diff --git a/R/kmNuggets.init.R b/R/kmNuggets.init.R old mode 100755 new mode 100644 diff --git a/R/kmStruct.R b/R/kmStruct.R old mode 100755 new mode 100644 index dce8f5c..fbebc43 --- a/R/kmStruct.R +++ b/R/kmStruct.R @@ -429,7 +429,7 @@ update.km <- function(object, if (is.null(kmcontrol$penalty)) kmcontrol$penalty <- object@penalty if (length(object@penalty == 0)) kmcontrol$penalty <- NULL if (is.null(kmcontrol$optim.method)) kmcontrol$optim.method <- object@optim.method - if (length(kmcontrol$optim.method) == 0) kmcontrol$optim.method <- "BFGS" + if (length(kmcontrol$optim.method) == 0) kmcontrol$optim.method <- "bfgs" if (is.null(kmcontrol$control)) kmcontrol$control <- object@control if(length(object@gr) == 0) object@gr <- TRUE knots <- NULL; if(TheClass == "covScaling") knots <- object@covariance@knots diff --git a/R/leaveOneOut.km.R b/R/leaveOneOut.km.R old mode 100755 new mode 100644 diff --git a/R/logLik.km.R b/R/logLik.km.R old mode 100755 new mode 100644 diff --git a/R/logLikFun.R b/R/logLikFun.R old mode 100755 new mode 100644 index b24ec2a..03068bf --- a/R/logLikFun.R +++ b/R/logLikFun.R @@ -1,12 +1,10 @@ `logLikFun` <- function(param, model, envir=NULL) { - if (model@known.param=="Trend") { beta <- model@trend.coef } else { beta <- NULL } - if (identical(model@case, "LLconcentration_beta_sigma2")) { model@covariance <- vect2covparam(model@covariance, param) @@ -31,13 +29,12 @@ function(param, model, envir=NULL) { } } else if (identical(model@case, "LLconcentration_beta")) { - nparam <- length(param) - if (class(model@covariance) != "covAdditive0") { + if (class(model@covariance) != "covAdditive0") { model@covariance <- vect2covparam(model@covariance, param[1:(nparam-1)]) model@covariance@sd2 <- param[nparam] } else { - model@covariance <- vect2covparam(model@covariance, param) + model@covariance <- vect2covparam(model@covariance, param) } # if (model@covariance@nugget.estim) { @@ -50,7 +47,6 @@ function(param, model, envir=NULL) { C <- aux[[1]] vn <- aux[[2]] - T <- chol(C) x <- backsolve(t(T), model@y, upper.tri = FALSE) M <- backsolve(t(T), model@F, upper.tri = FALSE) diff --git a/R/logLikGrad.R b/R/logLikGrad.R old mode 100755 new mode 100644 diff --git a/R/scalingFun.R b/R/scalingFun.R old mode 100755 new mode 100644 index 47601a7..1a15cc9 --- a/R/scalingFun.R +++ b/R/scalingFun.R @@ -111,9 +111,9 @@ scalingFun <- function(X, knots, eta, plot = FALSE) { transX <- matrix(NA, nrow = NROW(X), ncol = NCOL(X)) dimnames(transX) <- dimnames(X) X <- as.matrix(X) - for (i in 1:d){ - transX[ , i] <- scalingFun1d(x = X[ , i], knots = knots[[colnames(X)[i]]], eta = eta[[colnames(X)[i]]]) + transX[ , i] <- scalingFun1d(x = X[ , i], knots = knots[[i]], + eta = eta[[i]]) } if (plot) { diff --git a/R/scalingGrad.R b/R/scalingGrad.R old mode 100755 new mode 100644 diff --git a/R/show.km.R b/R/show.km.R old mode 100755 new mode 100644 diff --git a/R/trendMatrix.update.R b/R/trendMatrix.update.R old mode 100755 new mode 100644 diff --git a/man/km.Rd b/man/km.Rd index f8db796..275821f 100755 --- a/man/km.Rd +++ b/man/km.Rd @@ -23,15 +23,21 @@ km(formula=~1, design, response, covtype="matern5_2", \item{response}{ a vector (or 1-column matrix or data frame) containing the values of the 1-dimensional output given by the objective function at the \code{design} points. } \item{covtype}{ an optional character string specifying the covariance structure to be used, to be chosen between \code{"gauss"}, \code{"matern5_2"}, \code{"matern3_2"}, \code{"exp"} or \code{"powexp"}. See a full description of available covariance kernels in \code{\link{covTensorProduct-class}}. Default is \code{"matern5_2"}. See also the argument \code{kernel} that allows the user to build its own covariance structure.} \item{coef.trend,}{ (see below)} - \item{coef.cov,}{ (see below)} + \item{coef.cov}{ (see below)} \item{coef.var}{ optional vectors containing the values for the trend, covariance and variance parameters. For estimation, 4 cases are implemented: 1. (All unknown) If all are missing, all are estimated. 2. (All known) If all are provided, no estimation is performed; 3. (Known trend) If \code{coef.trend} is provided but at least one of \code{coef.cov} or \code{coef.var} is missing, then BOTH \code{coef.cov} and \code{coef.var} are estimated; 4. (Unknown trend) If \code{coef.cov} and \code{coef.var} are provided but \code{coef.trend} is missing, then \code{coef.trend} is estimated (GLS formula).} \item{nugget}{ an optional variance value standing for the homogeneous nugget effect.} \item{nugget.estim}{ an optional boolean indicating whether the nugget effect should be estimated. Note that this option does not concern the case of heterogeneous noisy observations (see \code{noise.var} below). If \code{nugget} is given, it is used as an initial value. Default is \code{FALSE}.} \item{noise.var}{ for noisy observations : an optional vector containing the noise variance at each observation. This is useful for stochastic simulators. Default is \code{NULL}.} \item{estim.method}{ a character string specifying the method by which unknown parameters are estimated. Default is \code{"MLE"} (Maximum Likelihood). At this stage, a beta version of leave-One-Out estimation (\code{estim.method="LOO"}) is also implemented for noise-free observations.} \item{penalty}{ (beta version) an optional list suitable for Penalized Maximum Likelihood Estimation. The list must contain the item \code{fun} indicating the penalty function, and the item \code{value} equal to the value of the penalty parameter. At this stage the only available \code{fun} is \code{"SCAD"}, and \code{covtype} must be \code{"gauss"}. Default is \code{NULL}, corresponding to (un-penalized) Maximum Likelihood Estimation.} - \item{optim.method}{ an optional character string indicating which optimization method is chosen for the likelihood maximization. \code{"BFGS"} is the \code{optim} quasi-Newton procedure of package \code{stats}, with the method "L-BFGS-B". \code{"gen"} is the \code{genoud} genetic algorithm (using derivatives) from package \code{rgenoud} (>= 5.3.3). } - \item{lower, }{ (see below) } + \item{optim.method}{ an optional character string indicating which + optimization method to choose for the likelihood + maximization. \code{"BFGS"} is the \code{optim} quasi-Newton + procedure of package \code{stats}, with the method + "L-BFGS-B". \code{"gen"} is the \code{genoud} genetic algorithm + (using derivatives) from package \code{rgenoud} (>= 5.3.3). Lower + case forms also supported.} + \item{lower}{ (see below) } \item{upper}{ optional vectors containing the bounds of the correlation parameters for optimization. The default values are given by \code{\link{covParametersBounds}}. } \item{parinit}{ an optional vector containing the initial values for the variables to be optimized over. If no vector is given, an initial point is generated as follows. For method \code{"gen"}, the initial point is generated uniformly inside the hyper-rectangle domain defined by \code{lower} and \code{upper}. For method \code{"BFGS"}, some points (see \code{control} below) are generated uniformly in the domain. Then the best point with respect to the likelihood (or penalized likelihood, see \code{penalty}) criterion is chosen. } \item{multistart}{ an optional integer indicating the number of initial points from which running the BFGS optimizer. These points will be selected as the best \code{multistart} one(s) among those evaluated (see above \code{parinit}). The multiple optimizations will be performed in parallel provided that a parallel backend is registered (see package \code{foreach}).} diff --git a/man/scalingFun.Rd b/man/scalingFun.Rd index 67c978b..595c2a8 100755 --- a/man/scalingFun.Rd +++ b/man/scalingFun.Rd @@ -14,8 +14,10 @@ scalingFun(X, knots, eta, plot=FALSE) \arguments{ \item{X}{ an n*d matrix standing for a design of n experiments in d-dimensional space } \item{knots}{ a list of knots parametrizing the transformation. } - \item{eta}{ a list of coefficients parametrizing the d marginal transformations. Each element stands for a set of marginal density values at the knots defined above.} - \item{plot}{ if TRUE plots the image of the columns of X according to the corresponding marginal transformations.} + \item{eta}{a list of coefficients parametrizing the d marginal + transformations. Each element stands for a set of marginal density + values at the knots defined above.} +\item{plot}{ if TRUE plots the image of the columns of X according to the corresponding marginal transformations.} } \value{ diff --git a/src/CovFuns.c b/src/CovFuns.c old mode 100755 new mode 100644 diff --git a/src/Scaling.c b/src/Scaling.c old mode 100755 new mode 100644