R Under development (unstable) (2022-03-19 r81942) -- "Unsuffered Consequences" Copyright (C) 2022 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # File src/library/stats/tests/nls.R > # Part of the R package, https://www.R-project.org > # > # This program is free software; you can redistribute it and/or modify > # it under the terms of the GNU General Public License as published by > # the Free Software Foundation; either version 2 of the License, or > # (at your option) any later version. > # > # This program is distributed in the hope that it will be useful, > # but WITHOUT ANY WARRANTY; without even the implied warranty of > # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the > # GNU General Public License for more details. > # > # A copy of the GNU General Public License is available at > # https://www.R-project.org/Licenses/ > > ## tests of nls, especially of weighted fits > > library(stats) > options(digits = 5) # to avoid trivial printed differences > options(useFancyQuotes = FALSE) # avoid fancy quotes in o/p > options(show.nls.convergence = FALSE) # avoid non-diffable output > options(warn = 1) > > have_MASS <- requireNamespace('MASS', quietly = TRUE) > > pdf("nls-test.pdf") > > ## utility for comparing nls() results: [TODO: use more often below] > .n <- function(r) r[names(r) != "call"] > > ## selfStart.default() w/ no parameters: > logist <- deriv( ~Asym/(1+exp(-(x-xmid)/scal)), c("Asym", "xmid", "scal"), + function(x, Asym, xmid, scal){} ) > logistInit <- function(mCall, LHS, data) { + xy <- sortedXyData(mCall[["x"]], LHS, data) + if(nrow(xy) < 3) stop("Too few distinct input values to fit a logistic") + Asym <- max(abs(xy[,"y"])) + if (Asym != max(xy[,"y"])) Asym <- -Asym # negative asymptote + xmid <- NLSstClosestX(xy, 0.5 * Asym) + scal <- NLSstClosestX(xy, 0.75 * Asym) - xmid + setNames(c(Asym, xmid, scal), + mCall[c("Asym", "xmid", "scal")]) + } > logist <- selfStart(logist, initial = logistInit) ##-> Error in R 1.5.0 > str(logist) function (x, Asym, xmid, scal) - attr(*, "initial")=function (mCall, LHS, data) - attr(*, "class")= chr "selfStart" > ## with parameters and getInitial(): > logist <- selfStart(logist, initial = logistInit, + parameters = c("Asym", "xmid", "scal")) > tools::assertWarning(verbose = TRUE, + in1 <- getInitial(circumference ~ logist(age, Asym, xmid, scal), Orange) + ) # no warning previously Asserted warning: selfStart initializing functions should have a final '...' argument since R 4.1.0 > ## but this then failed, now gives the same warning: > tools::assertWarning(verbose = TRUE, + fm <- nls(circumference ~ logist(age, Asym, xmid, scal), Orange) + ) Asserted warning: selfStart initializing functions should have a final '...' argument since R 4.1.0 > ## in R 4.1.{0,1,2} gave > ## Error in (attr(object, "initial"))(mCall = mCall, data = data, LHS = LHS, : > ## unused arguments (control = list(.......), trace = FALSE) > ## IGNORE_RDIFF_BEGIN > coef(summary(fm)) Estimate Std. Error t value Pr(>|t|) Asym 192.69 20.244 9.5182 7.4824e-11 xmid 728.76 107.299 6.7919 1.1202e-07 scal 353.53 81.472 4.3393 1.3382e-04 > ## IGNORE_RDIFF_END > > ## getInitial.formula() gets selfStart function from environment(formula) > plogis <- stats::SSlogis > in2 <- getInitial(circumference ~ plogis(age, Asym, xmid, scal), Orange) > ## in R <= 4.2.x gave Error: > ## no 'getInitial' method found for "function" objects > rm(plogis) > > > ## lower and upper in algorithm="port" > set.seed(123) > x <- runif(200) > a <- b <- 1; c <- -0.1 > y <- a+b*x+c*x^2+rnorm(200, sd=0.05) > plot(x,y) > curve(a+b*x+c*x^2, add = TRUE) > ## IGNORE_RDIFF_BEGIN > nls(y ~ a+b*x+c*I(x^2), start = c(a=1, b=1, c=0.1), algorithm = "port") Nonlinear regression model model: y ~ a + b * x + c * I(x^2) data: parent.frame() a b c 1.0058 0.9824 -0.0897 residual sum-of-squares: 0.46 Algorithm "port", convergence message: relative convergence (4) > (fm <- nls(y ~ a+b*x+c*I(x^2), start = c(a=1, b=1, c=0.1), + algorithm = "port", lower = c(0, 0, 0))) Nonlinear regression model model: y ~ a + b * x + c * I(x^2) data: parent.frame() a b c 1.02 0.89 0.00 residual sum-of-squares: 0.468 Algorithm "port", convergence message: both X-convergence and relative convergence (5) > ## IGNORE_RDIFF_END > if(have_MASS) { + print(confint(fm)) + } else message("skipping tests requiring the MASS package") Waiting for profiling to be done... 2.5% 97.5% a 1.00875 1.037847 b 0.84138 0.914645 c NA 0.042807 > > ## weighted nls fit > set.seed(123) > y <- x <- 1:10 > yeps <- y + rnorm(length(y), sd = 0.01) > wts <- rep(c(1, 2), length = 10); wts[5] <- 0 > fit0 <- lm(yeps ~ x, weights = wts) > ## IGNORE_RDIFF_BEGIN > summary(fit0, cor = TRUE) Call: lm(formula = yeps ~ x, weights = wts) Weighted Residuals: Min 1Q Median 3Q Max -0.01562 -0.00723 -0.00158 0.00403 0.02413 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.00517 0.00764 0.68 0.52 x 0.99915 0.00119 841.38 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0132 on 7 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 7.08e+05 on 1 and 7 DF, p-value: <2e-16 Correlation of Coefficients: (Intercept) x -0.89 > cf0 <- coef(summary(fit0))[, 1:2] > fit <- nls(yeps ~ a + b*x, start = list(a = 0.12345, b = 0.54321), + weights = wts, trace = TRUE) 112.14 (3.04e+02): par = (0.12345 0.54321) 0.0012128 (2.75e-06): par = (0.0051705 0.99915) > summary(fit, cor = TRUE) Formula: yeps ~ a + b * x Parameters: Estimate Std. Error t value Pr(>|t|) a 0.00517 0.00764 0.68 0.52 b 0.99915 0.00119 841.37 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0132 on 7 degrees of freedom Correlation of Parameter Estimates: a b -0.89 > ## IGNORE_RDIFF_END > stopifnot(all.equal(residuals(fit), residuals(fit0), tolerance = 1e-5, + check.attributes = FALSE)) > stopifnot(df.residual(fit) == df.residual(fit0)) > stopifnot(all.equal(logLik(fit), logLik(fit0), tolerance = 1e-8)) > cf1 <- coef(summary(fit))[, 1:2] > ## IGNORE_RDIFF_BEGIN > fit2 <- nls(yeps ~ a + b*x, start = list(a = 0.12345, b = 0.54321), + weights = wts, trace = TRUE, algorithm = "port") 0: 56.070572: 0.123450 0.543210 1: 6.3964587: 1.34546 0.700840 2: 0.00060639084: 0.00517053 0.999153 3: 0.00060639084: 0.00517051 0.999153 > summary(fit2, cor = TRUE) Formula: yeps ~ a + b * x Parameters: Estimate Std. Error t value Pr(>|t|) a 0.00517 0.00764 0.68 0.52 b 0.99915 0.00119 841.38 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0132 on 7 degrees of freedom Correlation of Parameter Estimates: a b -0.89 Algorithm "port", convergence message: both X-convergence and relative convergence (5) > ## IGNORE_RDIFF_END > cf2 <- coef(summary(fit2))[, 1:2] > rownames(cf0) <- c("a", "b") > # expect relative errors ca 2e-08 > stopifnot(all.equal(cf1, cf0, tolerance = 1e-6), + all.equal(cf1, cf0, tolerance = 1e-6)) > stopifnot(all.equal(residuals(fit2), residuals(fit0), tolerance = 1e5, + check.attributes = FALSE)) > stopifnot(all.equal(logLik(fit2), logLik(fit0), tolerance = 1e-8)) > > > DNase1 <- subset(DNase, Run == 1) > DNase1$wts <- rep(8:1, each = 2) > fm1 <- nls(density ~ SSlogis(log(conc), Asym, xmid, scal), + data = DNase1, weights = wts) > summary(fm1) Formula: density ~ SSlogis(log(conc), Asym, xmid, scal) Parameters: Estimate Std. Error t value Pr(>|t|) Asym 2.3350 0.0966 24.2 3.5e-12 *** xmid 1.4731 0.0947 15.6 8.8e-10 *** scal 1.0385 0.0304 34.1 4.2e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0355 on 13 degrees of freedom > > ## directly > fm2 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), + data = DNase1, weights = wts, + start = list(Asym = 3, xmid = 0, scal = 1)) > summary(fm2) Formula: density ~ Asym/(1 + exp((xmid - log(conc))/scal)) Parameters: Estimate Std. Error t value Pr(>|t|) Asym 2.3350 0.0966 24.2 3.5e-12 *** xmid 1.4731 0.0947 15.6 8.8e-10 *** scal 1.0385 0.0304 34.1 4.2e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0355 on 13 degrees of freedom > stopifnot(all.equal(coef(summary(fm2)), coef(summary(fm1)), tolerance = 1e-6)) > stopifnot(all.equal(residuals(fm2), residuals(fm1), tolerance = 1e-5)) > stopifnot(all.equal(fitted(fm2), fitted(fm1), tolerance = 1e-6)) > fm2a <- nls(density ~ Asym/(1 + exp((xmid - log(conc)))), + data = DNase1, weights = wts, + start = list(Asym = 3, xmid = 0)) > anova(fm2a, fm2) Analysis of Variance Table Model 1: density ~ Asym/(1 + exp((xmid - log(conc)))) Model 2: density ~ Asym/(1 + exp((xmid - log(conc))/scal)) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 14 0.0186 2 13 0.0164 1 0.00212 1.68 0.22 > > ## and without using weights > fm3 <- nls(~ sqrt(wts) * (density - Asym/(1 + exp((xmid - log(conc))/scal))), + data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1)) > summary(fm3) Formula: 0 ~ sqrt(wts) * (density - Asym/(1 + exp((xmid - log(conc))/scal))) Parameters: Estimate Std. Error t value Pr(>|t|) Asym 2.3350 0.0966 24.2 3.5e-12 *** xmid 1.4731 0.0947 15.6 8.8e-10 *** scal 1.0385 0.0304 34.1 4.2e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0355 on 13 degrees of freedom > stopifnot(all.equal(coef(summary(fm3)), coef(summary(fm1)), tolerance = 1e-6)) > ft <- with(DNase1, density - fitted(fm3)/sqrt(wts)) > stopifnot(all.equal(ft, fitted(fm1), tolerance = 1e-6)) > # sign of residuals is reversed > r <- with(DNase1, -residuals(fm3)/sqrt(wts)) > all.equal(r, residuals(fm1), tolerance = 1e-5) [1] TRUE > fm3a <- nls(~ sqrt(wts) * (density - Asym/(1 + exp((xmid - log(conc))))), + data = DNase1, start = list(Asym = 3, xmid = 0)) > anova(fm3a, fm3) Analysis of Variance Table Model 1: 0 ~ sqrt(wts) * (density - Asym/(1 + exp((xmid - log(conc))))) Model 2: 0 ~ sqrt(wts) * (density - Asym/(1 + exp((xmid - log(conc))/scal))) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 14 0.0186 2 13 0.0164 1 0.00212 1.68 0.22 > > ## using conditional linearity > fm4 <- nls(density ~ 1/(1 + exp((xmid - log(conc))/scal)), + data = DNase1, weights = wts, + start = list(xmid = 0, scal = 1), algorithm = "plinear") > summary(fm4) Formula: density ~ 1/(1 + exp((xmid - log(conc))/scal)) Parameters: Estimate Std. Error t value Pr(>|t|) xmid 1.4731 0.0947 15.6 8.8e-10 *** scal 1.0385 0.0304 34.1 4.2e-14 *** .lin 2.3350 0.0966 24.2 3.5e-12 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0355 on 13 degrees of freedom > cf <- coef(summary(fm4))[c(3,1,2), ] > rownames(cf)[2] <- "Asym" > stopifnot(all.equal(cf, coef(summary(fm1)), tolerance = 1e-6, + check.attributes = FALSE)) > stopifnot(all.equal(residuals(fm4), residuals(fm1), tolerance = 1e-5)) > stopifnot(all.equal(fitted(fm4), fitted(fm1), tolerance = 1e-6)) > fm4a <- nls(density ~ 1/(1 + exp((xmid - log(conc)))), + data = DNase1, weights = wts, + start = list(xmid = 0), algorithm = "plinear") > anova(fm4a, fm4) Analysis of Variance Table Model 1: density ~ 1/(1 + exp((xmid - log(conc)))) Model 2: density ~ 1/(1 + exp((xmid - log(conc))/scal)) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 14 0.0186 2 13 0.0164 1 0.00212 1.68 0.22 > > ## using 'port' > fm5 <- nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), + data = DNase1, weights = wts, + start = list(Asym = 3, xmid = 0, scal = 1), + algorithm = "port") > summary(fm5) Formula: density ~ Asym/(1 + exp((xmid - log(conc))/scal)) Parameters: Estimate Std. Error t value Pr(>|t|) Asym 2.3350 0.0966 24.2 3.5e-12 *** xmid 1.4731 0.0947 15.6 8.8e-10 *** scal 1.0385 0.0304 34.1 4.2e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.0355 on 13 degrees of freedom Algorithm "port", convergence message: relative convergence (4) > stopifnot(all.equal(coef(summary(fm5)), coef(summary(fm1)), tolerance = 1e-6)) > stopifnot(all.equal(residuals(fm5), residuals(fm1), tolerance = 1e-5)) > stopifnot(all.equal(fitted(fm5), fitted(fm1), tolerance = 1e-6)) > > ## check profiling > pfm1 <- profile(fm1) > pfm3 <- profile(fm3) > for(m in names(pfm1)) + stopifnot(all.equal(pfm1[[m]], pfm3[[m]], tolerance = 1e-5)) > pfm5 <- profile(fm5) > for(m in names(pfm1)) + stopifnot(all.equal(pfm1[[m]], pfm5[[m]], tolerance = 1e-5)) > if(have_MASS) { + print(c1 <- confint(fm1)) + print(c4 <- confint(fm4, 1:2)) + stopifnot(all.equal(c1[2:3, ], c4, tolerance = 1e-3)) + } Waiting for profiling to be done... 2.5% 97.5% Asym 2.14936 2.5724 xmid 1.28535 1.6966 scal 0.97526 1.1068 Waiting for profiling to be done... 2.5% 97.5% xmid 1.2866 1.6949 scal 0.9757 1.1063 > > ## some low-dimensional examples > npts <- 1000 > set.seed(1001) > x <- runif(npts) > b <- 0.7 > y <- x^b+rnorm(npts, sd=0.05) > a <- 0.5 > y2 <- a*x^b+rnorm(npts, sd=0.05) > c <- 1.0 > y3 <- a*(x+c)^b+rnorm(npts, sd=0.05) > d <- 0.5 > y4 <- a*(x^d+c)^b+rnorm(npts, sd=0.05) > m1 <- c(y ~ x^b, y2 ~ a*x^b, y3 ~ a*(x+exp(logc))^b) > s1 <- list(c(b=1), c(a=1,b=1), c(a=1,b=1,logc=0)) > for(p in 1:3) { + fm <- nls(m1[[p]], start = s1[[p]]) + print(fm) + if(have_MASS) print(confint(fm)) + fm <- nls(m1[[p]], start = s1[[p]], algorithm = "port") + print(fm) + if(have_MASS) print(confint(fm)) + } Nonlinear regression model model: y ~ x^b data: parent.frame() b 0.695 residual sum-of-squares: 2.39 Waiting for profiling to be done... 2.5% 97.5% 0.68704 0.70281 Nonlinear regression model model: y ~ x^b data: parent.frame() b 0.695 residual sum-of-squares: 2.39 Algorithm "port", convergence message: relative convergence (4) Waiting for profiling to be done... 2.5% 97.5% 0.68704 0.70281 Nonlinear regression model model: y2 ~ a * x^b data: parent.frame() a b 0.502 0.724 residual sum-of-squares: 2.51 Waiting for profiling to be done... 2.5% 97.5% a 0.49494 0.50893 b 0.70019 0.74767 Nonlinear regression model model: y2 ~ a * x^b data: parent.frame() a b 0.502 0.724 residual sum-of-squares: 2.51 Algorithm "port", convergence message: relative convergence (4) Waiting for profiling to be done... 2.5% 97.5% a 0.49494 0.50893 b 0.70019 0.74767 Nonlinear regression model model: y3 ~ a * (x + exp(logc))^b data: parent.frame() a b logc 0.558 0.603 -0.176 residual sum-of-squares: 2.44 Waiting for profiling to be done... 2.5% 97.5% a 0.35006 0.66057 b 0.45107 0.91473 logc -0.64627 0.40946 Nonlinear regression model model: y3 ~ a * (x + exp(logc))^b data: parent.frame() a b logc 0.558 0.603 -0.176 residual sum-of-squares: 2.44 Algorithm "port", convergence message: relative convergence (4) Waiting for profiling to be done... 2.5% 97.5% a 0.35006 0.66057 b 0.45107 0.91473 logc -0.64627 0.40946 > > if(have_MASS) { + fm <- nls(y2~x^b, start=c(b=1), algorithm="plinear") + print(confint(profile(fm))) + fm <- nls(y3 ~ (x+exp(logc))^b, start=c(b=1, logc=0), algorithm="plinear") + print(confint(profile(fm))) + } 2.5% 97.5% 0.70019 0.74767 2.5% 97.5% b 0.45105 0.91471 logc -0.64625 0.40933 > > > ## more profiling with bounds > op <- options(digits=3) > npts <- 10 > set.seed(1001) > a <- 2 > b <- 0.5 > x <- runif(npts) > y <- a*x/(1+a*b*x) + rnorm(npts, sd=0.2) > gfun <- function(a,b,x) { + if(a < 0 || b < 0) stop("bounds violated") + a*x/(1+a*b*x) + } > m1 <- nls(y ~ gfun(a,b,x), algorithm = "port", + lower = c(0,0), start = c(a=1, b=1)) > (pr1 <- profile(m1)) $a tau par.vals.a par.vals.b 1 -3.869 0.706 0.000 2 -3.114 0.802 0.000 3 -0.863 1.124 0.000 4 0.000 1.538 0.263 5 0.590 1.952 0.446 6 1.070 2.423 0.592 7 1.534 3.082 0.737 8 1.969 4.034 0.878 9 2.376 5.502 1.014 10 2.751 7.929 1.144 11 3.090 12.263 1.264 12 3.375 20.845 1.373 $b tau par.vals.a par.vals.b 1 -0.673 1.2087 0.0272 2 0.000 1.5381 0.2633 3 0.707 2.0026 0.4994 4 1.365 2.6295 0.7236 5 1.994 3.5762 0.9522 6 2.611 5.1820 1.1962 7 3.225 8.2162 1.4614 8 3.820 17.3946 1.7512 attr(,"original.fit") Nonlinear regression model model: y ~ gfun(a, b, x) data: parent.frame() a b 1.538 0.263 residual sum-of-squares: 0.389 Algorithm "port", convergence message: relative convergence (4) attr(,"summary") Formula: y ~ gfun(a, b, x) Parameters: Estimate Std. Error t value Pr(>|t|) a 1.538 0.617 2.49 0.037 * b 0.263 0.352 0.75 0.476 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.221 on 8 degrees of freedom Algorithm "port", convergence message: relative convergence (4) attr(,"class") [1] "profile.nls" "profile" > if(have_MASS) print(confint(pr1)) 2.5% 97.5% a 0.96 5.20 b NA 1.07 > > gfun <- function(a,b,x) { + if(a < 0 || b < 0 || a > 1.5 || b > 1) stop("bounds violated") + a*x/(1+a*b*x) + } > m2 <- nls(y ~ gfun(a,b,x), algorithm = "port", + lower = c(0, 0), upper=c(1.5, 1), start = c(a=1, b=1)) > profile(m2) $a tau par.vals.a par.vals.b 1 -3.681 0.729 0.000 2 -2.945 0.823 0.000 3 -0.977 1.099 0.000 4 0.000 1.500 0.243 $b tau par.vals.a par.vals.b 1 -0.733 1.18200 0.00395 2 0.000 1.50000 0.24263 3 1.645 1.50000 0.48132 4 2.154 1.50000 0.57869 5 2.727 1.50000 0.70706 6 3.288 1.50000 0.85748 attr(,"original.fit") Nonlinear regression model model: y ~ gfun(a, b, x) data: parent.frame() a b 1.500 0.243 residual sum-of-squares: 0.39 Algorithm "port", convergence message: relative convergence (4) attr(,"summary") Formula: y ~ gfun(a, b, x) Parameters: Estimate Std. Error t value Pr(>|t|) a 1.500 0.598 2.51 0.036 * b 0.243 0.356 0.68 0.514 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.221 on 8 degrees of freedom Algorithm "port", convergence message: relative convergence (4) attr(,"class") [1] "profile.nls" "profile" > if(have_MASS) print(confint(m2)) Waiting for profiling to be done... 2.5% 97.5% a 0.907 NA b NA 0.611 > options(op) > > ## scoping problems > test <- function(trace=TRUE) + { + x <- seq(0,5,len=20) + n <- 1 + y <- 2*x^2 + n + rnorm(x) + xy <- data.frame(x=x,y=y) + myf <- function(x,a,b,c) a*x^b+c + list(with.start= + nls(y ~ myf(x,a,b,n), data=xy, start=c(a=1,b=1), trace=trace), + no.start= ## cheap auto-init to 1 + suppressWarnings( + nls(y ~ myf(x,A,B,n), data=xy))) + } > ## IGNORE_RDIFF_BEGIN > t1 <- test() 8291.9 (1.06e+01): par = (1 1) 726.02 (5.04e+00): par = (0.80544 2.4297) 552.85 (4.45e+00): par = (1.29 2.129) 70.431 (1.29e+00): par = (1.9565 1.967) 26.555 (4.44e-02): par = (1.9788 2.0064) 26.503 (9.35e-05): par = (1.9798 2.0046) 26.503 (7.17e-07): par = (1.9799 2.0046) > ## IGNORE_RDIFF_END > t1$with.start Nonlinear regression model model: y ~ myf(x, a, b, n) data: xy a b 1.98 2.00 residual sum-of-squares: 26.5 > ##__with.start: > ## failed to find n in 2.2.x > ## found wrong n in 2.3.x > ## finally worked in 2.4.0 > ##__no.start: failed in 3.0.2 > ## 2018-09 fails on macOS with Accelerate framework. > stopifnot(all.equal(.n(t1[[1]]), .n(t1[[2]]), check.environment = FALSE)) > rm(a,b) > t2 <- test(FALSE) > stopifnot(all.equal(lapply(t1, .n), + lapply(t2, .n), tolerance = 0.16, # different random error + check.environment = FALSE)) > > > ## list 'start' > set.seed(101)# (remain independent of above) > getExpmat <- function(theta, t) + { + conc <- matrix(nrow = length(t), ncol = length(theta)) + for(i in 1:length(theta)) conc[, i] <- exp(-theta[i] * t) + conc + } > expsum <- as.vector(getExpmat(c(.05,.005), 1:100) %*% c(1,1)) > expsumNoisy <- expsum + max(expsum) *.001 * rnorm(100) > expsum.df <-data.frame(expsumNoisy) > > ## estimate decay rates, amplitudes with default Gauss-Newton > summary (nls(expsumNoisy ~ getExpmat(k, 1:100) %*% sp, expsum.df, + start = list(k = c(.6,.02), sp = c(1,2)))) Formula: expsumNoisy ~ getExpmat(k, 1:100) %*% sp Parameters: Estimate Std. Error t value Pr(>|t|) k1 5.00e-02 2.73e-04 183 <2e-16 *** k2 4.97e-03 4.77e-05 104 <2e-16 *** sp1 1.00e+00 3.96e-03 253 <2e-16 *** sp2 9.98e-01 4.43e-03 225 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.00182 on 96 degrees of freedom > > ## didn't work with port in 2.4.1 > summary (nls(expsumNoisy ~ getExpmat(k, 1:100) %*% sp, expsum.df, + start = list(k = c(.6,.02), sp = c(1,2)), + algorithm = "port")) Formula: expsumNoisy ~ getExpmat(k, 1:100) %*% sp Parameters: Estimate Std. Error t value Pr(>|t|) k1 5.00e-02 2.73e-04 183 <2e-16 *** k2 4.97e-03 4.77e-05 104 <2e-16 *** sp1 1.00e+00 3.96e-03 253 <2e-16 *** sp2 9.98e-01 4.43e-03 225 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.00182 on 96 degrees of freedom Algorithm "port", convergence message: both X-convergence and relative convergence (5) > > > ## PR13540 > > x <- runif(200) > b0 <- c(rep(0,100),runif(100)) > b1 <- 1 > fac <- as.factor(rep(c(0,1), each = 100)) > y <- b0 + b1*x + rnorm(200, sd=0.05) > # next failed in 2.8.1 > fit <- nls(y~b0[fac] + b1*x, start = list(b0=c(1,1), b1=1), + algorithm ="port", upper = c(100, 100, 100)) > # next did not "fail" in proposed fix: > fiB <- nls(y~b0[fac] + b1*x, start = list(b0=c(1,1), b1=101), + algorithm ="port", upper = c(100, 100, 100), + control = list(warnOnly=TRUE))# warning .. Warning in nls(y ~ b0[fac] + b1 * x, start = list(b0 = c(1, 1), b1 = 101), : Convergence failure: initial par violates constraints > with(fiB$convInfo, ## start par. violates constraints + stopifnot(isConv == FALSE, stopCode == 300)) > > > ## PR#17367 -- nls() quoting non-syntactical variable names > ## > op <- options(warn = 2)# no warnings allowed from here > ## > dN <- data.frame('NO [µmol/l]' = c(1,3,8,17), t = 1:4, check.names=FALSE) > fnN <- `NO [µmol/l]` ~ a + k* exp(t) > ## lm() works, nls() should too > lm.N <- lm(`NO [µmol/l]` ~ exp(t) , data = dN) > summary(lm.N) -> slmN > nm. <- nls(`NO [µmol/l]` ~ a + k*exp(t), start=list(a=0,k=1), data = dN) > ## In R <= 3.4.x : Error in eval(predvars, data, env) : object 'NO' not found > nmf <- nls(fnN, start=list(a=0,k=1), data = dN) > ## (ditto; gave identical error) > noC <- function(L) L[-match("call", names(L))] > stopifnot(all.equal(noC (nm.), noC (nmf))) > ## > ## with list for which as.data.frame() does not work [-> different branch, not using model.frame!] > ## list version (has been valid "forever", still doubtful, rather give error [FIXME] ?) > lsN <- c(as.list(dN), list(foo="bar")); lsN[["t"]] <- 1:8 > nmL <- nls(`NO [µmol/l]` ~ a + k*exp(t), start=list(a=0,k=1), data = lsN) > stopifnot(all.equal(coef(nmL), c(a = 5.069866, k = 0.003699669), tol = 4e-7))# seen 4.2e-8 > > ## trivial RHS -- should work even w/o 'start=' > fi1 <- nls(y ~ a, start = list(a=1)) > ## -> 2 deprecation warnings "length 1 in vector-arithmetic" from nlsModel() in R 3.4.x .. > options(op) # warnings about missing 'start' ok: > f.1 <- nls(y ~ a) # failed in R 3.4.x Warning in nls(y ~ a) : No starting values specified for some parameters. Initializing 'a' to '1.'. Consider specifying 'start' or using a selfStart model > stopifnot(all.equal(noC(f.1), noC(fi1)), + all.equal(coef(f.1), c(a = mean(y)))) > > > ##--- New option 'central' for numericDeriv() : > > ## Continuing the pnorm() example from example(numericDeriv): > > mkEnv <- function(n, from = -3, to = 3) { + stopifnot(is.numeric(n), n >= 2) + E <- new.env() + E$mean <- 0. + E$sd <- 1. + E$x <- seq(from, to, length.out = n) + E + } > > pnEnv <- mkEnv(65) # is used inside errE() : > > ## varying eps (very platform dependent?): > errE <- Vectorize(function(eps, central=FALSE) { + grad <- attr(numericDeriv(quote(pnorm(x, mean, sd)), c("mean", "sd"), + pnEnv, eps=eps, central=central), "gradient") + target <- with(pnEnv, -dnorm(x) * cbind(1, x, deparse.level=0L)) + ## return relative error {in the same sense as in all.equal()} : + sum(abs(target - grad)) / sum(abs(target)) + }) > > curve(errE(x), 1e-9, 1e-4, log="xy", n=512, ylim = c(1.5e-11, 5e-7), + xlab = quote(epsilon), ylab=quote(errE(epsilon))) -> rex > axis(1, at = 2^-(52/2), label = quote(sqrt(epsilon[c])), col=4, col.axis=4, line=-1/2) > axis(1, at = 2^-(52/3), label = quote(epsilon[c]^{1/3}), col=4, col.axis=4, line=-1/2) > curve(errE(x, central=TRUE), n=512, col=2, add = TRUE) -> rexC > ## IGNORE_RDIFF_BEGIN > str(xy1 <- approx(rex , xout= sqrt(2^-52)) ) List of 2 $ x: num 1.49e-08 $ y: num 1.56e-08 > str(xy2 <- approx(rexC, xout=(2^-52)^(1/3))) List of 2 $ x: num 6.06e-06 $ y: num 2.4e-11 > ## IGNORE_RDIFF_END > lines(xy1, type="h", col=4) > lines(xy2, type="h", col=4) > > proc.time() user system elapsed 3.481 0.135 3.613