| |
| R version 3.6.2 Patched (2020-02-12 r77795) -- "Dark and Stormy Night" |
| Copyright (C) 2020 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" |
| > |
| > ## 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) |
| > 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) |
| > if(have_MASS) print(confint(fm)) |
| 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: unsupported < 2.3.0 |
| > 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) |
| > 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 : 0.12345 0.54321 |
| 0.0012128 : 0.0051705 0.9991529 |
| > 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 |
| |
| > stopifnot(all.equal(residuals(fit), residuals(fit0), tolerance = 1e-5, |
| + check.attributes = FALSE)) |
| > stopifnot(df.residual(fit) == df.residual(fit0)) |
| > cf1 <- coef(summary(fit))[, 1:2] |
| > 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) |
| |
| > 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)) |
| > |
| > |
| > 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))) |
| + } |
| > t1 <- test(); t1$with.start |
| 8291.9 : 1 1 |
| 726.02 : 0.80544 2.42971 |
| 552.85 : 1.290 2.129 |
| 70.431 : 1.9565 1.9670 |
| 26.555 : 1.9788 2.0064 |
| 26.503 : 1.9798 2.0046 |
| 26.503 : 1.9799 2.0046 |
| 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]]))) |
| > rm(a,b) |
| > t2 <- test(FALSE) |
| > stopifnot(all.equal(lapply(t1, .n), |
| + lapply(t2, .n), tolerance = 0.16))# different random error |
| > |
| > |
| > ## 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)))) |
| > |
| > proc.time() |
| user system elapsed |
| 1.437 0.067 1.488 |