[Rd] nls profiling with algorithm="port" may violate bounds (PR#8508)
Spencer Graves
spencer.graves at pdf.com
Sat Jan 21 04:16:08 CET 2006
Hi, Ben, et al.:
The issue Ben identified with confint(nls(... )) generates a hard
failure for me. Specifically the command "confint(m1)" in his script
below under Rgui 2.2.1 first says, "Waiting for profiling to be done..."
then forces a screen to pop up with heading "R for Windows GUI
front-end" reading, "R for Windows GUI front-end has encountered a
problem and needs to close. We are sorry for the inconvenience. If you
were in the middle of something, the information you were working on
might be lost... ." When I try the same thing running R under XEmacs
with ESS, I get essentially the same response, exceppt "R for Windows
GUI" is replaced by "R for Windows terminal". In both cases, it kills
R. In both cases, sessionInfo() before this command is as follows:
> sessionInfo()
R version 2.2.1, 2005-12-20, i386-pc-mingw32
attached base packages:
[1] "stats4" "methods" "stats" "graphics" "grDevices" "utils"
[7] "datasets" "base"
I'm running Windows XP professional version 5.1 on an IBM T30
notebook computer.
Thanks to all of the R Core team for all your hard work to make R
what it is today, with these kinds of unpleasant surprises to rare.
Best Wishes,
spencer graves
bolker at zoo.ufl.edu wrote:
> [posted to R-devel, no discussion:
> resubmitting it as a bug, just so it gets
> logged appropriately]
>
> Sorry to report further difficulties with
> nls and profiling and constraints ... the problem
> this time (which I didn't check for in my last
> round of testing) is that the nls profiler doesn't
> seem to respect constraints that have been
> set when using the port algorithm.
> See test code below ...
> If I can I will try to hack the code, but I will
> probably start by redefining my function with
> some workarounds to make the fit quadratically "bad" (but well-defined)
> when the parameters are negative ...
> As always, please don't hesitate to correct me
> if I'm being an idiot ...
>
> cheers
> Ben Bolker
>
> -----------------------
> rm(list=ls())
>
> 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))
>
> try(confint(m1))
> ----------------
>
> for what it's worth, the logic appears to be OK in mle in the stats4
> library:
> --------------
> library(stats4)
>
> mfun <- function(a,b,s) {
> if (a<0 || b<0 || s<0) stop("bounds violated")
> -sum(dnorm(y,a*x/(1+a*b*x),sd=s,log=TRUE))
> }
>
> m2 = mle(minuslogl=mfun,
> start=list(a=1,b=1,s=0.1),
> method="L-BFGS-B",lower=c(0.002,0.002,0.002))
>
> confint(m2)
>
> m2b = mle(minuslogl=mfun,
> fixed=list(b=0),start=list(a=1,s=0.1),
> method="L-BFGS-B",lower=c(0.002,0.002,0.002))
> ## set boundary slightly above zero to avoid
> ## boundary cases
>
> dev <- 2*(-logLik(m2b)+logLik(m2))
> as.numeric(pchisq(dev,lower.tail=FALSE,df=1))
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
More information about the R-devel
mailing list