diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 049d60f..e9d8b55 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -38,7 +38,7 @@ jobs: - {os: macOS-latest, r: 'release'} - {os: windows-latest, r: 'devel', rtools: '44'} - {os: windows-latest, r: 'next', rtools: '44'} - - {os: windows-latest, r: 'release', rtools: '43'} + - {os: windows-latest, r: 'release', rtools: '44'} - {os: ubuntu-latest, r: 'devel'} - {os: ubuntu-latest, r: 'next'} - {os: ubuntu-latest, r: 'release'} diff --git a/DESCRIPTION b/DESCRIPTION index 89f6bcf..59f9b4a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,9 +2,9 @@ Package: gslnls Type: Package Title: GSL Multi-Start Nonlinear Least-Squares Fitting Version: 1.4.0 -Date: 2024-12-02 +Date: 2025-01-01 Authors@R: person("Joris", "Chau", email = "joris.chau@openanalytics.eu", role = c("aut", "cre")) -Description: An R interface to nonlinear least-squares optimization with the GNU Scientific Library (GSL), see M. Galassi et al. (2009, ISBN:0954612078). The available trust region methods include the Levenberg-Marquardt algorithm with and without geodesic acceleration, the Steihaug-Toint conjugate gradient algorithm for large systems and several variants of Powell's dogleg algorithm. Multi-start optimization based on quasi-random samples is implemented using a modified version of the algorithm in Hickernell and Yuan (1997, OR Transactions). Bindings are provided to tune a number of parameters affecting the low-level aspects of the trust region algorithms. The interface mimics R's nls() function and returns model objects inheriting from the same class. +Description: An R interface to nonlinear least-squares optimization with the GNU Scientific Library (GSL), see M. Galassi et al. (2009, ISBN:0954612078). The available trust region methods include the Levenberg-Marquardt algorithm with and without geodesic acceleration, the Steihaug-Toint conjugate gradient algorithm for large systems and several variants of Powell's dogleg algorithm. Multi-start optimization based on quasi-random samples is implemented using a modified version of the algorithm in Hickernell and Yuan (1997, OR Transactions). Robust nonlinear regression can be performed using various robust loss functions, in which case the optimization problem is solved by iterative reweighted least squares (IRLS). Bindings are provided to tune a number of parameters affecting the low-level aspects of the trust region algorithms. The interface mimics R's nls() function and returns model objects inheriting from the same class. BugReports: https://github.com/JorisChau/gslnls/issues URL: https://github.com/JorisChau/gslnls Depends: diff --git a/NAMESPACE b/NAMESPACE index c011f06..d239d9c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -42,6 +42,7 @@ importFrom(stats,formula) importFrom(stats,getInitial) importFrom(stats,hatvalues) importFrom(stats,model.weights) +importFrom(stats,naprint) importFrom(stats,nls) importFrom(stats,nls.control) importFrom(stats,nobs) diff --git a/NEWS.md b/NEWS.md index 4f2d9a2..dc53fd2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,10 @@ +# gslnls 1.4.0 + +* Robust loss optimization added in `gsl_nls()` via argument `loss` +* Added new function `gsl_nls_loss()` +* Added new method `cooks.distance()` +* Minor changes in `predict()` and `hatvalues()` for weighted NLS + # gslnls 1.3.3 * Fix standard errors `predict()` when using `newdata` diff --git a/R/nls.R b/R/nls.R index 5680a77..3a331d8 100644 --- a/R/nls.R +++ b/R/nls.R @@ -31,9 +31,6 @@ #' the values of \code{mstart_n}, \code{mstart_r} or \code{mstart_minsp} to avoid early termination of the algorithm at the cost of #' increased computational effort. #' -#' @section Robust loss functions: -#' If \code{loss} -#' #' @param fn a nonlinear model defined either as a two-sided \link{formula} including variables and parameters, #' or as a \link{function} returning a numeric vector, with first argument the vector of parameters to be estimated. #' See the individual method descriptions below. @@ -72,16 +69,18 @@ #' \itemize{ #' \item \code{"default"} default squared loss function. #' \item \code{"huber"} Huber loss function. -#' \item \code{"barron"} Smooth family of robust loss functions from (Barron (2019)). -#' \item \code{"biweight"} Tukey's bisquare loss function. -#' \item \code{"welsh"} Welsh loss function. -#' \item \code{"optimal"} Optimal loss function as given by (Maronna et al. (2006)). -#' \item \code{"hampel"} Hampel loss function (Hampel et al. (1986)). -#' \item \code{"ggw"} Generalized Gauss-Weight loss function (Koller and Stahel (2011)). -#' \item \code{"lqq"} Linear Quadratic Quadratic loss function (Koller and Stahel (2011)). -#' }. +#' \item \code{"barron"} Barron's smooth family of loss functions. +#' \item \code{"biweight"} Tukey's biweight/bisquare loss function. +#' \item \code{"welsh"} Welsh/Leclerc loss function. +#' \item \code{"optimal"} Optimal loss function (Maronna et al. (2006), Section 5.9.1). +#' \item \code{"hampel"} Hampel loss function. +#' \item \code{"ggw"} Generalized Gauss-Weight loss function. +#' \item \code{"lqq"} Linear Quadratic Quadratic loss function. +#' } #' If a character string, the default tuning parameters as specified by \code{\link{gsl_nls_loss}} are used. -#' Instead, a list as returned by \code{\link{gsl_nls_loss}} with non-default tuning parmaeters is also accepted. +#' Instead, a list as returned by \code{\link{gsl_nls_loss}} with non-default tuning parameters is also accepted. +#' For all choices other than \code{rho = "default"}, iterative reweighted least squares (IRLS) is used to +#' solve the MM-estimation problem. #' @param control an optional list of control parameters to tune the least squares iterations and multistart algorithm. #' See \code{\link{gsl_nls_control}} for the available control parameters and their default values. #' @param lower a named list or named numeric vector of parameter lower bounds, or an unnamed numeric @@ -126,7 +125,7 @@ #' applicable to objects of both classes. #' @useDynLib gslnls, .registration = TRUE #' @importFrom stats nls numericDeriv deriv as.formula coef deviance df.residual fitted vcov formula getInitial model.weights -#' @importFrom stats pf pt qt setNames sigma nobs hatvalues printCoefmat symnum cooks.distance +#' @importFrom stats pf pt qt setNames sigma nobs hatvalues printCoefmat symnum cooks.distance naprint #' @seealso \code{\link[stats]{nls}} #' @seealso \url{https://www.gnu.org/software/gsl/doc/html/nls.html} #' @seealso \url{https://CRAN.R-project.org/package=robustbase/vignettes/psi_functions.pdf} @@ -165,6 +164,14 @@ #' start = c(A = NA, lam = NA, b = NA) ## dynamic starting ranges #' ) #' +#' ## robust regression +#' gsl_nls( +#' fn = y ~ A * exp(-lam * x) + b, ## model formula +#' data = data.frame(x = x, y = y), ## model fit data +#' start = c(A = 0, lam = 0, b = 0), ## starting values +#' loss = "barron" ## L1-L2 loss +#' ) +#' #' ## analytic Jacobian 1 #' gsl_nls( #' fn = y ~ A * exp(-lam * x) + b, ## model formula @@ -303,7 +310,7 @@ gsl_nls <- function (fn, ...) { #' If \code{fn} is a \code{formula}, the returned list object is of classes \code{gsl_nls} and \code{nls}. #' Therefore, all generic functions applicable to objects of class \code{nls}, such as \code{anova}, \code{coef}, \code{confint}, #' \code{deviance}, \code{df.residual}, \code{fitted}, \code{formula}, \code{logLik}, \code{nobs}, \code{predict}, \code{print}, \code{profile}, -#' \code{residuals}, \code{summary}, \code{vcov} and \code{weights} are also applicable to the returned list object. +#' \code{residuals}, \code{summary}, \code{vcov}, \code{hatvalues}, \code{cooks.distance} and \code{weights} are also applicable to the returned list object. #' In addition, a method \code{confintd} is available for inference of derived parameters. #' @export gsl_nls.formula <- function(fn, data = parent.frame(), start, @@ -696,7 +703,7 @@ gsl_nls.formula <- function(fn, data = parent.frame(), start, cFit <- .Call(C_nls, .fn, .lhs, .jac, .fvv, environment(), start, wts, .lupars, .ctrl_int, .ctrl_dbl, .has_start, .loss_config, PACKAGE = "gslnls") ## convert to nls object - m <- nlsModel(formula, mf, cFit, wts, jac, upper = NULL) + m <- nlsModel(formula, mf, cFit, wts, jac) convInfo <- list( isConv = as.logical(!cFit$conv), @@ -748,8 +755,8 @@ gsl_nls.formula <- function(fn, data = parent.frame(), start, #' Although the returned object is not of class \code{nls}, the following generic functions remain #' applicable for an object of class \code{gsl_nls}: \code{anova}, \code{coef}, \code{confint}, \code{deviance}, #' \code{df.residual}, \code{fitted}, \code{formula}, \code{logLik}, \code{nobs}, \code{predict}, \code{print}, -#' \code{residuals}, \code{summary}, \code{vcov} and \code{weights}. In addition, a method \code{confintd} -#' is available for inference of derived parameters. +#' \code{residuals}, \code{summary}, \code{vcov}, \code{hatvalues}, \code{cooks.distance} and \code{weights}. +#' In addition, a method \code{confintd} is available for inference of derived parameters. #' @export gsl_nls.function <- function(fn, y, start, algorithm = c("lm", "lmaccel", "dogleg", "ddogleg", "subspace2D"), @@ -1097,7 +1104,7 @@ gsl_nls.function <- function(fn, y, start, #' @param mstart_maxiter positive integer, maximum number of iterations in the efficient local search algorithm (Algorithm B, Hickernell and Yuan (1997)), defaults to 10. #' @param mstart_maxstart positive integer, minimum number of major iterations (Algorithm 2.1, Hickernell and Yuan (1997)) before the multi-start algorithm terminates, defaults to 250. #' @param mstart_minsp positive integer, minimum number of detected stationary points before the multi-start algorithm terminates, defaults to 1. -#' @param irls_maxiter postivive integer, maximum number of IRLS iterations, defaults to 50. Only used in case of a non-default loss function (\code{loss != "default"}) optimized by IRLS. +#' @param irls_maxiter posivive integer, maximum number of IRLS iterations, defaults to 50. Only used in case of a non-default loss function (\code{loss != "default"}) optimized by IRLS. #' @param irls_xtol numeric value, termination of the IRLS procedure occurs when the relative change in parameters between IRLS iterations is \code{<= irls_xtol}, defaults to \code{.Machine$double.eps^(1/4)}. #' Only used in case of a non-default loss function (\code{loss != "default"}) optimized by IRLS. #' @param ... any additional arguments (currently not used). @@ -1378,18 +1385,24 @@ gslModel <- function(fn, lhs, cFit, start, wts, jac, ...) { gcall <- fcall } gcall[[2]] <- do.call(call, args = c(ifelse(is.list(start), "list", "c"), sapply(names(start), as.name)), quote = TRUE) - if(!any(is.na(cFit$grad) | is.infinite(cFit$grad))) { - QR <- qr(cFit$grad) + if(is.null(cFit$irls)) { + gr <- cFit$grad + } else { + .swts <- if(!missing(wts) && length(wts)) sqrt(wts) else rep_len(1, dim(cFit$grad)[1]) + gr <- .swts * cFit$grad / sqrt(cFit$irls$irls_weights) ## gradient w/o irls weights + } + if(!any(is.na(gr) | is.infinite(gr))) { + QR <- qr(gr) qrDim <- min(dim(QR$qr)) if(QR$rank < qrDim) warning("singular gradient matrix at parameter estimates") if(!is.null(cFit$irls)) { ## MM correction [Yohai, 1987. Thrm. 4.1] - tau <- mean(cFit$irlsirls_psi^2) / mean(cFit$irls$irls_dpsi)^2 + tau <- mean(cFit$irls$irls_psi^2) / mean(cFit$irls$irls_dpsi)^2 QR$qr <- QR$qr / sqrt(tau) } } else { - QR <- structure(list(qr = cFit$grad), class = "qr") + QR <- structure(list(qr = gr), class = "qr") QR$qr[] <- NA_real_ qrDim <- 0L warning("NA/NaN/Inf in gradient matrix at parameter estimates") @@ -1400,7 +1413,7 @@ gslModel <- function(fn, lhs, cFit, start, wts, jac, ...) { formula = function() fn, deviance = function() cFit$ssr, lhs = function() lhs, - gradient = function() cFit$grad, + gradient = function() gr, gradient1 = function(newdata = list()) { rho <- new.env(hash = TRUE, parent = env) for(i in names(newdata)) rho[[i]] <- newdata[[i]] diff --git a/R/nls_large.R b/R/nls_large.R index f7f2b1e..d08dbb0 100644 --- a/R/nls_large.R +++ b/R/nls_large.R @@ -127,8 +127,8 @@ gsl_nls_large <- function (fn, ...) { #' If \code{fn} is a \code{formula}, the returned list object is of classes \code{gsl_nls} and \code{nls}. #' Therefore, all generic functions applicable to objects of class \code{nls}, such as \code{anova}, \code{coef}, \code{confint}, #' \code{deviance}, \code{df.residual}, \code{fitted}, \code{formula}, \code{logLik}, \code{nobs}, \code{predict}, \code{print}, \code{profile}, -#' \code{residuals}, \code{summary}, \code{vcov} and \code{weights} are also applicable to the returned list object. -#' In addition, a method \code{confintd} is available for inference of derived parameters. +#' \code{residuals}, \code{summary}, \code{vcov}, \code{hatvalues}, \code{cooks.distance} and \code{weights} are also applicable to the returned +#' list object. In addition, a method \code{confintd} is available for inference of derived parameters. #' @export gsl_nls_large.formula <- function(fn, data = parent.frame(), start, algorithm = c("lm", "lmaccel", "dogleg", "ddogleg", "subspace2D", "cgst"), @@ -409,7 +409,7 @@ gsl_nls_large.formula <- function(fn, data = parent.frame(), start, cFit <- .Call(C_nls_large, .fn, .lhs, .jac, .fvv, environment(), start, wts, .ctrl_int, .ctrl_dbl, PACKAGE = "gslnls") ## convert to nls object - m <- nlsModel(formula, mf, cFit$par, wts, jac) + m <- nlsModel(formula, mf, cFit, wts, jac) convInfo <- list( isConv = as.logical(!cFit$conv), @@ -451,8 +451,8 @@ gsl_nls_large.formula <- function(fn, data = parent.frame(), start, #' Although the returned object is not of class \code{nls}, the following generic functions remain #' applicable for an object of class \code{gsl_nls}: \code{anova}, \code{coef}, \code{confint}, \code{deviance}, #' \code{df.residual}, \code{fitted}, \code{formula}, \code{logLik}, \code{nobs}, \code{predict}, \code{print}, -#' \code{residuals}, \code{summary}, \code{vcov} and \code{weights}. In addition, a method \code{confintd} -#' is available for inference of derived parameters. +#' \code{residuals}, \code{summary}, \code{vcov}, \code{hatvalues}, \code{cooks.distance} and \code{weights}. +#' In addition, a method \code{confintd} is available for inference of derived parameters. #' @export gsl_nls_large.function <- function(fn, y, start, algorithm = c("lm", "lmaccel", "dogleg", "ddogleg", "subspace2D", "cgst"), diff --git a/R/nls_methods.R b/R/nls_methods.R index e38e3ff..aa8f5ce 100644 --- a/R/nls_methods.R +++ b/R/nls_methods.R @@ -374,8 +374,7 @@ residuals.gsl_nls <- function(object, type = c("response", "pearson"), ...) { attr(val, "label") <- "Standardized residuals" } else { val <- as.vector(object$m$lhs() - object$m$fitted()) - lab <- "Residuals" - attr(val, "label") <- lab + attr(val, "label") <- "Residuals" } val } diff --git a/R/nls_rho.R b/R/nls_rho.R index fe92fa1..8ae8ff5 100644 --- a/R/nls_rho.R +++ b/R/nls_rho.R @@ -1,21 +1,107 @@ -#' Loss functions with tunable parameters +#' Robust loss functions with tunable parameters #' -#' Allow the user to tune the coefficient(s) of the loss functions -#' available in \code{\link{gsl_nls}} and \code{\link{gsl_nls_large}}. +#' @description +#' Allow the user to tune the coefficient(s) of the robust loss functions +#' supported by \code{\link{gsl_nls}}. For all choices other than \code{rho = "default"}, the MM-estimation +#' problem is optimized by means of iterative reweighted least squares (IRLS). #' -#' @param rho character loss function name, one of \code{"default", "huber", "barron", "absolute", "bisquare", "welsh", "optimal", "hampel", "ggw", "lqq"}. -#' @param cc numeric vector of loss function tuning parameters, with length depends on the chosen loss function. +#' @section Loss function details: +#' \describe{ +#' \item{\code{default}}{ +#' Default squared loss, no iterative reweighted least squares (IRLS) is required in this case. +#' \deqn{\rho(x) = x^2} +#' } +#' \item{\code{huber}}{ +#' Huber loss function with scaling constant \eqn{k}, defaulting to \eqn{k = 1.345} for 95\% efficiency of the regression estimator. +#' \deqn{\rho(x, k) = \left\{ \begin{array}{ll} \frac{1}{2} x^2 &\quad \text{ if } |x| \leq k \\[2pt] k(|x| - \frac{k}{2}) &\quad \text{ if } |x| > k \end{array} \right. } +#' } +#' \item{\code{barron}}{ +#' Barron's smooth family of loss functions with robustness parameter \eqn{\alpha \leq 2} (default \eqn{\alpha = 1}) and scaling constant \eqn{k} (default \eqn{k = 1.345}). +#' Special cases include: (scaled) squared loss for \eqn{\alpha = 2}; L1-L2 loss for \eqn{\alpha = 1}; Cauchy loss for \eqn{\alpha = 0}; Geman-McClure loss for \eqn{\alpha = -2}; +#' and Welsch/Leclerc loss for \eqn{\alpha = -\infty}. See Barron (2019) for additional details. +#' \deqn{\rho(x, \alpha, k) = \left\{ \begin{array}{ll} +#' \frac{1}{2} (x / k)^2 &\quad \text{ if } \alpha = 2 \\[2pt] +#' \log\left(\frac{1}{2} (x / k)^2 + 1 \right) &\quad \text{ if } \alpha = 0 \\[2pt] +#' 1 - \exp\left( -\frac{1}{2} (x / k)^2 \right) &\quad \text{ if } \alpha = -\infty \\[2pt] +#' \frac{|\alpha - 2|}{\alpha} \left( \left( \frac{(x / k)^2}{|\alpha-2|} + 1 \right)^{\alpha / 2} - 1 \right) &\quad \text{ otherwise } +#' \end{array} \right. +#' } +#' } +#' \item{\code{bisquare}}{ +#' Tukey's biweight/bisquare loss function with scaling constant \eqn{k}, defaulting to \eqn{k = 4.685} for 95\% efficiency of the regression estimator, +#' (\eqn{k = 1.548} gives a breakdown point of 0.5 of the S-estimator). +#' \deqn{\rho(x, k) = \left\{ \begin{array}{ll} 1 - (1 - (x / k)^2)^3 &\quad \text{ if } |x| \leq k \\[2pt] 1 &\quad \text{ if } |x| > k \end{array} \right. } +#' } +#' \item{\code{welsh}}{ +#' Welsh/Leclerc loss function with scaling constant \eqn{k}, defaulting to \eqn{k = 2.11} for 95\% efficiency of the regression estimator, (\eqn{k = 0.577} gives a +#' breakdown point of 0.5 of the S-estimator). This is equivalent to the Barron loss function with robustness parameter \eqn{\alpha = -\infty}. +#' \deqn{\rho(x, k) = 1 - \exp\left( -\frac{1}{2} (x / k)^2 \right) } +#' } +#' \item{\code{optimal}}{ +#' Optimal loss function as in Section 5.9.1 of Maronna et al. (2006) with scaling constant \eqn{k}, defaulting to \eqn{k = 1.060} for 95\% efficiency of the regression estimator, +#' (\eqn{k = 0.405} gives a breakdown point of 0.5 of the S-estimator). +#' \deqn{\rho(x, k) = \int_0^x \text{sign}(u) \left( - \dfrac{\phi'(|u|) + k}{\phi(|u|)} \right)^+ du} +#' with \eqn{\phi} the standard normal density. +#' } +#' \item{\code{hampel}}{ +#' Hampel loss function with a single scaling constant \eqn{k}, setting \eqn{a = 1.5k}, \eqn{b = 3.5k} and \eqn{r = 8k}. \eqn{k = 0.902} by default, +#' resulting in 95\% efficiency of the regression estimator, (\eqn{k = 0.212} gives a breakdown point of 0.5 of the S-estimator). +#' \deqn{\rho(x, a, b, r) = \left\{ \begin{array}{ll} +#' \frac{1}{2} x^2 / C &\quad \text{ if } |x| \leq a \\[2pt] +#' \left( \frac{1}{2}a^2 + a(|x|-a)\right) / C &\quad \text{ if } a < |x| \leq b \\[2pt] +#' \frac{a}{2}\left( 2b - a + (|x| - b) \left(1 + \frac{r - |x|}{r-b}\right) \right) / C &\quad \text{ if } b < |x| \leq r \\[2pt] +#' 1 &\quad \text{ if } r < |x| +#' \end{array} \right. } +#' with \eqn{C = \rho(\infty) = \frac{a}{2} (b - a + r)}. +#' } +#' \item{\code{ggw}}{ +#' Generalized Gauss-Weight loss function, a generalization of the Welsh/Leclerc loss with tuning constants \eqn{a, b, c}, defaulting to \eqn{a = 1.387}, +#' \eqn{b = 1.5}, \eqn{c = 1.063} for 95\% efficiency of the regression estimator, (\eqn{a = 0.204, b = 1.5, c = 0.296} gives a breakdown point of 0.5 of the S-estimator). +#' \deqn{\rho(x, a, b, c) = \int_0^x \psi(u, a, b, c) du} +#' with, +#' \deqn{\psi(x, a, b, c) = \left\{ \begin{array}{ll} +#' x &\quad \text{ if } |x| \leq c \\[2pt] +#' \exp\left( -\frac{1}{2} \frac{(|x| - c)^b}{a} \right) x &\quad \text{ if } |x| > c +#' \end{array} \right. } +#' } +#' \item{\code{lqq}}{ +#' Linear Quadratic Quadratic loss function with tuning constants \eqn{b, c, s}, defaulting to \eqn{b = 1.473}, \eqn{c = 0.982}, \eqn{s = 1.5} for +#' 95\% efficiency of the regression estimator, (\eqn{b = 0.402, c = 0.268, s = 1.5} gives a breakdown point of 0.5 of the S-estimator). +#' \deqn{\rho(x, b, c, s) = \int_0^x \psi(u, b, c, s) du} +#' with, +#' \deqn{\psi(x, b, c, s) = \left\{ \begin{array}{ll} +#' x &\quad \text{ if } |x| \leq c \\[2pt] +#' \text{sign}(x) \left( |x| - \frac{s}{2b} (|x| - c)^2 \right) &\quad \text{ if } c < |x| \leq b + c \\[2pt] +#' \text{sign}(x) \left( c + b - \frac{bs}{2} + \frac{s-1}{a} \left( \frac{1}{2} \tilde{x}^2 - a\tilde{x}\right) \right) &\quad \text{ if } b + c < |x| \leq a + b + c \\[2pt] +#' 0 &\quad \text{ otherwise } +#' \end{array} \right. } +#' where \eqn{\tilde{x} = |x| - b - c} and \eqn{a = (2c + 2b - bs) / (s - 1)}. +#' } +#' } +#' +#' @param rho character loss function, one of \code{"default", "huber", "barron", "bisquare", "welsh", "optimal", "hampel", "ggw", "lqq"}. +#' @param cc named or unnamed numeric vector of tuning parameters. The length of this argument depends on the selected \code{rho} function (see \sQuote{Details}). #' If \code{NULL}, the default tuning parameters are returned. -#' @seealso \url{https://CRAN.R-project.org/package=robustbase/vignettes/psi_functions.pdf} +#' @seealso \url{https://CRAN.R-project.org/package=robustbase} +#' @references J.T. Barron (2019). \emph{A general and adaptive robust loss function}. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition (pp. 4331-4339). +#' @references M. Galassi et al., \emph{GNU Scientific Library Reference Manual (3rd Ed.)}, ISBN 0954612078. +#' @references R.A. Maronna et al., \emph{Robust Statistics: Theory and Methods}, ISBN 0470010924. #' @examples -#' ## huber loss function with default tuning parameter -#' huber_loss() +#' ## Huber loss with default scale k +#' gsl_nls_loss(rho = "huber") +#' +#' ## L1-L2 loss (alpha = 1) +#' gsl_nls_loss(rho = "barron", cc = c(1, 1.345)) +#' +#' ## Cauchy loss (alpha = 0) +#' gsl_nls_loss(rho = "barron", cc = c(k = 1.345, alpha = 0)) +#' #' @return A \code{list} with two components: #' \itemize{ #' \item rho #' \item cc #' } -#' with meanings as explained under 'Arguments'. +#' with meanings as explained under \sQuote{Arguments}. #' @export gsl_nls_loss <- function(rho = c("default", "huber", "barron", "bisquare", "welsh", "optimal", "hampel", "ggw", "lqq"), cc = NULL) { @@ -23,29 +109,36 @@ gsl_nls_loss <- function(rho = c("default", "huber", "barron", "bisquare", "wels ## default tuning parameters cc_default <- switch(rho, - default = numeric(0), ## not used - huber = 1.345, - barron = c(1.0, 1.345), ## L1-L2 loss - bisquare = 4.685061, - welsh = 2.11, - optimal = 1.060158, - hampel = c(1.352413, 3.155630, 7.212868), - ggw = c(0.648, 1.0, 1.694), - lqq = c(1.473, 0.982, 1.5) + default = numeric(0), + huber = c(k = 1.345), + barron = c(alpha = 1.0, k = 1.345), + bisquare = c(k = 4.685061), + welsh = c(k = 2.11), + optimal = c(k = 1.060158), + hampel = c(k = 0.9016085), + ggw = c(a = 1.387, b = 1.5, c = 1.063), + lqq = c(b = 1.473, c = 0.982, s = 1.5) ) + ## check tuning parameters if(is.null(cc) || identical(rho, "default")) { cc <- cc_default } else if(length(cc) != length(cc_default)) { stop(sprintf("'cc' must be of length %d for function '%s'", length(cc_default), rho)) - } else { + } else if(is.null(names(cc))) { cc <- as.numeric(cc) - if(identical(rho, "barron") && cc[1] > 2) { - warning("Robustness parameter (alpha) in Barron loss function cannot be larger than 2") - cc[1] <- 2 + names(cc) <- names(cc_default) + } else { + if(!all(is.element(names(cc_default), names(cc)))) { + stop(sprintf("'cc' must be unnamed or include names %s", paste(names(cc_default), collapse = ", "))) } + cc <- as.numeric(cc)[match(names(cc_default), names(cc))] + names(cc) <- names(cc_default) + } + if(identical(rho, "barron") && cc[["alpha"]] > 2) { + warning("Robustness parameter (alpha) in Barron loss function cannot be larger than 2") + cc[["alpha"]] <- 2 } - list(rho = rho, cc = cc) - + return(list(rho = rho, cc = cc)) } diff --git a/README.md b/README.md index 2eb4ab8..ceadfa8 100644 --- a/README.md +++ b/README.md @@ -18,11 +18,12 @@ optimization with the [GNU Scientific Library (GSL)](https://www.gnu.org/software/gsl/). The function `gsl_nls()` solves small to moderate sized nonlinear least-squares problems with the `gsl_multifit_nlinear` interface with built-in support for multi-start -optimization. For large problems, where factoring the full Jacobian -matrix becomes prohibitively expensive, the `gsl_nls_large()` function -can be used to solve the system with the `gsl_multilarge_nlinear` -interface. The `gsl_nls_large()` function is also appropriate for -systems with sparse structure in the Jacobian matrix. +optimization and nonlinear robust regression. For large problems, where +factoring the full Jacobian matrix becomes prohibitively expensive, the +`gsl_nls_large()` function can be used to solve the system with the +`gsl_multilarge_nlinear` interface. The `gsl_nls_large()` function is +also appropriate for systems with sparse structure in the Jacobian +matrix. The following trust region methods to solve nonlinear least-squares problems are available in `gsl_nls()` (and `gsl_nls_large()`): @@ -81,8 +82,9 @@ instructions in the included README and INSTALL files. ##### Windows -A binary version of GSL (2.7) can be installed using the Rtools package -manager (see +Using Rtools \>= 42, GSL is available with the Rtools installation. For +earlier versions of Rtools, a binary version of GSL can be installed +using the Rtools package manager (see e.g. ): pacman -S mingw-w64-{i686,x86_64}-gsl @@ -126,7 +128,7 @@ install.packages("gslnls") #### Data -The code below simulates +Below, we simulate ![n = 25](https://latex.codecogs.com/png.latex?n%20%3D%2025 "n = 25") noisy observations ![y_1,\ldots,y_n](https://latex.codecogs.com/png.latex?y_1%2C%5Cldots%2Cy_n "y_1,\ldots,y_n") @@ -149,7 +151,7 @@ The exponential model parameters are set to ![A = 5](https://latex.codecogs.com/png.latex?A%20%3D%205 "A = 5"), ![\lambda = 1.5](https://latex.codecogs.com/png.latex?%5Clambda%20%3D%201.5 "\lambda = 1.5"), ![b = 1](https://latex.codecogs.com/png.latex?b%20%3D%201 "b = 1"), with -a noise standard deviation of +noise standard deviation ![\sigma = 0.25](https://latex.codecogs.com/png.latex?%5Csigma%20%3D%200.25 "\sigma = 0.25"). ``` r @@ -164,11 +166,11 @@ y <- f(A = 5, lam = 1.5, b = 1, x) + rnorm(n, sd = 0.25) #### Model fit -The exponential model is fitted to the data using the `gsl_nls()` -function by passing the nonlinear model as a two-sided `formula` and +The exponential model is fitted to the data using the function +`gsl_nls()` by passing the nonlinear model as a two-sided `formula` and providing starting parameters for the model parameters ![A, \lambda, b](https://latex.codecogs.com/png.latex?A%2C%20%5Clambda%2C%20b "A, \lambda, b") -analogous to an `nls()` function call. +analogous to a standard `nls()` call. ``` r library(gslnls) @@ -182,7 +184,6 @@ ex1_fit <- gsl_nls( ex1_fit #> Nonlinear regression model #> model: y ~ A * exp(-lam * x) + b -#> data: data.frame(x = x, y = y) #> A lam b #> 4.893 1.417 1.010 #> residual sum-of-squares: 1.316 @@ -190,17 +191,17 @@ ex1_fit #> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr) #> #> Number of iterations to convergence: 9 -#> Achieved convergence tolerance: 0 +#> Achieved convergence tolerance: 4.441e-16 ``` -Here, the nonlinear least squares problem has been solved with the +Here, the nonlinear least squares problem is solved with the Levenberg-Marquardt algorithm (default) in the `gsl_multifit_nlinear` -interface using the following `control` parameters: +interface with `control` parameters: ``` r ## default control parameters gsl_nls_control() |> str() -#> List of 21 +#> List of 23 #> $ maxiter : int 100 #> $ scale : chr "more" #> $ solver : chr "qr" @@ -222,12 +223,14 @@ gsl_nls_control() |> str() #> $ mstart_maxiter : int 10 #> $ mstart_maxstart: int 250 #> $ mstart_minsp : int 1 +#> $ irls_maxiter : int 50 +#> $ irls_xtol : num 0.000122 ``` -Run `?gsl_nls_control` or check the [GSL reference +Check `?gsl_nls_control` or the [GSL reference manual](https://www.gnu.org/software/gsl/doc/html/nls.html#tunable-parameters) -for further details on the available tuning parameters to control the -trust region algorithms. +for details on the available tuning parameters to control the trust +region algorithms. #### Object methods @@ -235,8 +238,8 @@ The fitted model object returned by `gsl_nls()` is of class `"gsl_nls"`, which inherits from class `"nls"`. For this reason, generic functions such as `anova`, `coef`, `confint`, `deviance`, `df.residual`, `fitted`, `formula`, `logLik`, `predict`, `print` `profile`, `residuals`, -`summary`, `vcov` and `weights` are also applicable for models fitted -with `gsl_nls()`. +`summary`, `vcov`, `hatvalues`, `cooks.distance` and `weights` are also +applicable for models fitted with `gsl_nls()`. ``` r ## model summary @@ -255,7 +258,7 @@ summary(ex1_fit) #> Residual standard error: 0.2446 on 22 degrees of freedom #> #> Number of iterations to convergence: 9 -#> Achieved convergence tolerance: 0 +#> Achieved convergence tolerance: 4.441e-16 ## asymptotic confidence intervals confint(ex1_fit) @@ -300,26 +303,25 @@ confintd(ex1_fit, expr = c("b", "A + b", "log(lam)"), level = 0.95) #### Jacobian calculation -If the `jac` argument in `gsl_nls()` is undefined, the Jacobian matrix -used to solve the [trust region +If the `jac` argument in `gsl_nls()` is undefined, the Jacobian to solve +the [trust region subproblem](https://www.gnu.org/software/gsl/doc/html/nls.html#solving-the-trust-region-subproblem-trs) -is approximated by forward (or centered) finite differences. Instead, an -analytic Jacobian can be passed to `jac` by defining a function that -returns the +is approximated numerically by (forward or centered) finite differences. +Instead, an analytic Jacobian can be passed to `jac` by defining a +function that returns the ![(n \times p)](https://latex.codecogs.com/png.latex?%28n%20%5Ctimes%20p%29 "(n \times p)")-dimensional Jacobian matrix of the nonlinear model `fn`, where the first argument must be the vector of parameters of length ![p](https://latex.codecogs.com/png.latex?p "p"). -In the exponential model example, the Jacobian matrix is a -![(50 \times 3)](https://latex.codecogs.com/png.latex?%2850%20%5Ctimes%203%29 "(50 \times 3)")-dimensional -matrix +In the above example, the Jacobian is a +![(50 \times 3)](https://latex.codecogs.com/png.latex?%2850%20%5Ctimes%203%29 "(50 \times 3)")-matrix ![\[\boldsymbol{J}\_{ij}\]\_{ij}](https://latex.codecogs.com/png.latex?%5B%5Cboldsymbol%7BJ%7D_%7Bij%7D%5D_%7Bij%7D "[\boldsymbol{J}_{ij}]_{ij}") with rows: ![\boldsymbol{J}\_i \\= \\\left\[ \frac{\partial f_i}{\partial A}, \frac{\partial f_i}{\partial \lambda}, \frac{\partial f_i}{\partial b} \right\] \\= \\\left\[ \exp(-\lambda \cdot x_i), -A \cdot \exp(-\lambda \cdot x_i) \cdot x_i, 1 \right\]](https://latex.codecogs.com/png.latex?%5Cboldsymbol%7BJ%7D_i%20%5C%20%3D%20%5C%20%5Cleft%5B%20%5Cfrac%7B%5Cpartial%20f_i%7D%7B%5Cpartial%20A%7D%2C%20%5Cfrac%7B%5Cpartial%20f_i%7D%7B%5Cpartial%20%5Clambda%7D%2C%20%5Cfrac%7B%5Cpartial%20f_i%7D%7B%5Cpartial%20b%7D%20%5Cright%5D%20%5C%20%3D%20%5C%20%5Cleft%5B%20%5Cexp%28-%5Clambda%20%5Ccdot%20x_i%29%2C%20-A%20%5Ccdot%20%5Cexp%28-%5Clambda%20%5Ccdot%20x_i%29%20%5Ccdot%20x_i%2C%201%20%5Cright%5D "\boldsymbol{J}_i \ = \ \left[ \frac{\partial f_i}{\partial A}, \frac{\partial f_i}{\partial \lambda}, \frac{\partial f_i}{\partial b} \right] \ = \ \left[ \exp(-\lambda \cdot x_i), -A \cdot \exp(-\lambda \cdot x_i) \cdot x_i, 1 \right]") -which is encoded in the following call to `gsl_nls()`: +which can be encoded as: ``` r ## analytic Jacobian (1) @@ -332,7 +334,6 @@ gsl_nls( ) #> Nonlinear regression model #> model: y ~ A * exp(-lam * x) + b -#> data: data.frame(x = x, y = y) #> A lam b #> 4.893 1.417 1.010 #> residual sum-of-squares: 1.316 @@ -345,8 +346,7 @@ gsl_nls( If the model formula `fn` can be derived with `stats::deriv()`, then the analytic Jacobian in `jac` can be computed automatically using symbolic -differentiation and no manual calculations are necessary. To evaluate -`jac` by means of symbolic differentiation, set `jac = TRUE`: +differentiation and it suffices to set `jac = TRUE`: ``` r ## analytic Jacobian (2) @@ -358,7 +358,6 @@ gsl_nls( ) #> Nonlinear regression model #> model: y ~ A * exp(-lam * x) + b -#> data: data.frame(x = x, y = y) #> A lam b #> 4.893 1.417 1.010 #> residual sum-of-squares: 1.316 @@ -383,7 +382,6 @@ ss_fit <- gsl_nls( ss_fit #> Nonlinear regression model #> model: y ~ SSasymp(x, Asym, R0, lrc) -#> data: data.frame(x = x, y = y) #> Asym R0 lrc #> 1.0097 5.9028 0.3484 #> residual sum-of-squares: 1.316 @@ -391,7 +389,7 @@ ss_fit #> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr) #> #> Number of iterations to convergence: 1 -#> Achieved convergence tolerance: 1.068e-13 +#> Achieved convergence tolerance: 1.181e-13 ``` The self-starting model `SSasymp()` uses a different model @@ -427,15 +425,14 @@ gsl_nls( ) #> Nonlinear regression model #> model: y ~ A * exp(-lam * x) + b -#> data: data.frame(x = x, y = y) #> A lam b #> 4.893 1.417 1.010 #> residual sum-of-squares: 1.316 #> #> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr) #> -#> Number of iterations to convergence: 2 -#> Achieved convergence tolerance: 6.106e-14 +#> Number of iterations to convergence: 3 +#> Achieved convergence tolerance: 8.882e-16 ``` The multi-start procedure is a modified version of the algorithm @@ -454,7 +451,6 @@ gsl_nls( ) #> Nonlinear regression model #> model: y ~ A * exp(-lam * x) + b -#> data: data.frame(x = x, y = y) #> A lam b #> 4.893 1.417 1.010 #> residual sum-of-squares: 1.316 @@ -462,7 +458,7 @@ gsl_nls( #> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr) #> #> Number of iterations to convergence: 3 -#> Achieved convergence tolerance: 0 +#> Achieved convergence tolerance: 2.22e-16 ``` **Remark**: the dynamic multi-start procedure is not expected to always @@ -471,6 +467,44 @@ objective function contains many local optima, the multi-start algorithm may be unable to select parameter ranges that include the global minimizing solution. +#### Robust loss functions + +The `loss` argument in `gsl_nls()` allows to specify a robust loss +function +![\rho(x)](https://latex.codecogs.com/png.latex?%5Crho%28x%29 "\rho(x)") +other than the default squared loss +![\rho(x) = \frac{1}{2}x^2](https://latex.codecogs.com/png.latex?%5Crho%28x%29%20%3D%20%5Cfrac%7B1%7D%7B2%7Dx%5E2 "\rho(x) = \frac{1}{2}x^2") +in the optimization objective: + +![\arg \min\_{\boldsymbol{\theta}} \sum\_{i=1}^n \rho\left(y_i - f(x_i, \boldsymbol{\theta})\right)](https://latex.codecogs.com/png.latex?%5Carg%20%5Cmin_%7B%5Cboldsymbol%7B%5Ctheta%7D%7D%20%5Csum_%7Bi%3D1%7D%5En%20%5Crho%5Cleft%28y_i%20-%20f%28x_i%2C%20%5Cboldsymbol%7B%5Ctheta%7D%29%5Cright%29 "\arg \min_{\boldsymbol{\theta}} \sum_{i=1}^n \rho\left(y_i - f(x_i, \boldsymbol{\theta})\right)") + +The (MM-)estimates in the robust regression problem are obtained by +iterative reweighted least squares (IRLS), see `?gsl_nls_loss` for the +available loss functions and their (optional) tuning parameters. + +``` r +## Huber loss +gsl_nls( + fn = y ~ A * exp(-lam * x) + b, ## model formula + data = data.frame(x = x, y = y), ## model fit data + start = c(A = 0, lam = 0, b = 0), ## starting values + loss = "huber" +) +#> Nonlinear regression model +#> model: y ~ A * exp(-lam * x) + b +#> A lam b +#> 4.796 1.463 1.092 +#> weighted residual sum-of-squares: 0.8127 +#> +#> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr) +#> +#> Number of IRLS iterations to convergence: 8 +#> Achieved IRLS tolerance: 0.0001023 +#> +#> Number of NLS iterations to convergence: 9 +#> Achieved NLS tolerance: 1.11e-16 +``` + ### Example 2: Gaussian function #### Data @@ -535,23 +569,23 @@ ex2a_fit <- gsl_nls( #> iter 2: ssr = 171.29, par = (1.85555, 2.84851, -5.66642) #> iter 3: ssr = 168.703, par = (1.93316, 2.34745, -6.33662) #> iter 4: ssr = 167.546, par = (1.84665, 1.24131, -7.399) -#> iter 5: ssr = 166.799, par = (1.87948, 0.380488, -7.61723) -#> iter 6: ssr = 165.623, par = (1.90658, -2.00515, -7.19306) -#> iter 7: ssr = 163.944, par = (2.18675, -3.19622, -5.38804) +#> iter 5: ssr = 166.799, par = (1.87948, 0.380487, -7.61723) +#> iter 6: ssr = 165.623, par = (1.90658, -2.00516, -7.19306) +#> iter 7: ssr = 163.944, par = (2.18675, -3.19622, -5.38805) #> iter 8: ssr = 161.679, par = (2.41824, -3.1487, -5.03149) -#> iter 9: ssr = 159.955, par = (2.67372, -3.18703, -4.20169) -#> iter 10: ssr = 157.391, par = (3.1083, -2.84066, -3.10574) -#> iter 11: ssr = 153.131, par = (3.54614, -2.38071, -2.53824) -#> iter 12: ssr = 150.129, par = (4.02799, -1.97822, -1.96376) -#> iter 13: ssr = 146.735, par = (4.49887, -1.35879, -1.39554) -#> iter 14: ssr = 142.928, par = (3.14455, -0.573017, -1.08117) -#> iter 15: ssr = 124.562, par = (2.14743, 0.465351, -0.443335) -#> iter 16: ssr = 102.183, par = (3.74451, 0.263346, 0.353849) -#> iter 17: ssr = 35.4869, par = (3.49962, 0.437339, 0.18613) -#> iter 18: ssr = 9.24985, par = (4.41273, 0.381446, 0.16023) -#> iter 19: ssr = 3.13669, par = (4.94382, 0.40041, 0.150317) -#> iter 20: ssr = 2.7623, par = (5.11741, 0.397815, 0.14729) -#> iter 21: ssr = 2.75831, par = (5.13786, 0.397886, 0.146865) +#> iter 9: ssr = 159.955, par = (2.67373, -3.18704, -4.20168) +#> iter 10: ssr = 157.391, par = (3.10831, -2.84066, -3.10573) +#> iter 11: ssr = 153.131, par = (3.54614, -2.3807, -2.53824) +#> iter 12: ssr = 150.129, par = (4.028, -1.97822, -1.96376) +#> iter 13: ssr = 146.735, par = (4.49887, -1.35878, -1.39554) +#> iter 14: ssr = 142.928, par = (3.14454, -0.573011, -1.08117) +#> iter 15: ssr = 124.564, par = (2.14741, 0.465371, -0.443323) +#> iter 16: ssr = 102.171, par = (3.74447, 0.263284, 0.353778) +#> iter 17: ssr = 35.4992, par = (3.4995, 0.437365, 0.186149) +#> iter 18: ssr = 9.25455, par = (4.41259, 0.381433, 0.160233) +#> iter 19: ssr = 3.13721, par = (4.94369, 0.400412, 0.15032) +#> iter 20: ssr = 2.76231, par = (5.11738, 0.397815, 0.14729) +#> iter 21: ssr = 2.75831, par = (5.13785, 0.397886, 0.146865) #> iter 22: ssr = 2.7583, par = (5.1389, 0.397884, 0.146831) #> iter 23: ssr = 2.7583, par = (5.13894, 0.397884, 0.146829) #> iter 24: ssr = 2.7583, par = (5.13894, 0.397884, 0.146829) @@ -560,21 +594,19 @@ ex2a_fit <- gsl_nls( #> ******************* #> summary from method 'multifit/levenberg-marquardt' #> number of iterations: 26 -#> reason for stopping: input domain error -#> initial ssr = 210.146 -#> final ssr = 2.7583 -#> ssr/dof = 0.0586872 -#> ssr achieved tolerance = 8.88178e-16 -#> function evaluations: 125 +#> initial ssr: 210.146 +#> final ssr: 2.7583 +#> ssr/dof: 0.0586872 +#> ssr achieved tolerance: 1.33227e-15 +#> function evaluations: 124 #> jacobian evaluations: 0 #> fvv evaluations: 0 -#> status = success +#> status: success #> ******************* ex2a_fit #> Nonlinear regression model #> model: y ~ a * exp(-(x - b)^2/(2 * c^2)) -#> data: data.frame(x = x, y = y) #> a b c #> 5.1389 0.3979 0.1468 #> residual sum-of-squares: 2.758 @@ -582,7 +614,7 @@ ex2a_fit #> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr) #> #> Number of iterations to convergence: 26 -#> Achieved convergence tolerance: 8.882e-16 +#> Achieved convergence tolerance: 1.332e-15 ``` #### Geodesic acceleration @@ -604,8 +636,8 @@ ex2b_fit <- gsl_nls( #> iter 1: ssr = 158.039, par = (1.58476, 0.502555, 0.511498) #> iter 2: ssr = 126.469, par = (1.8444, 0.366374, 0.403898) #> iter 3: ssr = 77.0115, par = (2.5025, 0.392374, 0.272717) -#> iter 4: ssr = 26.4036, par = (3.64063, 0.394564, 0.205619) -#> iter 5: ssr = 5.35492, par = (4.63506, 0.396865, 0.16366) +#> iter 4: ssr = 26.4036, par = (3.64063, 0.394564, 0.205618) +#> iter 5: ssr = 5.35491, par = (4.63506, 0.396865, 0.163659) #> iter 6: ssr = 2.82529, par = (5.05578, 0.397909, 0.149333) #> iter 7: ssr = 2.75877, par = (5.13246, 0.397896, 0.147057) #> iter 8: ssr = 2.7583, par = (5.13862, 0.397885, 0.146843) @@ -616,21 +648,19 @@ ex2b_fit <- gsl_nls( #> ******************* #> summary from method 'multifit/levenberg-marquardt+accel' #> number of iterations: 12 -#> reason for stopping: input domain error -#> initial ssr = 210.146 -#> final ssr = 2.7583 -#> ssr/dof = 0.0586872 -#> ssr achieved tolerance = 3.24185e-14 +#> initial ssr: 210.146 +#> final ssr: 2.7583 +#> ssr/dof: 0.0586872 +#> ssr achieved tolerance: 3.19744e-14 #> function evaluations: 76 #> jacobian evaluations: 0 #> fvv evaluations: 0 -#> status = success +#> status: success #> ******************* ex2b_fit #> Nonlinear regression model #> model: y ~ a * exp(-(x - b)^2/(2 * c^2)) -#> data: data.frame(x = x, y = y) #> a b c #> 5.1389 0.3979 0.1468 #> residual sum-of-squares: 2.758 @@ -638,7 +668,7 @@ ex2b_fit #> Algorithm: multifit/levenberg-marquardt+accel, (scaling: more, solver: qr) #> #> Number of iterations to convergence: 12 -#> Achieved convergence tolerance: 3.242e-14 +#> Achieved convergence tolerance: 3.197e-14 ``` With geodesic acceleration enabled the method converges after 12 @@ -754,19 +784,17 @@ gsl_nls( #> ******************* #> summary from method 'multifit/levenberg-marquardt+accel' #> number of iterations: 12 -#> reason for stopping: input domain error -#> initial ssr = 210.146 -#> final ssr = 2.7583 -#> ssr/dof = 0.0586872 -#> ssr achieved tolerance = 3.28626e-14 +#> initial ssr: 210.146 +#> final ssr: 2.7583 +#> ssr/dof: 0.0586872 +#> ssr achieved tolerance: 3.10862e-14 #> function evaluations: 58 #> jacobian evaluations: 0 #> fvv evaluations: 18 -#> status = success +#> status: success #> ******************* #> Nonlinear regression model #> model: y ~ a * exp(-(x - b)^2/(2 * c^2)) -#> data: data.frame(x = x, y = y) #> a b c #> 5.1389 0.3979 0.1468 #> residual sum-of-squares: 2.758 @@ -774,7 +802,7 @@ gsl_nls( #> Algorithm: multifit/levenberg-marquardt+accel, (scaling: more, solver: qr) #> #> Number of iterations to convergence: 12 -#> Achieved convergence tolerance: 3.286e-14 +#> Achieved convergence tolerance: 3.109e-14 ``` If the model formula `fn` can be derived with `stats::deriv()`, then the @@ -808,19 +836,17 @@ gsl_nls( #> ******************* #> summary from method 'multifit/levenberg-marquardt+accel' #> number of iterations: 12 -#> reason for stopping: input domain error -#> initial ssr = 210.146 -#> final ssr = 2.7583 -#> ssr/dof = 0.0586872 -#> ssr achieved tolerance = 2.93099e-14 +#> initial ssr: 210.146 +#> final ssr: 2.7583 +#> ssr/dof: 0.0586872 +#> ssr achieved tolerance: 3.10862e-14 #> function evaluations: 58 #> jacobian evaluations: 0 #> fvv evaluations: 18 -#> status = success +#> status: success #> ******************* #> Nonlinear regression model #> model: y ~ a * exp(-(x - b)^2/(2 * c^2)) -#> data: data.frame(x = x, y = y) #> a b c #> 5.1389 0.3979 0.1468 #> residual sum-of-squares: 2.758 @@ -828,7 +854,7 @@ gsl_nls( #> Algorithm: multifit/levenberg-marquardt+accel, (scaling: more, solver: qr) #> #> Number of iterations to convergence: 12 -#> Achieved convergence tolerance: 2.931e-14 +#> Achieved convergence tolerance: 3.109e-14 ``` ### Example 3: Branin function @@ -914,7 +940,7 @@ ex3_fit #> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr) #> #> Number of iterations to convergence: 20 -#> Achieved convergence tolerance: 0 +#> Achieved convergence tolerance: 1.665e-16 ``` **Note**: When using the `function` method of `gsl_nls()`, the returned @@ -1047,7 +1073,7 @@ system.time({ ) }) #> user system elapsed -#> 33.720 0.084 33.812 +#> 57.948 0.229 58.200 cat("Residual sum-of-squares:", deviance(ex4_fit_lm), "\n") #> Residual sum-of-squares: 0.004778845 @@ -1069,7 +1095,7 @@ system.time({ ) }) #> user system elapsed -#> 1.214 0.312 1.527 +#> 2.154 0.664 2.819 cat("Residual sum-of-squares:", deviance(ex4_fit_cgst), "\n") #> Residual sum-of-squares: 0.004778845 @@ -1114,10 +1140,10 @@ bench::mark( #> # A tibble: 4 × 6 #> expression min median `itr/sec` mem_alloc `gc/sec` #> -#> 1 Dense LM 6.33s 6.39s 0.157 1.8GB 6.61 -#> 2 Dense CGST 1.29s 1.36s 0.726 1.02GB 16.4 -#> 3 Sparse LM 4.87s 4.91s 0.203 33.19MB 0.122 -#> 4 Sparse CGST 166.09ms 176.09ms 4.66 23.04MB 3.73 +#> 1 Dense LM 11.91s 12.36s 0.0816 1.82GB 3.39 +#> 2 Dense CGST 1.95s 2.03s 0.467 1.02GB 9.89 +#> 3 Sparse LM 9.61s 9.62s 0.104 33.41MB 0.0624 +#> 4 Sparse CGST 252.96ms 262.35ms 3.83 23.04MB 2.30 ``` ## NLS test problems @@ -1215,7 +1241,7 @@ solutions and a vector of suggested starting values for the parameters: #> #> $fn #> y ~ b1/(1 + exp(b2 - b3 * x)) -#> +#> #> #> $start #> b1 b2 b3 @@ -1245,7 +1271,7 @@ with(ratkowsky2, #> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr) #> #> Number of iterations to convergence: 10 -#> Achieved convergence tolerance: 4.619e-14 +#> Achieved convergence tolerance: 5.151e-14 ## example optimization problem madsen <- nls_test_problem(name = "Madsen") @@ -1266,7 +1292,7 @@ with(madsen, #> Algorithm: multifit/levenberg-marquardt, (scaling: more, solver: qr) #> #> Number of iterations to convergence: 42 -#> Achieved convergence tolerance: 1.11e-16 +#> Achieved convergence tolerance: 2.22e-16 ``` ## Other R-packages @@ -1281,9 +1307,19 @@ for this package include: - [RcppZiggurat](https://cran.r-project.org/web/packages/RcppZiggurat/index.html) by Dirk Eddelbuettel +Except for Barron’s family of loss functions, all robust loss functions +are credited to +[robustbase](https://cran.r-project.org/web/packages/robustbase/index.html) +by Martin Maechler and others. + +## License + +LGPL-3 + # References -
+
diff --git a/configure b/configure index 274e0c7..e062daa 100755 --- a/configure +++ b/configure @@ -1,6 +1,6 @@ #! /bin/sh # Guess values for system-dependent variables and create Makefiles. -# Generated by GNU Autoconf 2.71 for gslnls 1.3.2. +# Generated by GNU Autoconf 2.71 for gslnls 1.4.0. # # Report bugs to . # @@ -610,8 +610,8 @@ MAKEFLAGS= # Identity of this package. PACKAGE_NAME='gslnls' PACKAGE_TARNAME='gslnls' -PACKAGE_VERSION='1.3.2' -PACKAGE_STRING='gslnls 1.3.2' +PACKAGE_VERSION='1.4.0' +PACKAGE_STRING='gslnls 1.4.0' PACKAGE_BUGREPORT='https://github.com/JorisChau/gslnls/issues/' PACKAGE_URL='' @@ -1231,7 +1231,7 @@ if test "$ac_init_help" = "long"; then # Omit some internal or obsolete options to make the list less imposing. # This message is too long to be a string in the A/UX 3.1 sh. cat <<_ACEOF -\`configure' configures gslnls 1.3.2 to adapt to many kinds of systems. +\`configure' configures gslnls 1.4.0 to adapt to many kinds of systems. Usage: $0 [OPTION]... [VAR=VALUE]... @@ -1293,7 +1293,7 @@ fi if test -n "$ac_init_help"; then case $ac_init_help in - short | recursive ) echo "Configuration of gslnls 1.3.2:";; + short | recursive ) echo "Configuration of gslnls 1.4.0:";; esac cat <<\_ACEOF @@ -1375,7 +1375,7 @@ fi test -n "$ac_init_help" && exit $ac_status if $ac_init_version; then cat <<\_ACEOF -gslnls configure 1.3.2 +gslnls configure 1.4.0 generated by GNU Autoconf 2.71 Copyright (C) 2021 Free Software Foundation, Inc. @@ -1533,7 +1533,7 @@ cat >config.log <<_ACEOF This file contains any messages produced by compilers while running configure, to aid debugging if configure makes a mistake. -It was created by gslnls $as_me 1.3.2, which was +It was created by gslnls $as_me 1.4.0, which was generated by GNU Autoconf 2.71. Invocation command line was $ $0$ac_configure_args_raw @@ -1819,7 +1819,9 @@ struct stat; /* Most of the following tests are stolen from RCS 5.7 src/conf.sh. */ struct buf { int x; }; struct buf * (*rcsopen) (struct buf *, struct stat *, int); -static char *e (char **p, int i) +static char *e (p, i) + char **p; + int i; { return p[i]; } @@ -4063,7 +4065,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1 # report actual input values of CONFIG_FILES etc. instead of their # values after options handling. ac_log=" -This file was extended by gslnls $as_me 1.3.2, which was +This file was extended by gslnls $as_me 1.4.0, which was generated by GNU Autoconf 2.71. Invocation command line was CONFIG_FILES = $CONFIG_FILES @@ -4118,7 +4120,7 @@ ac_cs_config_escaped=`printf "%s\n" "$ac_cs_config" | sed "s/^ //; s/'/'\\\\\\\\ cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1 ac_cs_config='$ac_cs_config_escaped' ac_cs_version="\\ -gslnls config.status 1.3.2 +gslnls config.status 1.4.0 configured by $0, generated by GNU Autoconf 2.71, with options \\"\$ac_cs_config\\" diff --git a/configure.ac b/configure.ac index 4533467..5968a35 100644 --- a/configure.ac +++ b/configure.ac @@ -1,5 +1,5 @@ -AC_INIT([gslnls],[1.3.2],[https://github.com/JorisChau/gslnls/issues/]) +AC_INIT([gslnls],[1.4.0],[https://github.com/JorisChau/gslnls/issues/]) : ${R_HOME=`R RHOME`} if test -z "${R_HOME}"; then diff --git a/cran-comments.md b/cran-comments.md index 75b1cb5..3b18edb 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,14 +1,7 @@ -## CRAN package version 1.3.2 +## CRAN package version 1.4.0 * System requirements: GSL (>= 2.2) -## Comments - -* Fix errors cran check results (fedora builds) -* Removed configure.win + Makevars.win.in -* Added static Makevars.win as recommended by T. Kalibera -* Fix prototype warning in config.log - ## Test environments * ubuntu gcc R-oldrel, R-release, R-next, R-devel (rhub, gh-actions) @@ -21,4 +14,4 @@ * ubuntu-rchk R-devel (docker) * fedora gcc R-devel --use-valgrind (rhub) -* ubuntu gcc R-4.3 --use-gct (local install) +* ubuntu gcc R-4.4 --use-gct (local install) diff --git a/gslnls.Rproj b/gslnls.Rproj index 859e2ba..6c5e160 100644 --- a/gslnls.Rproj +++ b/gslnls.Rproj @@ -1,4 +1,5 @@ Version: 1.0 +ProjectId: 2083418d-d95f-4128-80b0-3f75a4f38c1a RestoreWorkspace: Default SaveWorkspace: Default @@ -17,4 +18,3 @@ StripTrailingWhitespace: Yes BuildType: Package PackageUseDevtools: Yes -PackageInstallArgs: --no-multiarch --with-keep.source --clean diff --git a/inst/unit_tests/unit_tests_gslnls.R b/inst/unit_tests/unit_tests_gslnls.R index f81c24a..1e564e8 100644 --- a/inst/unit_tests/unit_tests_gslnls.R +++ b/inst/unit_tests/unit_tests_gslnls.R @@ -48,19 +48,19 @@ for(i in 34:59) { misra1a <- nls_test_problem("Misra1a") misra1a_fit1 <- with(misra1a, gsl_nls(fn = fn, data = data, start = as.list(start), trace = TRUE, model = TRUE)) -dotest_tol("1.5.1", coef(misra1a_fit1), misra1a[["target"]]) +dotest_tol("2.1.1", coef(misra1a_fit1), misra1a[["target"]]) misra1a_fit2 <- with(misra1a, gsl_nls(fn = fn, data = data, start = start, algorithm = "dogleg", jac = TRUE, control = list(scale = "levenberg"))) -dotest_tol("1.5.2", coef(misra1a_fit2), misra1a[["target"]]) +dotest_tol("2.1.2", coef(misra1a_fit2), misra1a[["target"]]) misra1a_fit3 <- with(misra1a, gsl_nls(fn = fn, data = data, start = start, algorithm = "lmaccel", fvv = TRUE, control = list(scale = "marquardt"))) -dotest_tol("1.5.3", coef(misra1a_fit3), misra1a[["target"]]) +dotest_tol("2.1.3", coef(misra1a_fit3), misra1a[["target"]]) misra1a_fit4 <- with(misra1a, gsl_nls(fn = fn, data = data, start = start, algorithm = "lm", weights = rep(100.0, 14L), control = list(solver = "cholesky"))) -dotest_tol("1.5.4", coef(misra1a_fit4), misra1a[["target"]]) +dotest_tol("2.1.4", coef(misra1a_fit4), misra1a[["target"]]) misra1a_fit5 <- with(misra1a, gsl_nls(fn = fn, data = data, start = start, algorithm = "ddogleg", lower = c(b1 = 100, b2 = 0), control = list(solver = "svd"))) -dotest_tol("1.5.5", coef(misra1a_fit5), misra1a[["target"]]) +dotest_tol("2.1.5", coef(misra1a_fit5), misra1a[["target"]]) misra1a_fit6 <- with(misra1a, gsl_nls(fn = fn, data = data, start = start, algorithm = "subspace2D", upper = list(b1 = 500, b2 = 1), control = list(fdtype = "center"))) -dotest_tol("1.5.6", coef(misra1a_fit6), misra1a[["target"]]) +dotest_tol("2.1.6", coef(misra1a_fit6), misra1a[["target"]]) misra1a_fit7 <- with(misra1a, gsl_nls(fn = fn, data = data, start = c(b1 = 300, b2 = 0), algorithm = "lm", jac = TRUE, lower = list(b1 = 250), upper = c(b2 = 1))) -dotest_tol("1.5.7", coef(misra1a_fit7), c(b1 = 250, misra1a[["target"]]["b2"])) +dotest_tol("2.1.7", coef(misra1a_fit7), c(b1 = 250, misra1a[["target"]]["b2"])) tools::assertError(gsl_nls(fn = ~ b1 * (1 - exp(-b2 * x)), data = misra1a$data)) tools::assertError(gsl_nls(fn = misra1a$fn, data = misra1a$data, start = c(b1 = 0))) @@ -72,11 +72,11 @@ linear1[["start"]] <- c(x1 = 0, x2 = 0, x3 = 0, x4 = 0, x5 = 0) linear1[["target"]] <- c(x1 = -1, x2 = -1, x3 = -1, x4 = -1, x5 = -1) linear1_fit1 <- with(linear1, gsl_nls(fn = fn, y = y, start = start, algorithm = "lmaccel", jac = jac, fvv = TRUE, trace = TRUE)) -dotest_tol("1.5.8", coef(linear1_fit1), linear1[["target"]]) +dotest_tol("2.2.1", coef(linear1_fit1), linear1[["target"]]) linear1_fit2 <- with(linear1, gsl_nls(fn = fn, y = y, start = start, algorithm = "dogleg", weights = rep(100.0, 5L))) -dotest_tol("1.5.9", coef(linear1_fit2), linear1[["target"]]) +dotest_tol("2.2.2", coef(linear1_fit2), linear1[["target"]]) tools::assertWarning(linear1_fit3 <- with(linear1, gsl_nls(fn = fn, y = y, start = start, algorithm = "ddogleg", lower = start, upper = as.list(start)))) -dotest_tol("1.5.10", coef(linear1_fit3), linear1[["start"]]) +dotest_tol("2.2.3", coef(linear1_fit3), linear1[["start"]]) linear1_fn <- function(pars) { ## access parameters as list fn <- numeric(5L) @@ -90,21 +90,21 @@ linear1_fn <- function(pars) { ## access parameters as list } linear1_fit4 <- with(linear1, gsl_nls(fn = linear1_fn, y = y, start = as.list(start), algorithm = "lmaccel")) -dotest_tol("1.5.11", coef(linear1_fit4), linear1[["target"]]) +dotest_tol("2.3.1", coef(linear1_fit4), linear1[["target"]]) linear1_fit5 <- with(linear1, gsl_nls(fn = linear1_fn, y = y, start = as.list(start), algorithm = "subspace2D", lower = target / 2, upper = as.list(start))) -dotest_tol("1.5.12", coef(linear1_fit5), linear1[["target"]] / 2) +dotest_tol("2.3.2", coef(linear1_fit5), linear1[["target"]] / 2) ## gsl_nls_large ## formula input -misra1a_fit6 <- with(misra1a, gsl_nls_large(fn = fn, data = data, start = as.list(start), jac = TRUE, trace = TRUE)) -dotest_tol("1.6.1", coef(misra1a_fit6), misra1a[["target"]]) -misra1a_fit7 <- with(misra1a, gsl_nls_large(fn = fn, data = data, start = start, algorithm = "dogleg", jac = TRUE, control = list(scale = "levenberg"))) -dotest_tol("1.6.2", coef(misra1a_fit7), misra1a[["target"]]) -misra1a_fit8 <- with(misra1a, gsl_nls_large(fn = fn, data = data, start = start, algorithm = "lmaccel", jac = TRUE, fvv = TRUE, control = list(scale = "marquardt"))) -dotest_tol("1.6.3", coef(misra1a_fit8), misra1a[["target"]]) -misra1a_fit9 <- with(misra1a, gsl_nls_large(fn = fn, data = data, start = start, algorithm = "lm", weights = rep(1.0, 14L), jac = TRUE)) -dotest_tol("1.6.4", coef(misra1a_fit9), misra1a[["target"]]) +misra1a_fit8 <- with(misra1a, gsl_nls_large(fn = fn, data = data, start = as.list(start), jac = TRUE, trace = TRUE)) +dotest_tol("3.1.1", coef(misra1a_fit8), misra1a[["target"]]) +misra1a_fit9 <- with(misra1a, gsl_nls_large(fn = fn, data = data, start = start, algorithm = "dogleg", jac = TRUE, control = list(scale = "levenberg"))) +dotest_tol("3.1.2", coef(misra1a_fit9), misra1a[["target"]]) +misra1a_fit10 <- with(misra1a, gsl_nls_large(fn = fn, data = data, start = start, algorithm = "lmaccel", jac = TRUE, fvv = TRUE, control = list(scale = "marquardt"))) +dotest_tol("3.1.3", coef(misra1a_fit10), misra1a[["target"]]) +misra1a_fit11 <- with(misra1a, gsl_nls_large(fn = fn, data = data, start = start, algorithm = "lm", weights = rep(1.0, 14L), jac = TRUE)) +dotest_tol("3.1.4", coef(misra1a_fit11), misra1a[["target"]]) tools::assertError(gsl_nls_large(fn = ~ b1 * (1 - exp(-b2 * x)), data = misra1a$data)) tools::assertError(gsl_nls_large(fn = misra1a$fn, data = misra1a$data, start = c(b1 = 0))) @@ -112,50 +112,113 @@ tools::assertError(gsl_nls_large(fn = misra1a$fn, data = list(x = misra1a$data$x ## function input linear1_fit6 <- with(linear1, gsl_nls_large(fn = fn, y = y, start = start, jac = jac, trace = TRUE, control = list(solver = "cholesky"))) -dotest_tol("1.6.5", coef(linear1_fit6), linear1[["target"]]) +dotest_tol("3.2.1", coef(linear1_fit6), linear1[["target"]]) linear1_fit7 <- with(linear1, gsl_nls_large(fn = fn, y = y, start = start, jac = jac, algorithm = "subspace2D", control = list(solver = "svd"))) -dotest_tol("1.6.6", coef(linear1_fit7), linear1[["target"]]) +dotest_tol("3.2.2", coef(linear1_fit7), linear1[["target"]]) linear1_fit8 <- with(linear1, gsl_nls_large(fn = fn, y = y, start = start, jac = jac, algorithm = "cgst")) -dotest_tol("1.6.7", coef(linear1_fit8), linear1[["target"]]) +dotest_tol("3.2.3", coef(linear1_fit8), linear1[["target"]]) linear1_fit9 <- with(linear1, gsl_nls_large(fn = linear1_fn, y = y, start = as.list(start), algorithm = "lmaccel")) -dotest_tol("1.6.8", coef(linear1_fit9), linear1[["target"]]) +dotest_tol("3.2.4", coef(linear1_fit9), linear1[["target"]]) linear1_fit10 <- with(linear1, gsl_nls_large(fn = fn, y = y, start = start, jac = jac, weights = rep(1.0, 5L))) -dotest_tol("1.6.9", coef(linear1_fit10), linear1[["target"]]) +dotest_tol("3.2.5", coef(linear1_fit10), linear1[["target"]]) ## multistart boxbod <- nls_test_problem("BoxBOD") madsen <- nls_test_problem("Madsen example") boxbod_fit1 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = c(200, 250), b2 = c(0, 1)), trace = TRUE, control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.1", coef(boxbod_fit1), boxbod[["target"]]) +dotest_tol("4.1.1", coef(boxbod_fit1), boxbod[["target"]]) boxbod_fit2 <- with(boxbod, gsl_nls(fn = fn, data = data, start = cbind(b1 = c(200, 250), b2 = c(0, 1)), algorithm = "lmaccel", fvv = TRUE, control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.2", coef(boxbod_fit2), boxbod[["target"]]) +dotest_tol("4.1.2", coef(boxbod_fit2), boxbod[["target"]]) boxbod_fit3 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = c(200, 250), b2 = 1), weights = rep(10.0, 6L), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.3", coef(boxbod_fit3), boxbod[["target"]]) +dotest_tol("4.1.3", coef(boxbod_fit3), boxbod[["target"]]) boxbod_fit4 <- with(boxbod, gsl_nls(fn = fn, data = data, start = c(b1 = 200, b2 = NA), lower = c(b2 = 0), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.4", coef(boxbod_fit4), boxbod[["target"]]) +dotest_tol("4.1.4", coef(boxbod_fit4), boxbod[["target"]]) boxbod_fit5 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = NA, b2 = 0.5), jac = TRUE, upper = c(b2 = 1), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.5", coef(boxbod_fit5), boxbod[["target"]]) +dotest_tol("4.1.5", coef(boxbod_fit5), boxbod[["target"]]) boxbod_fit6 <- with(boxbod, gsl_nls(fn = fn, data = data, start = cbind(b1 = NA, b2 = c(0, 1)), lower = 0, upper = 250, control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.6", coef(boxbod_fit6), boxbod[["target"]]) +dotest_tol("4.1.6", coef(boxbod_fit6), boxbod[["target"]]) boxbod_fit7 <- with(boxbod, gsl_nls(fn = fn, data = data, start = list(b1 = c(200, 250), b2 = NA), jac = TRUE, lower = list(b2 = 0), upper = c(b2 = 1), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.7", coef(boxbod_fit7), boxbod[["target"]]) +dotest_tol("4.1.7", coef(boxbod_fit7), boxbod[["target"]]) madsen_fit1 <- with(madsen, gsl_nls(fn = fn, y = y, start = cbind(x1 = c(-1, 1), x2 = c(0, 1)), trace = TRUE, control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.8", coef(madsen_fit1), madsen[["target"]]) +dotest_tol("4.2.1", coef(madsen_fit1), madsen[["target"]]) madsen_fit2 <- with(madsen, gsl_nls(fn = fn, y = y, start = c(x1 = 0, x2 = NA), jac = jac, lower = c(x2 = 0), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.9", coef(madsen_fit2), madsen[["target"]]) +dotest_tol("4.2.2", coef(madsen_fit2), madsen[["target"]]) madsen_fit3 <- with(madsen, gsl_nls(fn = fn, y = y, start = c(x1 = NA, x2 = 0), lower = -1, upper = 1, weights = rep(10.0, 3L), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.10", coef(madsen_fit3), madsen[["target"]]) +dotest_tol("4.2.3", coef(madsen_fit3), madsen[["target"]]) madsen_fit4 <- with(madsen, gsl_nls(fn = fn, y = y, start = c(x1 = NA, x2 = NA), jac = jac, control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.11", coef(madsen_fit4), madsen[["target"]]) +dotest_tol("4.2.4", coef(madsen_fit4), madsen[["target"]]) madsen_fit5 <- with(madsen, gsl_nls(fn = fn, y = y, start = cbind(x1 = c(-0.5, -0.5), x2 = c(1, 1)), jac = jac, lower = c(x1 = -Inf))) -dotest_tol("1.7.12", coef(madsen_fit4), madsen[["target"]]) +dotest_tol("4.2.5", coef(madsen_fit4), madsen[["target"]]) linear1_fit10 <- with(linear1, gsl_nls(fn = linear1_fn, y = y, start = list(x1 = NA, x2 = 0, x3 = 0, x4 = 0, x5 = 0), algorithm = "lmaccel", lower = list(x1 = -5), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.13", coef(linear1_fit10), linear1[["target"]]) +dotest_tol("4.3.1", coef(linear1_fit10), linear1[["target"]]) linear1_fit11 <- with(linear1, gsl_nls(fn = linear1_fn, y = y, start = list(x1 = c(-5, 0), x2 = 0, x3 = 0, x4 = 0, x5 = NA), lower = list(x1 = -5, x5 = -5), control = list(mstart_n = 5, mstart_q = 1, mstart_r = 1.1))) -dotest_tol("1.7.14", coef(linear1_fit11), linear1[["target"]]) +dotest_tol("4.3.2", coef(linear1_fit11), linear1[["target"]]) + +## robust loss functions + +misra1a_data <- misra1a$data +misra1a_data[1, "y"] <- 25 ## outlier + +## formula input +misra1a_fit12 <- with(misra1a, gsl_nls(fn = fn, data = misra1a_data, start = as.list(start), loss = "huber")) +dotest_tol("5.1.1", misra1a_fit12$convInfo$isConv && misra1a_fit12$irls$irls_conv && max(abs(1 - coef(misra1a_fit12) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit13 <- with(misra1a, gsl_nls(fn = fn, data = data, start = as.list(start), loss = gsl_nls_loss("huber", cc = 1.0))) +dotest_tol("5.1.2", misra1a_fit13$convInfo$isConv && misra1a_fit13$irls$irls_conv && max(abs(1 - coef(misra1a_fit13) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit14 <- with(misra1a, gsl_nls(fn = fn, data = misra1a_data, start = as.list(start), loss = "barron")) +dotest_tol("5.1.3", misra1a_fit14$convInfo$isConv && misra1a_fit14$irls$irls_conv && max(abs(1 - coef(misra1a_fit14) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit15 <- with(misra1a, gsl_nls(fn = fn, data = data, start = as.list(start), loss = gsl_nls_loss("barron", cc = c(1.0, 1.345)))) +dotest_tol("5.1.4", misra1a_fit15$convInfo$isConv && misra1a_fit15$irls$irls_conv && max(abs(1 - coef(misra1a_fit15) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit16 <- with(misra1a, gsl_nls(fn = fn, data = misra1a_data, start = as.list(start), loss = "bisquare")) +dotest_tol("5.1.5", misra1a_fit16$convInfo$isConv && misra1a_fit16$irls$irls_conv && max(abs(1 - coef(misra1a_fit16) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit17 <- with(misra1a, gsl_nls(fn = fn, data = data, start = as.list(start), loss = "bisquare", algorithm = "lmaccel", jac = TRUE)) +dotest_tol("5.1.6", misra1a_fit17$convInfo$isConv && misra1a_fit17$irls$irls_conv && max(abs(1 - coef(misra1a_fit17) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit18 <- with(misra1a, gsl_nls(fn = fn, data = misra1a_data, start = as.list(start), loss = "welsh")) +dotest_tol("5.1.7", misra1a_fit18$convInfo$isConv && misra1a_fit18$irls$irls_conv && max(abs(1 - coef(misra1a_fit18) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit19 <- with(misra1a, gsl_nls(fn = fn, data = data, start = as.list(start), loss = "welsh", weights = rep(10.0, 14L))) +dotest_tol("5.1.8", misra1a_fit19$convInfo$isConv && misra1a_fit19$irls$irls_conv && max(abs(1 - coef(misra1a_fit19) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit20 <- with(misra1a, gsl_nls(fn = fn, data = misra1a_data, start = as.list(start), loss = "optimal")) +dotest_tol("5.1.9", misra1a_fit20$convInfo$isConv && misra1a_fit20$irls$irls_conv && max(abs(1 - coef(misra1a_fit20) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit21 <- with(misra1a, gsl_nls(fn = fn, data = data, start = as.list(start), loss = "optimal", fvv = TRUE, trace = TRUE)) +dotest_tol("5.1.10", misra1a_fit21$convInfo$isConv && misra1a_fit21$irls$irls_conv && max(abs(1 - coef(misra1a_fit21) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit22 <- with(misra1a, gsl_nls(fn = fn, data = misra1a_data, start = as.list(start), loss = "hampel")) +dotest_tol("5.1.11", misra1a_fit22$convInfo$isConv && misra1a_fit22$irls$irls_conv && max(abs(1 - coef(misra1a_fit22) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit23 <- with(misra1a, gsl_nls(fn = fn, data = data, start = c(b1 = NA, b2 = 1e-4), loss = "hampel")) +dotest_tol("5.1.12", misra1a_fit23$convInfo$isConv && misra1a_fit23$irls$irls_conv && max(abs(1 - coef(misra1a_fit23) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit24 <- with(misra1a, gsl_nls(fn = fn, data = misra1a_data, start = as.list(start), loss = "ggw")) +dotest_tol("5.1.13", misra1a_fit24$convInfo$isConv && misra1a_fit24$irls$irls_conv && max(abs(1 - coef(misra1a_fit24) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit25 <- with(misra1a, gsl_nls(fn = fn, data = data, start = c(b1 = NA, b2 = NA), loss = "ggw")) +dotest_tol("5.1.14", misra1a_fit25$convInfo$isConv && misra1a_fit25$irls$irls_conv && max(abs(1 - coef(misra1a_fit25) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit26 <- with(misra1a, gsl_nls(fn = fn, data = misra1a_data, start = as.list(start), loss = "lqq")) +dotest_tol("5.1.15", misra1a_fit26$convInfo$isConv && misra1a_fit26$irls$irls_conv && max(abs(1 - coef(misra1a_fit26) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit27 <- with(misra1a, gsl_nls(fn = fn, data = data, start = as.list(start), loss = "lqq", lower = c(b1 = 0, b2 = 0), upper = c(b1 = 1000, b2 = 1))) +dotest_tol("5.1.16", misra1a_fit27$convInfo$isConv && misra1a_fit27$irls$irls_conv && max(abs(1 - coef(misra1a_fit27) / misra1a$target)) < 1e-2, TRUE) +misra1a_fit28 <- with(misra1a, gsl_nls(fn = fn, data = data, start = as.list(start), loss = gsl_nls_loss("barron", cc = c(-Inf, 1.345)), control = list(irls_xtol = 1e-20))) +dotest_tol("5.1.17", !misra1a_fit28$convInfo$isConv && !misra1a_fit28$irls$irls_conv, TRUE) + +## function input +misra1a_fn <- function(b, x) b[[1]] * (1 - exp(-b[[2]] * x)) + +misra1a_fit12_1 <- with(misra1a, gsl_nls(fn = misra1a_fn, x = misra1a_data$x, y = misra1a_data$y, start = as.list(start), loss = "huber")) +dotest_tol("5.2.1", coef(misra1a_fit12_1), coef(misra1a_fit12)) +misra1a_fit14_1 <- with(misra1a, gsl_nls(fn = misra1a_fn, x = misra1a_data$x, y = misra1a_data$y, start = as.list(start), loss = "barron")) +dotest_tol("5.2.2", coef(misra1a_fit14_1), coef(misra1a_fit14)) +misra1a_fit16_1 <- with(misra1a, gsl_nls(fn = misra1a_fn, x = misra1a_data$x, y = misra1a_data$y, start = as.list(start), loss = "bisquare")) +dotest_tol("5.2.3", coef(misra1a_fit16_1), coef(misra1a_fit16)) +misra1a_fit18_1 <- with(misra1a, gsl_nls(fn = misra1a_fn, x = misra1a_data$x, y = misra1a_data$y, start = as.list(start), loss = "welsh")) +dotest_tol("5.2.4", coef(misra1a_fit18_1), coef(misra1a_fit18)) +misra1a_fit20_1 <- with(misra1a, gsl_nls(fn = misra1a_fn, x = misra1a_data$x, y = misra1a_data$y, start = as.list(start), loss = "optimal")) +dotest_tol("5.2.5", coef(misra1a_fit20_1), coef(misra1a_fit20)) +misra1a_fit22_1 <- with(misra1a, gsl_nls(fn = misra1a_fn, x = misra1a_data$x, y = misra1a_data$y, start = as.list(start), loss = "hampel")) +dotest_tol("5.2.6", coef(misra1a_fit22_1), coef(misra1a_fit22)) +misra1a_fit24_1 <- with(misra1a, gsl_nls(fn = misra1a_fn, x = misra1a_data$x, y = misra1a_data$y, start = as.list(start), loss = "ggw")) +dotest_tol("5.2.7", coef(misra1a_fit24_1), coef(misra1a_fit24)) +misra1a_fit26_1 <- with(misra1a, gsl_nls(fn = misra1a_fn, x = misra1a_data$x, y = misra1a_data$y, start = as.list(start), loss = "lqq")) +dotest_tol("5.2.8", coef(misra1a_fit26_1), coef(misra1a_fit26)) + +tools::assertWarning(gsl_nls_loss("barron", cc = c(2.5, 1.345))) ## miscellaneous @@ -199,10 +262,10 @@ ss_fit2_large <- gsl_nls_large( algorithm = "lmaccel" ) -dotest("1.8.1", inherits(ss_fit, "gsl_nls"), TRUE) -dotest("1.8.2", inherits(ss_fit_large, "gsl_nls"), TRUE) -dotest("1.8.3", inherits(ss_fit2, "gsl_nls"), TRUE) -dotest("1.8.4", inherits(ss_fit2_large, "gsl_nls"), TRUE) +dotest("6.1.1", inherits(ss_fit, "gsl_nls"), TRUE) +dotest("6.1.2", inherits(ss_fit_large, "gsl_nls"), TRUE) +dotest("6.1.3", inherits(ss_fit2, "gsl_nls"), TRUE) +dotest("6.1.4", inherits(ss_fit2_large, "gsl_nls"), TRUE) ## sparse jacobian misra1a_fit_sp <- gsl_nls_large( @@ -219,7 +282,7 @@ misra1a_fit_sp <- gsl_nls_large( x = misra1a$data$x ) -dotest_tol("1.8.5", coef(misra1a_fit_sp), misra1a[["target"]]) +dotest_tol("6.2.1", coef(misra1a_fit_sp), misra1a[["target"]]) penalty_fit_dgC <- gsl_nls_large( fn = function(theta) c(sqrt(1e-5) * (theta - 1), sum(theta^2) - 0.25), @@ -246,51 +309,72 @@ penalty_fit_dge <- gsl_nls_large( jac = function(theta) as(rbind(diag(rep(sqrt(1e-5), length(theta))), 2 * t(theta)), Class = "unpackedMatrix") ) -dotest_tol("1.8.6", inherits(penalty_fit_dgC, "gsl_nls"), TRUE) -dotest_tol("1.8.7", inherits(penalty_fit_dgR, "gsl_nls"), TRUE) -dotest_tol("1.8.8", inherits(penalty_fit_dgT, "gsl_nls"), TRUE) -dotest_tol("1.8.9", inherits(penalty_fit_dge, "gsl_nls"), TRUE) +dotest_tol("6.2.2", inherits(penalty_fit_dgC, "gsl_nls"), TRUE) +dotest_tol("6.2.3", inherits(penalty_fit_dgR, "gsl_nls"), TRUE) +dotest_tol("6.2.4", inherits(penalty_fit_dgT, "gsl_nls"), TRUE) +dotest_tol("6.2.5", inherits(penalty_fit_dge, "gsl_nls"), TRUE) ## gsl_nls methods -dotest("1.9.1", length(fitted(misra1a_fit1)), 14L) -dotest("1.9.2", nobs(misra1a_fit1), 14L) -dotest_tol("1.9.3", deviance(misra1a_fit1), .1245514) -dotest_tol("1.9.4a", sigma(misra1a_fit1), .1018788) -dotest_tol("1.9.4b", sigma(misra1a_fit4), .1018788) -dotest("1.9.5", names(summary(misra1a_fit1)), c("formula", "residuals", "sigma", "df", "cov.unscaled", "call", - "convInfo", "control", "na.action", "coefficients", "parameters")) -dotest("1.9.6", dim(predict(misra1a_fit1, interval = "prediction")), c(14L, 3L)) -dotest("1.9.7", dim(predict(misra1a_fit1, newdata = misra1a$data, interval = "confidence")), c(14L, 3L)) -dotest("1.9.8", length(residuals(misra1a_fit1)), 14L) -dotest("1.9.9", inherits(logLik(misra1a_fit1), "logLik"), TRUE) -dotest("1.9.10", df.residual(misra1a_fit1), 12L) -dotest("1.9.11", dim(vcov(misra1a_fit1)), c(2L, 2L)) -dotest("1.9.12", names(anova(misra1a_fit1, misra1a_fit2)), c("Res.Df", "Res.Sum Sq", "Df", "Sum Sq", "F value", "Pr(>F)")) +## formula input +dotest("7.1.1", length(fitted(misra1a_fit1)), 14L) +dotest("7.1.2", nobs(misra1a_fit1), 14L) +dotest_tol("7.1.3", deviance(misra1a_fit1), .1245514) +dotest_tol("7.1.4", sigma(misra1a_fit1), .1018788) +dotest_tol("7.1.5", sigma(misra1a_fit4), 1.018788) +dotest("7.1.6", length(capture.output(summary(misra1a_fit1))), 15L) +dotest("7.1.7", names(summary(misra1a_fit1)), c("formula", "residuals", "sigma", "df", "cov.unscaled", "call", + "convInfo", "control", "na.action", "coefficients", "parameters", "irls")) +dotest("7.1.8", length(capture.output(summary(misra1a_fit19, correlation = TRUE))), 22L) +dotest("7.1.9", names(summary(misra1a_fit19, correlation = TRUE)), c("formula", "residuals", "sigma", "df", "cov.unscaled", "call", "convInfo", "control", "na.action", + "coefficients", "parameters", "irls", "correlation", "symbolic.cor")) +dotest("7.1.10", length(predict(misra1a_fit1)), 14L) +dotest("7.1.11", dim(predict(misra1a_fit4, interval = "prediction")), c(14L, 3L)) +dotest("7.1.12", dim(predict(misra1a_fit12, newdata = misra1a$data, interval = "confidence")), c(14L, 3L)) +dotest("7.1.13", length(residuals(misra1a_fit1)), 14L) +dotest("7.1.14", inherits(logLik(misra1a_fit1), "logLik"), TRUE) +dotest("7.1.15", df.residual(misra1a_fit1), 12L) +dotest("7.1.16", dim(vcov(misra1a_fit1)), c(2L, 2L)) +dotest("7.1.17", names(anova(misra1a_fit1, misra1a_fit2)), c("Res.Df", "Res.Sum Sq", "Df", "Sum Sq", "F value", "Pr(>F)")) if(requireNamespace("MASS")) { - dotest("1.9.13", dim(confint(misra1a_fit1, method = "profile")), c(2L, 2L)) + dotest("7.1.18", dim(confint(misra1a_fit1, method = "profile")), c(2L, 2L)) } -dotest("1.9.14", dim(confintd(misra1a_fit1, expr = "b1 + b2")), c(1L, 3L)) -dotest("1.9.15", capture.output(misra1a_fit1)[c(1, 2)], c("Nonlinear regression model", " model: y ~ b1 * (1 - exp(-b2 * x))")) -dotest("1.9.16", length(hatvalues(misra1a_fit1)), 14L) - -dotest("1.10.1", length(fitted(madsen_fit1)), 3L) -dotest("1.10.2", nobs(madsen_fit1), 3L) -dotest_tol("1.10.3", deviance(madsen_fit1), 0.7731991) -dotest_tol("1.10.4", sigma(madsen_fit1), 0.8793174) -dotest("1.10.5", names(summary(madsen_fit1)), c("formula", "residuals", "sigma", "df", "cov.unscaled", "call", - "convInfo", "control", "na.action", "coefficients", "parameters")) -dotest("1.10.6", dim(predict(madsen_fit1, interval = "prediction")), c(3L, 3L)) -dotest("1.10.7", dim(predict(madsen_fit1, newdata = data.frame(), interval = "confidence")), c(3L, 3L)) -dotest("1.10.8a", length(residuals(madsen_fit1, type = "pearson")), 3L) -dotest("1.10.8b", length(residuals(madsen_fit1, type = "response")), 3L) -dotest("1.10.9", inherits(logLik(madsen_fit1), "logLik"), TRUE) -dotest("1.10.10", df.residual(madsen_fit1), 1L) -dotest("1.10.11", dim(vcov(madsen_fit1)), c(2L, 2L)) -dotest("1.10.12", names(anova(madsen_fit1, madsen_fit2)), c("Res.Df", "Res.Sum Sq", "Df", "Sum Sq", "F value", "Pr(>F)")) -dotest("1.10.13", dim(confint(madsen_fit1, parm = c(1L, 2L))), c(2L, 2L)) -dotest("1.10.14", dim(confintd(madsen_fit1, expr = quote(x1 - x2), dtype = "numeric")), c(1L, 3L)) -dotest("1.10.15", capture.output(madsen_fit1)[c(1, 2)], c("Nonlinear regression model", " model: y ~ fn(x)")) -dotest("1.10.16", length(hatvalues(madsen_fit1)), 3L) +dotest("7.1.19", dim(confintd(misra1a_fit1, expr = "b1 + b2")), c(1L, 3L)) +dotest("7.1.20", capture.output(misra1a_fit1)[c(1, 2)], c("Nonlinear regression model", " model: y ~ b1 * (1 - exp(-b2 * x))")) +dotest("7.1.21", length(hatvalues(misra1a_fit1)), 14L) +dotest_tol("7.1.22", sigma(misra1a_fit12), .1342963) +dotest_tol("7.1.23", sigma(misra1a_fit12_1), .1342963) +dotest("7.1.24", length(cooks.distance(misra1a_fit1)), 14L) +dotest_tol("7.1.25", misra1a_fit4$m$gradient() / sqrt(rep(100.0, 14L)), misra1a_fit4$m$gradient1()) +dotest_tol("7.1.26", misra1a_fit19$m$gradient() / sqrt(rep(10.0, 14L)), misra1a_fit19$m$gradient1()) +dotest("7.1.27", length(capture.output(misra1a_fit1)), 11L) +dotest("7.1.28", length(capture.output(misra1a_fit12)), 14L) +dotest("7.1.29", length(capture.output(misra1a_fit28)), 12L) + +tools::assertError(logLik(misra1a_fit1, REML = TRUE)) + +## function input +dotest("7.2.1", length(fitted(madsen_fit1)), 3L) +dotest("7.2.2", nobs(madsen_fit1), 3L) +dotest_tol("7.2.3", deviance(madsen_fit1), 0.7731991) +dotest_tol("7.2.4", sigma(madsen_fit1), 0.8793174) +dotest("7.2.5", names(summary(madsen_fit1)), c("formula", "residuals", "sigma", "df", "cov.unscaled", "call", + "convInfo", "control", "na.action", "coefficients", "parameters", "irls")) +dotest("7.2.6", length(predict(madsen_fit1)), 3L) +dotest("7.2.7", dim(predict(madsen_fit1, interval = "prediction")), c(3L, 3L)) +dotest("7.2.8", dim(predict(madsen_fit1, newdata = data.frame(), interval = "confidence")), c(3L, 3L)) +dotest("7.2.9", length(residuals(madsen_fit1, type = "pearson")), 3L) +dotest("7.2.10", length(residuals(madsen_fit1, type = "response")), 3L) +dotest("7.2.11", inherits(logLik(madsen_fit1), "logLik"), TRUE) +dotest("7.2.12", df.residual(madsen_fit1), 1L) +dotest("7.2.13", dim(vcov(madsen_fit1)), c(2L, 2L)) +dotest("7.2.14", names(anova(madsen_fit1, madsen_fit2)), c("Res.Df", "Res.Sum Sq", "Df", "Sum Sq", "F value", "Pr(>F)")) +dotest("7.2.15", dim(confint(madsen_fit1, parm = c(1L, 2L))), c(2L, 2L)) +dotest("7.2.16", dim(confintd(madsen_fit1, expr = quote(x1 - x2), dtype = "numeric")), c(1L, 3L)) +dotest("7.2.17", capture.output(madsen_fit1)[c(1, 2)], c("Nonlinear regression model", " model: y ~ fn(x)")) +dotest("7.2.18", length(hatvalues(madsen_fit1)), 3L) +dotest("7.2.19", length(cooks.distance(madsen_fit1)), 3L) + +tools::assertError(confint(madsen_fit1, method = "profile")) cat("Completed gslnls unit tests\n") diff --git a/man/gsl_nls.Rd b/man/gsl_nls.Rd index 1922120..561b800 100644 --- a/man/gsl_nls.Rd +++ b/man/gsl_nls.Rd @@ -1,5 +1,3 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/nls.R \name{gsl_nls} \alias{gsl_nls} \alias{gsl_nls.formula} @@ -28,7 +26,7 @@ gsl_nls(fn, ...) ... ) -\method{gsl_nls}{`function`}( +\method{gsl_nls}{function}( fn, y, start, @@ -50,12 +48,12 @@ gsl_nls(fn, ...) or as a \link{function} returning a numeric vector, with first argument the vector of parameters to be estimated. See the individual method descriptions below.} -\item{...}{additional arguments passed to the calls of \code{fn}, \code{jac} and \code{fvv} if -defined as functions.} - \item{data}{an optional data frame in which to evaluate the variables in \code{fn} if defined as a \link{formula}. Can also be a list or an environment, but not a matrix.} +\item{y}{numeric response vector if \code{fn} is defined as a \link{function}, equal in +length to the vector returned by evaluation of the function \code{fn}.} + \item{start}{a vector, list or matrix of initial parameter values or parameter ranges. \code{start} is only allowed to be missing if \code{fn} is a \code{\link{selfStart}} model. The following choices are supported: \itemize{ @@ -89,16 +87,18 @@ on some problems. \itemize{ \item \code{"default"} default squared loss function. \item \code{"huber"} Huber loss function. -\item \code{"barron"} Smooth family of robust loss functions from (Barron (2019)). -\item \code{"biweight"} Tukey's bisquare loss function. -\item \code{"welsh"} Welsh loss function. -\item \code{"optimal"} Optimal loss function as given by (Maronna et al. (2006)). -\item \code{"hampel"} Hampel loss function (Hampel et al. (1986)). -\item \code{"ggw"} Generalized Gauss-Weight loss function (Koller and Stahel (2011)). -\item \code{"lqq"} Linear Quadratic Quadratic loss function (Koller and Stahel (2011)). -}. +\item \code{"barron"} Barron's smooth family of loss functions. +\item \code{"biweight"} Tukey's biweight/bisquare loss function. +\item \code{"welsh"} Welsh/Leclerc loss function. +\item \code{"optimal"} Optimal loss function (Maronna et al. (2006), Section 5.9.1). +\item \code{"hampel"} Hampel loss function. +\item \code{"ggw"} Generalized Gauss-Weight loss function. +\item \code{"lqq"} Linear Quadratic Quadratic loss function. +} If a character string, the default tuning parameters as specified by \code{\link{gsl_nls_loss}} are used. -Instead, a list as returned by \code{\link{gsl_nls_loss}} with non-default tuning parmaeters is also accepted.} +Instead, a list as returned by \code{\link{gsl_nls_loss}} with non-default tuning parameters is also accepted. +For all choices other than \code{rho = "default"}, iterative reweighted least squares (IRLS) is used to +solve the MM-estimation problem.} \item{control}{an optional list of control parameters to tune the least squares iterations and multistart algorithm. See \code{\link{gsl_nls_control}} for the available control parameters and their default values.} @@ -145,8 +145,8 @@ This argument is only used if \code{fn} is defined as a \link{formula}.} \item{model}{a logical value. If \code{TRUE}, the model frame is returned as part of the object. Defaults to \code{FALSE}. This argument is only used if \code{fn} is defined as a \link{formula}.} -\item{y}{numeric response vector if \code{fn} is defined as a \link{function}, equal in -length to the vector returned by evaluation of the function \code{fn}.} +\item{...}{additional arguments passed to the calls of \code{fn}, \code{jac} and \code{fvv} if +defined as functions.} } \value{ If \code{fn} is a \code{formula} returns a list object of class \code{nls}. @@ -164,10 +164,10 @@ the GNU Scientific Library (GSL). \item \code{gsl_nls(formula)}: If \code{fn} is a \code{formula}, the returned list object is of classes \code{gsl_nls} and \code{nls}. Therefore, all generic functions applicable to objects of class \code{nls}, such as \code{anova}, \code{coef}, \code{confint}, \code{deviance}, \code{df.residual}, \code{fitted}, \code{formula}, \code{logLik}, \code{nobs}, \code{predict}, \code{print}, \code{profile}, -\code{residuals}, \code{summary}, \code{vcov} and \code{weights} are also applicable to the returned list object. +\code{residuals}, \code{summary}, \code{vcov}, \code{hatvalues}, \code{cooks.distance} and \code{weights} are also applicable to the returned list object. In addition, a method \code{confintd} is available for inference of derived parameters. -\item \code{gsl_nls(`function`)}: If \code{fn} is a \code{function}, the first argument must be the vector of parameters and +\item \code{gsl_nls(function)}: If \code{fn} is a \code{function}, the first argument must be the vector of parameters and the function should return a numeric vector containing the nonlinear model evaluations at the provided parameter and predictor or covariate vectors. In addition, the argument \code{y} needs to contain the numeric vector of observed responses, equal in length to the numeric @@ -175,8 +175,8 @@ vector returned by \code{fn}. The returned list object is (only) of class \code{ Although the returned object is not of class \code{nls}, the following generic functions remain applicable for an object of class \code{gsl_nls}: \code{anova}, \code{coef}, \code{confint}, \code{deviance}, \code{df.residual}, \code{fitted}, \code{formula}, \code{logLik}, \code{nobs}, \code{predict}, \code{print}, -\code{residuals}, \code{summary}, \code{vcov} and \code{weights}. In addition, a method \code{confintd} -is available for inference of derived parameters. +\code{residuals}, \code{summary}, \code{vcov}, \code{hatvalues}, \code{cooks.distance} and \code{weights}. +In addition, a method \code{confintd} is available for inference of derived parameters. }} \section{Multi-start algorithm}{ @@ -209,11 +209,6 @@ the values of \code{mstart_n}, \code{mstart_r} or \code{mstart_minsp} to avoid e increased computational effort. } -\section{Robust loss functions}{ - -If \code{loss} -} - \examples{ # Example 1: exponential model # (https://www.gnu.org/software/gsl/doc/html/nls.html#exponential-fitting-example) @@ -247,6 +242,14 @@ gsl_nls( start = c(A = NA, lam = NA, b = NA) ## dynamic starting ranges ) +## robust regression +gsl_nls( + fn = y ~ A * exp(-lam * x) + b, ## model formula + data = data.frame(x = x, y = y), ## model fit data + start = c(A = 0, lam = 0, b = 0), ## starting values + loss = "barron" ## L1-L2 loss +) + ## analytic Jacobian 1 gsl_nls( fn = y ~ A * exp(-lam * x) + b, ## model formula diff --git a/man/gsl_nls_control.Rd b/man/gsl_nls_control.Rd index 89b2a88..1a6f032 100644 --- a/man/gsl_nls_control.Rd +++ b/man/gsl_nls_control.Rd @@ -108,7 +108,7 @@ If the starting ranges for one or more parameters are unbounded and updated dyna \item{mstart_minsp}{positive integer, minimum number of detected stationary points before the multi-start algorithm terminates, defaults to 1.} -\item{irls_maxiter}{postivive integer, maximum number of IRLS iterations, defaults to 50. Only used in case of a non-default loss function (\code{loss != "default"}) optimized by IRLS.} +\item{irls_maxiter}{posivive integer, maximum number of IRLS iterations, defaults to 50. Only used in case of a non-default loss function (\code{loss != "default"}) optimized by IRLS.} \item{irls_xtol}{numeric value, termination of the IRLS procedure occurs when the relative change in parameters between IRLS iterations is \code{<= irls_xtol}, defaults to \code{.Machine$double.eps^(1/4)}. Only used in case of a non-default loss function (\code{loss != "default"}) optimized by IRLS.} diff --git a/man/gsl_nls_large.Rd b/man/gsl_nls_large.Rd index 84ae9a2..5a40b06 100644 --- a/man/gsl_nls_large.Rd +++ b/man/gsl_nls_large.Rd @@ -1,5 +1,3 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/nls_large.R \name{gsl_nls_large} \alias{gsl_nls_large} \alias{gsl_nls_large.formula} @@ -24,7 +22,7 @@ gsl_nls_large(fn, ...) ... ) -\method{gsl_nls_large}{`function`}( +\method{gsl_nls_large}{function}( fn, y, start, @@ -42,12 +40,12 @@ gsl_nls_large(fn, ...) or as a \link{function} returning a numeric vector, with first argument the vector of parameters to be estimated. See the individual method descriptions below.} -\item{...}{additional arguments passed to the calls of \code{fn}, \code{jac} and \code{fvv} if -defined as functions.} - \item{data}{an optional data frame in which to evaluate the variables in \code{fn} if defined as a \link{formula}. Can also be a list or an environment, but not a matrix.} +\item{y}{numeric response vector if \code{fn} is defined as a \link{function}, equal in +length to the vector returned by evaluation of the function \code{fn}.} + \item{start}{a named list or named numeric vector of starting estimates. \code{start} is only allowed to be missing if \code{fn} is a \code{\link{selfStart}} model. If \code{fn} is a \code{formula}, a naive guess for \code{start} is tried, but this should not be relied on.} @@ -110,8 +108,8 @@ This argument is only used if \code{fn} is defined as a \link{formula}.} \item{model}{a logical value. If \code{TRUE}, the model frame is returned as part of the object. Defaults to \code{FALSE}. This argument is only used if \code{fn} is defined as a \link{formula}.} -\item{y}{numeric response vector if \code{fn} is defined as a \link{function}, equal in -length to the vector returned by evaluation of the function \code{fn}.} +\item{...}{additional arguments passed to the calls of \code{fn}, \code{jac} and \code{fvv} if +defined as functions.} } \value{ If \code{fn} is a \code{formula} returns a list object of class \code{nls}. @@ -129,10 +127,10 @@ the GNU Scientific Library (GSL). \item \code{gsl_nls_large(formula)}: If \code{fn} is a \code{formula}, the returned list object is of classes \code{gsl_nls} and \code{nls}. Therefore, all generic functions applicable to objects of class \code{nls}, such as \code{anova}, \code{coef}, \code{confint}, \code{deviance}, \code{df.residual}, \code{fitted}, \code{formula}, \code{logLik}, \code{nobs}, \code{predict}, \code{print}, \code{profile}, -\code{residuals}, \code{summary}, \code{vcov} and \code{weights} are also applicable to the returned list object. -In addition, a method \code{confintd} is available for inference of derived parameters. +\code{residuals}, \code{summary}, \code{vcov}, \code{hatvalues}, \code{cooks.distance} and \code{weights} are also applicable to the returned +list object. In addition, a method \code{confintd} is available for inference of derived parameters. -\item \code{gsl_nls_large(`function`)}: If \code{fn} is a \code{function}, the first argument must be the vector of parameters and +\item \code{gsl_nls_large(function)}: If \code{fn} is a \code{function}, the first argument must be the vector of parameters and the function should return a numeric vector containing the nonlinear model evaluations at the provided parameter and predictor or covariate vectors. In addition, the argument \code{y} needs to contain the numeric vector of observed responses, equal in length to the numeric @@ -140,8 +138,8 @@ vector returned by \code{fn}. The returned list object is (only) of class \code{ Although the returned object is not of class \code{nls}, the following generic functions remain applicable for an object of class \code{gsl_nls}: \code{anova}, \code{coef}, \code{confint}, \code{deviance}, \code{df.residual}, \code{fitted}, \code{formula}, \code{logLik}, \code{nobs}, \code{predict}, \code{print}, -\code{residuals}, \code{summary}, \code{vcov} and \code{weights}. In addition, a method \code{confintd} -is available for inference of derived parameters. +\code{residuals}, \code{summary}, \code{vcov}, \code{hatvalues}, \code{cooks.distance} and \code{weights}. +In addition, a method \code{confintd} is available for inference of derived parameters. }} \examples{ diff --git a/man/gsl_nls_loss.Rd b/man/gsl_nls_loss.Rd index 57b3aad..727ea51 100644 --- a/man/gsl_nls_loss.Rd +++ b/man/gsl_nls_loss.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/nls_rho.R \name{gsl_nls_loss} \alias{gsl_nls_loss} -\title{Loss functions with tunable parameters} +\title{Robust loss functions with tunable parameters} \usage{ gsl_nls_loss( rho = c("default", "huber", "barron", "bisquare", "welsh", "optimal", "hampel", "ggw", @@ -11,9 +11,9 @@ gsl_nls_loss( ) } \arguments{ -\item{rho}{character loss function name, one of \code{"default", "huber", "barron", "absolute", "bisquare", "welsh", "optimal", "hampel", "ggw", "lqq"}.} +\item{rho}{character loss function, one of \code{"default", "huber", "barron", "bisquare", "welsh", "optimal", "hampel", "ggw", "lqq"}.} -\item{cc}{numeric vector of loss function tuning parameters, with length depends on the chosen loss function. +\item{cc}{named or unnamed numeric vector of tuning parameters. The length of this argument depends on the selected \code{rho} function (see \sQuote{Details}). If \code{NULL}, the default tuning parameters are returned.} } \value{ @@ -22,16 +22,107 @@ A \code{list} with two components: \item rho \item cc } -with meanings as explained under 'Arguments'. +with meanings as explained under \sQuote{Arguments}. } \description{ -Allow the user to tune the coefficient(s) of the loss functions -available in \code{\link{gsl_nls}} and \code{\link{gsl_nls_large}}. +Allow the user to tune the coefficient(s) of the robust loss functions +supported by \code{\link{gsl_nls}}. For all choices other than \code{rho = "default"}, the MM-estimation +problem is optimized by means of iterative reweighted least squares (IRLS). } +\section{Loss function details}{ + +\describe{ + \item{\code{default}}{ + Default squared loss, no iterative reweighted least squares (IRLS) is required in this case. + \deqn{\rho(x) = x^2} + } + \item{\code{huber}}{ + Huber loss function with scaling constant \eqn{k}, defaulting to \eqn{k = 1.345} for 95\% efficiency of the regression estimator. + \deqn{\rho(x, k) = \left\{ \begin{array}{ll} \frac{1}{2} x^2 &\quad \text{ if } |x| \leq k \\[2pt] k(|x| - \frac{k}{2}) &\quad \text{ if } |x| > k \end{array} \right. } + } + \item{\code{barron}}{ + Barron's smooth family of loss functions with robustness parameter \eqn{\alpha \leq 2} (default \eqn{\alpha = 1}) and scaling constant \eqn{k} (default \eqn{k = 1.345}). + Special cases include: (scaled) squared loss for \eqn{\alpha = 2}; L1-L2 loss for \eqn{\alpha = 1}; Cauchy loss for \eqn{\alpha = 0}; Geman-McClure loss for \eqn{\alpha = -2}; + and Welsch/Leclerc loss for \eqn{\alpha = -\infty}. See Barron (2019) for additional details. + \deqn{\rho(x, \alpha, k) = \left\{ \begin{array}{ll} + \frac{1}{2} (x / k)^2 &\quad \text{ if } \alpha = 2 \\[2pt] + \log\left(\frac{1}{2} (x / k)^2 + 1 \right) &\quad \text{ if } \alpha = 0 \\[2pt] + 1 - \exp\left( -\frac{1}{2} (x / k)^2 \right) &\quad \text{ if } \alpha = -\infty \\[2pt] + \frac{|\alpha - 2|}{\alpha} \left( \left( \frac{(x / k)^2}{|\alpha-2|} + 1 \right)^{\alpha / 2} - 1 \right) &\quad \text{ otherwise } + \end{array} \right. + } + } + \item{\code{bisquare}}{ + Tukey's biweight/bisquare loss function with scaling constant \eqn{k}, defaulting to \eqn{k = 4.685} for 95\% efficiency of the regression estimator, + (\eqn{k = 1.548} gives a breakdown point of 0.5 of the S-estimator). + \deqn{\rho(x, k) = \left\{ \begin{array}{ll} 1 - (1 - (x / k)^2)^3 &\quad \text{ if } |x| \leq k \\[2pt] 1 &\quad \text{ if } |x| > k \end{array} \right. } + } + \item{\code{welsh}}{ + Welsh/Leclerc loss function with scaling constant \eqn{k}, defaulting to \eqn{k = 2.11} for 95\% efficiency of the regression estimator, (\eqn{k = 0.577} gives a + breakdown point of 0.5 of the S-estimator). This is equivalent to the Barron loss function with robustness parameter \eqn{\alpha = -\infty}. + \deqn{\rho(x, k) = 1 - \exp\left( -\frac{1}{2} (x / k)^2 \right) } + } + \item{\code{optimal}}{ + Optimal loss function as in Section 5.9.1 of Maronna et al. (2006) with scaling constant \eqn{k}, defaulting to \eqn{k = 1.060} for 95\% efficiency of the regression estimator, + (\eqn{k = 0.405} gives a breakdown point of 0.5 of the S-estimator). + \deqn{\rho(x, k) = \int_0^x \text{sign}(u) \left( - \dfrac{\phi'(|u|) + k}{\phi(|u|)} \right)^+ du} + with \eqn{\phi} the standard normal density. + } + \item{\code{hampel}}{ + Hampel loss function with a single scaling constant \eqn{k}, setting \eqn{a = 1.5k}, \eqn{b = 3.5k} and \eqn{r = 8k}. \eqn{k = 0.902} by default, + resulting in 95\% efficiency of the regression estimator, (\eqn{k = 0.212} gives a breakdown point of 0.5 of the S-estimator). + \deqn{\rho(x, a, b, r) = \left\{ \begin{array}{ll} + \frac{1}{2} x^2 / C &\quad \text{ if } |x| \leq a \\[2pt] + \left( \frac{1}{2}a^2 + a(|x|-a)\right) / C &\quad \text{ if } a < |x| \leq b \\[2pt] + \frac{a}{2}\left( 2b - a + (|x| - b) \left(1 + \frac{r - |x|}{r-b}\right) \right) / C &\quad \text{ if } b < |x| \leq r \\[2pt] + 1 &\quad \text{ if } r < |x| + \end{array} \right. } + with \eqn{C = \rho(\infty) = \frac{a}{2} (b - a + r)}. + } + \item{\code{ggw}}{ + Generalized Gauss-Weight loss function, a generalization of the Welsh/Leclerc loss with tuning constants \eqn{a, b, c}, defaulting to \eqn{a = 1.387}, + \eqn{b = 1.5}, \eqn{c = 1.063} for 95\% efficiency of the regression estimator, (\eqn{a = 0.204, b = 1.5, c = 0.296} gives a breakdown point of 0.5 of the S-estimator). + \deqn{\rho(x, a, b, c) = \int_0^x \psi(u, a, b, c) du} + with, + \deqn{\psi(x, a, b, c) = \left\{ \begin{array}{ll} + x &\quad \text{ if } |x| \leq c \\[2pt] + \exp\left( -\frac{1}{2} \frac{(|x| - c)^b}{a} \right) x &\quad \text{ if } |x| > c + \end{array} \right. } + } + \item{\code{lqq}}{ + Linear Quadratic Quadratic loss function with tuning constants \eqn{b, c, s}, defaulting to \eqn{b = 1.473}, \eqn{c = 0.982}, \eqn{s = 1.5} for + 95\% efficiency of the regression estimator, (\eqn{b = 0.402, c = 0.268, s = 1.5} gives a breakdown point of 0.5 of the S-estimator). + \deqn{\rho(x, b, c, s) = \int_0^x \psi(u, b, c, s) du} + with, + \deqn{\psi(x, b, c, s) = \left\{ \begin{array}{ll} + x &\quad \text{ if } |x| \leq c \\[2pt] + \text{sign}(x) \left( |x| - \frac{s}{2b} (|x| - c)^2 \right) &\quad \text{ if } c < |x| \leq b + c \\[2pt] + \text{sign}(x) \left( c + b - \frac{bs}{2} + \frac{s-1}{a} \left( \frac{1}{2} \tilde{x}^2 - a\tilde{x}\right) \right) &\quad \text{ if } b + c < |x| \leq a + b + c \\[2pt] + 0 &\quad \text{ otherwise } + \end{array} \right. } + where \eqn{\tilde{x} = |x| - b - c} and \eqn{a = (2c + 2b - bs) / (s - 1)}. + } +} +} + \examples{ -## huber loss function with default tuning parameter -huber_loss() +## Huber loss with default scale k +gsl_nls_loss(rho = "huber") + +## L1-L2 loss (alpha = 1) +gsl_nls_loss(rho = "barron", cc = c(1, 1.345)) + +## Cauchy loss (alpha = 0) +gsl_nls_loss(rho = "barron", cc = c(k = 1.345, alpha = 0)) + +} +\references{ +J.T. Barron (2019). \emph{A general and adaptive robust loss function}. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition (pp. 4331-4339). + +M. Galassi et al., \emph{GNU Scientific Library Reference Manual (3rd Ed.)}, ISBN 0954612078. + +R.A. Maronna et al., \emph{Robust Statistics: Theory and Methods}, ISBN 0470010924. } \seealso{ -\url{https://CRAN.R-project.org/package=robustbase/vignettes/psi_functions.pdf} +\url{https://CRAN.R-project.org/package=robustbase} } diff --git a/src/nls.c b/src/nls.c index ea338be..8abe189 100644 --- a/src/nls.c +++ b/src/nls.c @@ -448,6 +448,7 @@ SEXP C_nls_internal(void *data) status = gsl_multifit_nlinear_rho_driver(pars, &fdf, wgt_i, niter, xtol, gtol, ftol, ¶ms, &info, &chisq0, &chisq1, &irls_sigma, &irls_iter, &irls_status, verbose); + /* irls x tolerance */ for (R_len_t k = 0; k < p; k++) { @@ -455,6 +456,7 @@ SEXP C_nls_internal(void *data) double x1 = gsl_vector_get((pars->w)->x, k); irls_delta = gsl_max(irls_delta, fabs(x0 - x1)); } + } R_len_t iter = (R_len_t)gsl_multifit_nlinear_niter(pars->w); @@ -716,7 +718,7 @@ static int gsl_f(const gsl_vector *x, void *params, gsl_vector *f) if (R_IsNaN(fvalptr[i]) || !R_finite(fvalptr[i])) gsl_vector_set(f, i, (double)GSL_POSINF); else - gsl_vector_set(f, i, yptr[i] - fvalptr[i]); + gsl_vector_set(f, i, fvalptr[i] - yptr[i]); } UNPROTECT(2); diff --git a/src/nls_irls.c b/src/nls_irls.c index 6a22ec0..7047584 100644 --- a/src/nls_irls.c +++ b/src/nls_irls.c @@ -149,12 +149,15 @@ static double psip_opt(double x, double c) } -static double psi_hmpl(double x, double *k) +static double psi_hmpl(double x, double k) { /* * psi = rho' : First derivative of Hampel's loss function - * constants (a, b, r) == k[0:2] s.t. slope of psi is 1 in the center + * constants (a, b, r) s.t. slope of psi is 1 in the center */ + double a = 1.5 * k; + double b = 3.5 * k; + double r = 8.0 * k; double sx, u; if (x < 0) { @@ -166,30 +169,33 @@ static double psi_hmpl(double x, double *k) sx = +1; u = x; } - if (u <= k[0]) + if (u <= a) return (x); - else if (u <= k[1]) - return sx * k[0]; - else if (u <= k[2]) - return sx * k[0] * (k[2] - u) / (k[2] - k[1]); + else if (u <= b) + return sx * a; + else if (u <= r) + return sx * a * (r - u) / (r - b); else return 0.; } -static double psip_hmpl(double x, double *k) +static double psip_hmpl(double x, double k) { /* * psi' = rho'' : Second derivative of Hampel's loss function - * constants (a, b, r) == k[0:2] s.t. slope of psi is 1 in the center + * constants (a, b, r) s.t. slope of psi is 1 in the center */ + double a = 1.5 * k; + double b = 3.5 * k; + double r = 8.0 * k; double u = fabs(x); - if (u <= k[0]) + if (u <= a) return (1); - else if (u <= k[1]) + else if (u <= b) return (0); - else if (u <= k[2]) - return (k[0] / (k[1] - k[2])); + else if (u <= r) + return (a / (b - r)); else return (0); } @@ -299,7 +305,7 @@ double psi(double x, double *c, int i) case 5: return (psi_opt(x, c[0])); // Optimal case 6: - return (psi_hmpl(x, c)); // Hampel + return (psi_hmpl(x, c[0])); // Hampel case 7: return (psi_ggw(x, c)); // GGW case 8: @@ -326,7 +332,7 @@ double psip(double x, double *c, int i) case 5: return (psip_opt(x, c[0])); // Optimal case 6: - return (psip_hmpl(x, c)); // Hampel + return (psip_hmpl(x, c[0])); // Hampel case 7: return (psip_ggw(x, c)); // GGW case 8: @@ -466,32 +472,19 @@ int gsl_multifit_nlinear_rho_driver( */ if (status == GSL_EBADFUNC || (status == GSL_ENOPROG && *irls_iter == 1)) { + UNPROTECT(1); *info = status; return status; } /* update weights */ - // double rmin = GSL_POSINF; - // double rmax = GSL_NEGINF; for (R_len_t i = 0; i < n; i++) { resid[i] = gsl_vector_get((pars->w)->f, i) / gsl_vector_get((pars->w)->sqrt_wts, i); abs_resid[i] = fabs(resid[i]); - // if(abs_resid[i] < rmin) - // rmin = abs_resid[i]; - // if(abs_resid[i] > rmax) - // rmax = abs_resid[i]; } *irls_sigma = 1.48258 * gsl_stats_median(abs_resid, 1, n); - /* dynamic updates (Sigl 2015)*/ - // if (wgt_i == 3) - // { - // rmin /= *irls_sigma; - // rmax /= *irls_sigma; - // *wgt_cc_ptr = GSL_MIN(GSL_MIN(GSL_MAX(rmin, GSL_SQRT_DBL_EPSILON), *wgt_cc_ptr), rmax); - // } - double sum_wts = 0.0; for (R_len_t i = 0; i < n; i++) { @@ -518,6 +511,7 @@ int gsl_multifit_nlinear_rho_driver( if(*irls_status == GSL_SUCCESS) { + UNPROTECT(1); *info = status; return status; } @@ -532,8 +526,6 @@ int gsl_multifit_nlinear_rho_driver( } while (*irls_status == GSL_CONTINUE && *irls_iter < irls_maxiter); - UNPROTECT(1); - /* check if max iterations reached */ if (*irls_iter >= irls_maxiter && *irls_status != GSL_SUCCESS) { @@ -542,5 +534,6 @@ int gsl_multifit_nlinear_rho_driver( status = GSL_EMAXITER; } + UNPROTECT(1); return status; } /* gsl_multifit_nlinear_rho_driver() */ diff --git a/src/trust.c b/src/trust.c index 1ada355..23e77d5 100644 --- a/src/trust.c +++ b/src/trust.c @@ -259,6 +259,7 @@ int trust_iterate_lu(void *vstate, const gsl_vector *swts, status = gsl_multifit_nlinear_eval_df(x_trial, f_trial, swts, params->h_df, params->fdtype, fdf, J, state->workn); + if (status) return status;