nls2 <- function(formula, data = parent.frame(), start,
control = nls.control(),
algorithm = c("default", "plinear", "port", "brute-force",
"grid-search", "random-search", "lhs", "CPoptim",
"plinear-brute-force", "plinear-random", "plinear-lhs"),
trace = FALSE, weights, subset, ..., all = FALSE) {
qdata <- substitute(data)
# if formula is not of class "formula" convert it to "formula" class
if (!inherits(formula, "formula")) formula <- as.formula(formula,
env = parent.frame())
data.arg <- substitute(data)
subset.arg <- if (!missing(subset)) substitute(subset) else TRUE
if (!missing(subset)) {
if (is.data.frame(data)) data <- data[with(data, subset), ]
else stop("data must be a specified as a data.frame if subset specified")
}
# reconstruct arguments
L <- list(formula = formula, data = data, control = control, trace = trace)
if (!missing(start)) {
co <- try(coef(start), silent = TRUE)
if (!inherits(co, "try-error") && !is.null(co)) {
start <- co
# remove names beginning with dot (generated by plinear)
start <- start[grep("^[^.]", names(start))]
}
L$start <- start
}
finIter <- if (length(dim(L$start)) == 2) nrow(L$start) else 1
L <- append(L, list(...))
algorithm <- match.arg(algorithm)
if (algorithm == "grid-search") algorithm <- "brute-force"
if (algorithm == "CPoptim") {
L$data <- qdata
L$algorithm <- NULL
return(do.call("nls.CPoptim", L))
}
call <- match.call()
if (algorithm %in% c("brute-force", "random-search", "lhs",
"plinear-brute-force", "plinear-random", "plinear-lhs")) {
nls <- function(formula, data, start, weights, ...) {
nlsModel <- if (grepl("plinear", algorithm)) nlsModel.plinear
else nlsModel
environment(nlsModel) <- environment()
# disable nlsModel gradient error
stop <- function(...) {
msg <- "singular gradient matrix at initial parameter estimates"
if (list(...)[[1]] == msg) return()
stop(...)
}
m <- if (missing(weights)) {
nlsModel(formula, data, start)
} else {
wts <- eval(substitute(weights), data, environment(formula))
nlsModel(formula, data, start, wts)
}
structure(list(m = m,
call = call,
convInfo = list(isConv = TRUE, finIter = finIter,
finTol = NA)), class = "nls")
}
} else L$algorithm <- algorithm
if (!missing(weights)) L$weights <-
eval(substitute(weights), data, environment(formula))
# here nls may be nls, nlsModel or nlsModel.plinear
# depending on algorithm=
# nlsModel and nlsModel.plinear give the model evaluated at start
if (missing(start)) return(do.call(nls, L))
else L$start <- as.data.frame(as.list(start))
if (NROW(L$start) == 1) {
L$data <- data.arg
L$subset <- subset.arg
return(do.call(nls, L, envir = parent.frame()))
}
if (NROW(L$start) == 2) {
if (algorithm == "brute-force" || algorithm == "plinear-brute") {
rng <- as.data.frame(lapply(start, range))
mn <- rng[1,]
mx <- rng[2,]
# k is number of points in each dimension to take
# so that cross product is close to maxiter points.
k1 <- pmax(ceiling(sum(mx > mn)), 1)
k <- pmax(ceiling(control$maxiter ^ (1/k1)), 1)
DF <- as.data.frame(rbind(mn, mx, k))
names(DF) <- names(L$start)
# L$start <- expand.grid(by(DF, 1:NROW(DF), do.call, what =
#function(from, to, len) seq(from, to, length = len)))
finIter <- k^k1
L$start <- expand.grid(lapply(DF,
function(x) seq(x[1], x[2], length = x[3])))
} else if (algorithm == "lhs" || algorithm == "plinear-lhs") {
finIter <- control$maxiter
u <- t(lhs::randomLHS(finIter, NCOL(start)))
L$start <- t(u * unlist(start[1,]) + (1-u) * unlist(start[2,]))
L$start <- as.data.frame(L$start)
names(L$start) <- names(start)
} else { # if (algorithm == "random" || algorithm == "plinear-random") {
finIter <- control$maxiter
u_vec <- runif(finIter * NCOL(start))
u <- matrix(u_vec, NCOL(start))
L$start <- t(u * unlist(start[1,]) + (1-u) * unlist(start[2,]))
L$start <- as.data.frame(L$start)
names(L$start) <- names(start)
}
}
# each component of result is the result of one run
run <- function(start) {
L$start <- start
xx <- try(do.call(nls, L))
yy <- if (inherits(xx, "try-error")) NA else xx
if (trace) print(yy)
yy
}
result <- apply(L$start, 1, run)
if (all(is.na(result))) return(NULL)
# insert data argument and if !all take only result with minimum RSS
if (all) {
for(i in seq_along(result)) result[[i]]$data <- data.arg
} else {
ss <- lapply(result, function(x) if (identical(x, NA)) NA
else deviance(x)) # deviance is residual sum of squares
result <- result[[which.min(ss)]]
result$data <- data.arg
}
result
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.