The R Project SVN
R
(root)/branches/R-glmscore/src/library/stats/tests/nls.Rout.save – Rev 55359
Rev 54672 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R version 2.13.0 Under development (unstable) (2011-03-04 r54667)Copyright (C) 2011 The R Foundation for Statistical ComputingISBN 3-900051-07-0Platform: x86_64-unknown-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, http://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> # http://www.r-project.org/Licenses/>> ## tests of nls, especially of weighted fits>> .proctime00 <- proc.time()> library(stats)> options(digits=5) # to avoid trivial printed differences> options(show.signif.stars=FALSE) # avoid fancy quotes in o/p> options(show.nls.convergence=FALSE) # avoid non-diffable output> pdf("nls-test.pdf")>> ## 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+ value <- c(Asym, xmid, scal)+ names(value) <- mCall[c("Asym", "xmid", "scal")]+ value+ }> 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 modelmodel: y ~ a + b * x + c * I(x^2)data: parent.frame()a b c1.0058 0.9824 -0.0897residual sum-of-squares: 0.46Algorithm "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 modelmodel: y ~ a + b * x + c * I(x^2)data: parent.frame()a b c1.02 0.89 0.00residual sum-of-squares: 0.468Algorithm "port", convergence message: both X-convergence and relative convergence (5)> confint(fm)Waiting for profiling to be done...2.5% 97.5%a 1.00875 1.037847b 0.84138 0.914645c 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)Residuals:Min 1Q Median 3Q Max-0.01562 -0.00723 -0.00158 0.00403 0.02413Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 0.00517 0.00764 0.68 0.52x 0.99915 0.00119 841.38 <2e-16Residual standard error: 0.0132 on 7 degrees of freedomMultiple R-squared: 1, Adjusted R-squared: 1F-statistic: 7.08e+05 on 1 and 7 DF, p-value: <2e-16Correlation 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.543210.0012128 : 0.0051705 0.9991529> summary(fit, cor = TRUE)Formula: yeps ~ a + b * xParameters:Estimate Std. Error t value Pr(>|t|)a 0.00517 0.00764 0.68 0.52b 0.99915 0.00119 841.37 <2e-16Residual standard error: 0.0132 on 7 degrees of freedomCorrelation of Parameter Estimates:ab -0.89> stopifnot(all.equal(residuals(fit), residuals(fit0), 1e5,+ 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.5432101: 6.3964587: 1.34546 0.7008402: 0.00060639084: 0.00517053 0.9991533: 0.00060639084: 0.00517051 0.999153> summary(fit2, cor = TRUE)Formula: yeps ~ a + b * xParameters:Estimate Std. Error t value Pr(>|t|)a 0.00517 0.00764 0.68 0.52b 0.99915 0.00119 841.38 <2e-16Residual standard error: 0.0132 on 7 degrees of freedomCorrelation of Parameter Estimates:ab -0.89Algorithm "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, 1e-6), all.equal(cf1, cf0, 1e-6))> stopifnot(all.equal(residuals(fit2), residuals(fit0), 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-12xmid 1.4731 0.0947 15.6 8.8e-10scal 1.0385 0.0304 34.1 4.2e-14Residual 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-12xmid 1.4731 0.0947 15.6 8.8e-10scal 1.0385 0.0304 34.1 4.2e-14Residual standard error: 0.0355 on 13 degrees of freedom> stopifnot(all.equal(coef(summary(fm2)), coef(summary(fm1)), tol = 1e-6))> stopifnot(all.equal(residuals(fm2), residuals(fm1), tol = 1e-5))> stopifnot(all.equal(fitted(fm2), fitted(fm1), tol = 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 TableModel 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.01862 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-12xmid 1.4731 0.0947 15.6 8.8e-10scal 1.0385 0.0304 34.1 4.2e-14Residual standard error: 0.0355 on 13 degrees of freedom> stopifnot(all.equal(coef(summary(fm3)), coef(summary(fm1)), tol = 1e-6))> ft <- with(DNase1, density - fitted(fm3)/sqrt(wts))> stopifnot(all.equal(ft, fitted(fm1), tol = 1e-6))> # sign of residuals is reversed> r <- with(DNase1, -residuals(fm3)/sqrt(wts))> all.equal(r, residuals(fm1), tol = 1e05)[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 TableModel 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.01862 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-10scal 1.0385 0.0304 34.1 4.2e-14.lin 2.3350 0.0966 24.2 3.5e-12Residual 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)), tol = 1e-6,+ check.attributes = FALSE))> stopifnot(all.equal(residuals(fm4), residuals(fm1), tol = 1e-5))> stopifnot(all.equal(fitted(fm4), fitted(fm1), tol = 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 TableModel 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.01862 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-12xmid 1.4731 0.0947 15.6 8.8e-10scal 1.0385 0.0304 34.1 4.2e-14Residual standard error: 0.0355 on 13 degrees of freedomAlgorithm "port", convergence message: relative convergence (4)> stopifnot(all.equal(coef(summary(fm5)), coef(summary(fm1)), tol = 1e-6))> stopifnot(all.equal(residuals(fm5), residuals(fm1), tol = 1e-5))> stopifnot(all.equal(fitted(fm5), fitted(fm1), tol = 1e-6))>> ## check profiling> pfm1 <- profile(fm1)> pfm3 <- profile(fm3)> for(m in names(pfm1)) stopifnot(all.equal(pfm1[[m]], pfm3[[m]], tol=1e-5))> pfm5 <- profile(fm5)> for(m in names(pfm1))+ stopifnot(all.equal(pfm1[[m]], pfm5[[m]], tol=1e-5))> (c1 <- confint(fm1))Waiting for profiling to be done...2.5% 97.5%Asym 2.14936 2.5724xmid 1.28535 1.6966scal 0.97526 1.1068> (c4 <- confint(fm4, 1:2))Waiting for profiling to be done...2.5% 97.5%xmid 1.2866 1.6949scal 0.9757 1.1063> stopifnot(all.equal(c1[2:3, ], c4, tol = 1e-3))>> ## 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)+ print(confint(fm))+ fm <- nls(m1[[p]], start = s1[[p]], algorithm="port")+ print(fm)+ print(confint(fm))+ }Nonlinear regression modelmodel: y ~ x^bdata: parent.frame()b0.695residual sum-of-squares: 2.39Waiting for profiling to be done...2.5% 97.5%0.68704 0.70281Nonlinear regression modelmodel: y ~ x^bdata: parent.frame()b0.695residual sum-of-squares: 2.39Algorithm "port", convergence message: relative convergence (4)Waiting for profiling to be done...2.5% 97.5%0.68704 0.70281Nonlinear regression modelmodel: y2 ~ a * x^bdata: parent.frame()a b0.502 0.724residual sum-of-squares: 2.51Waiting for profiling to be done...2.5% 97.5%a 0.49494 0.50893b 0.70019 0.74767Nonlinear regression modelmodel: y2 ~ a * x^bdata: parent.frame()a b0.502 0.724residual sum-of-squares: 2.51Algorithm "port", convergence message: relative convergence (4)Waiting for profiling to be done...2.5% 97.5%a 0.49494 0.50893b 0.70019 0.74767Nonlinear regression modelmodel: y3 ~ a * (x + exp(logc))^bdata: parent.frame()a b logc0.558 0.603 -0.176residual sum-of-squares: 2.44Waiting for profiling to be done...2.5% 97.5%a 0.35006 0.66057b 0.45107 0.91473logc -0.64627 0.40946Nonlinear regression modelmodel: y3 ~ a * (x + exp(logc))^bdata: parent.frame()a b logc0.558 0.603 -0.176residual sum-of-squares: 2.44Algorithm "port", convergence message: relative convergence (4)Waiting for profiling to be done...2.5% 97.5%a 0.35006 0.66057b 0.45107 0.91473logc -0.64627 0.40946>> fm <- nls(y2~x^b, start=c(b=1), algorithm="plinear")> confint(profile(fm))2.5% 97.5%0.70019 0.74767> fm <- nls(y3 ~ (x+exp(logc))^b, start=c(b=1, logc=0), algorithm="plinear")> confint(profile(fm))2.5% 97.5%b 0.45105 0.91471logc -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))$atau par.vals.a par.vals.b1 -3.869 0.706 0.0002 -3.114 0.802 0.0003 -0.863 1.124 0.0004 0.000 1.538 0.2635 0.590 1.952 0.4466 1.070 2.423 0.5927 1.534 3.082 0.7378 1.969 4.034 0.8789 2.376 5.502 1.01410 2.751 7.929 1.14411 3.090 12.263 1.26412 3.375 20.845 1.373$btau par.vals.a par.vals.b1 -0.673 1.2087 0.02722 0.000 1.5381 0.26333 0.707 2.0026 0.49944 1.365 2.6295 0.72365 1.994 3.5762 0.95226 2.611 5.1820 1.19627 3.225 8.2162 1.46148 3.820 17.3946 1.7512attr(,"original.fit")Nonlinear regression modelmodel: y ~ gfun(a, b, x)data: parent.frame()a b1.538 0.263residual sum-of-squares: 0.389Algorithm "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.037b 0.263 0.352 0.75 0.476Residual standard error: 0.221 on 8 degrees of freedomAlgorithm "port", convergence message: relative convergence (4)attr(,"class")[1] "profile.nls" "profile"> confint(pr1)2.5% 97.5%a 0.96 5.20b 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)$atau par.vals.a par.vals.b1 -3.681 0.729 0.0002 -2.945 0.823 0.0003 -0.977 1.099 0.0004 0.000 1.500 0.243$btau par.vals.a par.vals.b1 -0.733 1.18200 0.003952 0.000 1.50000 0.242633 1.645 1.50000 0.481324 2.154 1.50000 0.578695 2.727 1.50000 0.707066 3.288 1.50000 0.85748attr(,"original.fit")Nonlinear regression modelmodel: y ~ gfun(a, b, x)data: parent.frame()a b1.500 0.243residual sum-of-squares: 0.39Algorithm "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.036b 0.243 0.356 0.68 0.514Residual standard error: 0.221 on 8 degrees of freedomAlgorithm "port", convergence message: relative convergence (4)attr(,"class")[1] "profile.nls" "profile"> confint(m2)Waiting for profiling to be done...2.5% 97.5%a 0.907 NAb NA 0.611> options(op)>> ## scoping problems> test <- function()+ {+ 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+ nls(y ~ myf(x,a,b,n), data=xy, start=c(a=1,b=1), trace=TRUE)+ }> test()8291.9 : 1 1726.02 : 0.80544 2.42971552.85 : 1.290 2.12970.431 : 1.9565 1.967026.555 : 1.9788 2.006426.503 : 1.9798 2.004626.503 : 1.9799 2.0046Nonlinear regression modelmodel: y ~ myf(x, a, b, n)data: xya b1.98 2.00residual sum-of-squares: 26.5> ## failed to find n in 2.2.x> ## found wrong n in 2.3.x> ## finally worked in 2.4.0>>> ## list 'start'> 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) %*% spParameters:Estimate Std. Error t value Pr(>|t|)k1 5.00e-02 3.07e-04 162.5 <2e-16k2 4.99e-03 5.38e-05 92.8 <2e-16sp1 1.00e+00 4.47e-03 224.0 <2e-16sp2 9.99e-01 5.00e-03 199.8 <2e-16Residual standard error: 0.00205 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) %*% spParameters:Estimate Std. Error t value Pr(>|t|)k1 5.00e-02 3.07e-04 162.5 <2e-16k2 4.99e-03 5.38e-05 92.8 <2e-16sp1 1.00e+00 4.47e-03 224.0 <2e-16sp2 9.99e-01 5.00e-03 199.8 <2e-16Residual standard error: 0.00205 on 96 degrees of freedomAlgorithm "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:> fit <- 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 message:In nls(y ~ b0[fac] + b1 * x, start = list(b0 = c(1, 1), b1 = 101), :Convergence failure: initial par violates constraints> with(fit$convInfo, ## start par. violates constraints+ stopifnot(isConv == FALSE, stopCode == 300))>> cat('Time elapsed: ', proc.time() - .proctime00,'\n')Time elapsed: 3.564 0.028 3.639 0 0>