[Rd] confint/nls
Prof Brian Ripley
ripley at stats.ox.ac.uk
Mon Jan 9 20:19:48 CET 2006
These are all solved, except those for "plinear", where you cannot profile
linear parameters and so you must specify parms (or call
confint(profile(fm)). And the third plinear model does not converge, so
isn't a useful test.
On Sat, 7 Jan 2006, Ben Bolker wrote:
>
> I have found some "issues" (bugs?) with nls confidence intervals ...
> some with the relatively new "port" algorithm, others more general
> (but possibly in the "well, don't do that" category). I have
> corresponded some with Prof. Ripley about them, but I thought I
> would just report how far I've gotten in case anyone else has
> thoughts. (I'm finding the code in stats/nls.R and stats/nle-profile.R
> quite dense & scary ...)
> All of this has been done with R-devel from 3 Jan 2006; the changes
> that Prof. Ripley already made to allow confint.nls not to crash
> when algorithm="port" are in R-devel, not R-patched.
>
> a synopsis of the problems with confint():
>
> with a 1-parameter model (is confint not appropriate for 1-parameter
> models? it doesn't say so in the docs [by the way, "normality" is
> misspelled as "nornality" in ?confint]):
>
> algorithm=default or plinear: get a complaint from qr.qty ('qr' and
> 'y' must have the same number of rows)
> port: "cannot allocate vector of size [large]" [caused by C code
> looking for dims when they aren't there]
>
> 2-parameter models:
> default OK
> port "cannot allocate vector"
> plinear "Error in xy.coords"
>
> 3-parameter models are OK
>
> I can fix the 2-parameter port case by adding drop=FALSE in
> appropriate places, but I wanted to check in just in case
> there are better/more efficient ways than my slogging through
> one case at a time ...
>
> apologies for the long message, but I am temporarily cut
> off from any way to post these files to the web.
>
> cheers
> Ben Bolker
>
> code that tests various combinations of numbers of parameters
> and algorithms:
> -----------
> resmat = array(dim=c(3,2,3),
> dimnames=list(npar=1:3,c("fit","confint"),c("default","plinear","port")))
> resmat.fix <- resmat
> ## sim. values
> 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)
>
> testfit <- function(model,start,alg) {
> tryfit <- try(fit <-
> nls(model,start=start,algorithm=alg,control=list(maxiter=200)))
> if (class(tryfit)!="try-error") {
> fitcode="OK"
> tryci <- try(confint(fit))
> if (class(tryci)!="try-error") {
> cicode="OK"
> } else cicode = as.character(tryci)
> } else {
> fitcode = as.character(tryfit)
> cicode="?"
> }
> c(fitcode,cicode)
> }
>
> m1 = c(y~x^b,y2~a*x^b,y3~a*(x+exp(logc))^b)
> m2 = c(y2~x^b,y3~(x+exp(logc))^b,y4~(x^d+exp(logc))^b)
> s1 = list(c(b=1),c(a=1,b=1),c(a=1,b=1,logc=0))
> s2 = list(c(b=1),c(b=1,logc=0),c(b=1,logc=0,d=0.5))
>
> for (p in 1:3) {
> resmat[p,,"default"] <- testfit(m1[[p]],start=s1[[p]],alg=NULL)
> resmat[p,,"port"] <- testfit(m1[[p]],start=s1[[p]],alg="port")
> }
>
> for (p in 1:3) {
> resmat[p,,"plinear"] <- testfit(m2[[p]],start=s2[[p]],alg="plinear")
> }
>
> print(resmat)
> set.seed(1002)
> example(nls,local=TRUE)
>
> diffs:
> --
> *** /usr/local/src/R/R-devel/src/library/stats/R/nls.R 2006-01-07
> 10:57:08.000000000 -0500
> --- nlsnew.R 2006-01-07 19:18:53.000000000 -0500
> ***************
> *** 266,277 ****
> gradSetArgs[[1]] <- (~attr(ans, "gradient"))[[2]]
> gradCall <-
> switch(length(gradSetArgs) - 1,
> ! call("[", gradSetArgs[[1]], gradSetArgs[[2]]),
> ! call("[", gradSetArgs[[1]], gradSetArgs[[2]],
> gradSetArgs[[2]]),
> call("[", gradSetArgs[[1]], gradSetArgs[[2]],
> gradSetArgs[[2]],
> ! gradSetArgs[[3]]),
> call("[", gradSetArgs[[1]], gradSetArgs[[2]],
> gradSetArgs[[2]],
> ! gradSetArgs[[3]], gradSetArgs[[4]]))
> getRHS.varying <- function()
> {
> ans <- getRHS.noVarying()
> --- 266,277 ----
> gradSetArgs[[1]] <- (~attr(ans, "gradient"))[[2]]
> gradCall <-
> switch(length(gradSetArgs) - 1,
> ! call("[", gradSetArgs[[1]], gradSetArgs[[2]],drop=FALSE),
> ! call("[", gradSetArgs[[1]], gradSetArgs[[2]],
> gradSetArgs[[2]],drop=FALSE),
> call("[", gradSetArgs[[1]], gradSetArgs[[2]],
> gradSetArgs[[2]],
> ! gradSetArgs[[3]],drop=FALSE),
> call("[", gradSetArgs[[1]], gradSetArgs[[2]],
> gradSetArgs[[2]],
> ! gradSetArgs[[3]], gradSetArgs[[4]],drop=FALSE))
> getRHS.varying <- function()
> {
> ans <- getRHS.noVarying()
> ***************
> *** 331,337 ****
> else {
> vary
> }, envir = thisEnv)
> ! gradCall[[length(gradCall)]] <<- useParams
> if(all(useParams)) {
> assign("setPars", setPars.noVarying, envir = thisEnv)
> assign("getPars", getPars.noVarying, envir = thisEnv)
> --- 331,337 ----
> else {
> vary
> }, envir = thisEnv)
> ! gradCall[[length(gradCall)-1]] <<- useParams
> if(all(useParams)) {
> assign("setPars", setPars.noVarying, envir = thisEnv)
> assign("getPars", getPars.noVarying, envir = thisEnv)
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
More information about the R-devel
mailing list